Note5: Random3

其实每天进展都还是听慢的，今天也只考虑了一个小问题，就是生成随机数。

Monte-Carlo方法产生特定分布的随机数

给定一个分布，或许不一定归一化：

P(x=x′)=f(x)P(x = x') = f(x) P(x=x′)=f(x)

方法一：

生成一个平均分布的随机数x∈[a1,a2]x \in [a_1,a_2]x∈[a1,a2]
生成一个平均分布的随机数y∈[b1,b2]y \in [b_1,b_2]y∈[b1,b2]
如果y<=P(x=x0)/(b2−b1)y<=P(x=x_0)/(b_2-b_1)y<=P(x=x0)/(b2−b1)，则保留xxx。否则，丢弃xxx。

这种方法非常直观，生成一个数，再以给定的概率决定是否保留。或者在平面内打一个点，仅落在曲线下方的点。
缺点是很浪费算力。浪费多于一半的随机数。当

方法二：

生成一个平均分布的随机数r
将它带入方程求解：

∫−∞x0f(x′)dx′=r\int_{-\infty }^{x_0} f(x')dx' = r ∫−∞x0f(x′)dx′=r

即得到要求分布的随机数。

这个方法的数学证明我还不会。很明显，节省算力，但有些分布函数在定积分方面就已经很复杂了，再计算反函数难以给出解析解。

不过好在，ROOT给出了大量的特殊函数供使用。

Section 1-2

Using random number generator TRandom3 to generate 10 random numbers andprint them out.
Define a histogram using TH1D, with 100 bins, with a range from 0 to 1. Using the random number generator to generate 1000000 random numbers and fill inside thishistogram. And then, draw the histogram.
Define a 1-D function f1, f(x) = 1=k · e−x=k, for 0 < x < 10, set parameter k = 2,and draw this function.
Define another histogram using TH1D, with 100 bins and range from 0 to 10. UsingMonte-Carlo method (with “throw away events” method), to generate the 1000000events according to the function f1 defined above.

首先又是非常简单的函数定义环节。

// define function f1 for h2 and f2 for h1
TF1* f1 = new TF1("f1","[0]*TMath::Gaus(x,[1],[2])+[3]*TMath::Gaus(x,[4],[5])",0,10);
TF1* f2 = new TF1("f1","-[0]*TMath::Log(1-x)",0,1);// give parameters
f1->SetParameters(0.3,3,0.7,0.5,5,0.7);
f2->SetParameter(0,0.5);// Draw
TCanvas* c1 = new TCanvas();
f1->DrawClone();
c1->Draw();

使用TRandom类生成随机数，首先要实例化这个类。并且设置seed。计算机伪随机数的产生算法是依据一个真随机数来产生的，因此使用不同的seed可以保证随机数的质量。

TRandm3 r(seed);

这里还可以提一句，TF1函数的值可以用Eval方法获得。

// difine a random number
TRandom3 rnd(1234);// define a histogram
TH1D* h1 = new TH1D("h1","h1",100,0,8);// fill with randoms
for (int i=0;i<100000;i++)
{h1->Fill(f2->Eval(rnd.Rndm()));
}// Draw
TCanvas * c2 = new TCanvas();
h1->Draw();
c2->Draw();

Section 2

Define a function f1, f(x) = a · G(x; µ1; σ1)+b · G(x; µ2; σ2), where G(x) is Gaussian function. Set parameters to be a = 1:3, b = 2, µ1 = 3, µ2 = 6, σ1 = 0:9, and σ2 = 1. And set the function range from 0 to 10.
Define a TRandom3 random number, and a histogram. Generate the events according to f1 using Monte-Carlo method (with “throw away events” method), and store them in the histogram. Then draw the histogram.

产生随机数之后，由于归一化的问题，直方图和PDF（概率密度分布函数）总是会在不同的scale上，把它们对比起来是有些麻烦的。这里我直接拟合（偷懒）了。

// defin a hist
TH1D* h2 = new TH1D("h2","h2",100,0,8);// fill
Double_t x = 0;
for(int i=0;i<100000;i++)
{x = 10*rnd.Rndm();if(rnd.Rndm()<f1->Eval(x)/f1->GetMaximum()){h2->Fill(x);}
}
TCanvas* c3 = new TCanvas();
h2->Fit(f1);
h2->Draw();
c3->Draw();

 FCN=107.343 FROM MIGRAD    STATUS=CONVERGED     543 CALLS         544 TOTALEDM=3.05335e-09    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   1.9 per centEXT PARAMETER                                   STEP         FIRST   NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 1  p0           7.95501e+02   8.02893e+00   3.05068e-03  -4.94309e-072  p1           5.00054e+00   1.12739e-02   6.54210e-06   9.63575e-043  p2           6.95568e-01   7.22364e-03  -2.87591e-06   1.55035e-024  p3           4.79497e+02   6.45497e+00   2.41912e-03  -7.68525e-065  p4           2.99113e+00   1.73288e-02   3.99869e-07   4.67192e-036  p5           7.05836e-01   1.01369e-02   2.87406e-06  -4.32032e-03

Section 3

Same as Exercise 2, but now, you use the Monte-Carlo transformation method to generate a Gaussian distribution for e.g. µ = 5, σ = 3.

正态分布的积分很难给出解析解，好在它的积分是一个特殊函数，叫做误差函数。ROOT中也给出了它的反函数。当然，还是有一定差别的，需要换算一下。

erf(z)=2π∫0ze−t2dterf(z) = \frac{2}{\sqrt{\pi}} \int_0^z e^{-t^2} dt erf(z)=π2∫0ze−t2dt

∫−∞xGaus(x′)dx′=−12[erf(μ−x2σ)]−∞x\int_{-\infty}^{x} Gaus(x') dx' = -\frac{1}{2} \left[erf\left(\frac{\mu-x}{\sqrt{2}\sigma }\right)\right]^{x}_{-\infty} ∫−∞xGaus(x′)dx′=−21[erf(2σμ−x)]−∞x

其他没有什么难点，都是曾经见过的。

// define a hist to store the data
TH1D* h3 = new TH1D("h3","h3",100,0,10);Double_t data = 0;  // last data for fill
for(int i=0;i<10000;i++)
{// get randoms in range of cdfdata = 5-sqrt(2)*3*TMath::ErfInverse(1-2*rnd.Rndm());h3->Fill(data);
}
TCanvas* c4 = new TCanvas();
h3->Fit("gaus");
h3->Draw();
c4->Draw();

 FCN=116.545 FROM MIGRAD    STATUS=CONVERGED      75 CALLS          76 TOTALEDM=2.18376e-08    STRATEGY= 1      ERROR MATRIX ACCURATE EXT PARAMETER                                   STEP         FIRST   NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 1  Constant     1.30797e+02   1.87669e+00   7.34889e-03   8.16277e-052  Mean         5.01755e+00   3.98001e-02   2.10686e-04  -3.17898e-033  Sigma        2.99027e+00   4.53944e-02   2.25064e-05   3.89702e-02

ROOT(a Data analysis Framework)-Note6: iSTEP day3-Random相关推荐

ROOT(a Data analysis Framework)-Note4: iSTEP day1-TH1TF1
Note4: TF1 & TH1 最近参加了一个叫做iSTEP的暑期学校,这里教授一些基本的ROOT使用方法.曾经用root的过程忘的差不多了,也没有记录.现在正好记录一下. TF1类 TF1 ...
ROOT(a Data analysis Framework)-Note5: iSTEP day2-TH1::Fit()
Note5 :数据拟合今天还是关于TH1和TF1的内容,使用Fit方法,使用不同的函数对数据进行拟合.下面以问题的形式做个小的总结记录. 摘要 TFile基本操作打开一个.root文件: TFil ...
Hi-C data analysis tools and papers
Hi-C data analysis tools and papers 全文链接如下: https://github.com/mdozmorov/HiC_tools Tools are sorted ...
python进行探索性数据分析EDA（Exploratory Data Analysis）分析
python进行探索性数据分析EDA(Exploratory Data Analysis)分析 show holy respect to python community, for there ded ...
Android11 无Root 访问data目录实现、Android11访问data目录、Android11解除data目录限制、Android11 data空白解决
Android11 无Root 访问data目录实现正文开始关于Android11权限变化作为普通安卓用户该如何方便快速地访问Android/data目录开发者该如何实现无ROOT访问Dat ...
R语言统计入门课程推荐——生物科学中的数据分析Data Analysis for the Life Sciences
Data Analysis for the Life Sciences是哈佛大学PH525x系列课程--生物医学中的数据分析(PH525x series - Biomedical Data Scien ...
R语言explore包进行探索性数据分析实战（EDA、exploratory data analysis）：基于iris数据集
R语言explore包进行探索性数据分析实战(EDA.exploratory data analysis):基于iris数据集目录
R探索新数据分析（Exploratory Data Analysis，EDA）
R探索新数据分析(Exploratory Data Analysis,EDA) 目录 R探索新数据分析(Exploratory Data Analysis,EDA) str方法进行数据概览及类型查看
ADO.NET Data Services Framework 基础概述
随着.NET Framework 4.0 及 Visual Studio 2010 的发布, ADO.NET Data Services Framework 2.0 的版本也将同时发布,新的.Net ...

ROOT(a Data analysis Framework)-Note6: iSTEP day3-Random