随机数空间喊麦制作

**粒子物理与核物理实验中的数据分析杨振伟清华大学第四讲:蒙特卡罗方法**上一讲回顾概率的基本概念随机变量与概率密度函数随机变量的平均值与方差能不通过实验对随机变量进行研究吗**本讲要点蒙特卡罗方法随机数产生子任意分布抽样之函数变换法与舍选法蒙特卡罗方法中的精度问题在粒子物理与核物理中的应用**蒙特卡罗方法简介蒙特卡罗方法就是利用一系列随机数来计算各种概率大小和随机变量均值等等的数值分析技术.
通常的步骤为:产生一系列在[0,1]之间均匀分布的随机数.
利用这些随机数按某些概率密度函数抽样生成我们感兴趣的另一随机序列.
利用这些值来估计的一些特性,例如:通过找到在区间的比例,给出积分值.
第一层面上的应用:蒙特卡罗计算=积分第二层面上的应用:蒙特卡罗变量="模拟的数据"**随机数的产生用物理方法产生真正的随机数不可重复产生速度慢用数学方法产生伪随机数可以重复产生的速度快**真随机数与伪随机数美国兰德(RAND)公司在1950年代,利用真空管中产生的噪音制作了一个含十万个真正的随机数表,并运用于其开展的所有模拟研究中.

真正的随机数与伪随机数之间的区别在于:数据串是否具有可压缩性,即能否用更短的形式来表示.

真正的随机数是不可压缩的,非常不规则,以至于无法用更短的形式来表示它.

在粒子物理与核物理研究中,随机数的可重复性经常也是非常有用的,尤其是程序的调试(debugging).

**随机数产生子目的是使在[0,1]范围内产生的伪随机数满足:均匀性;相互独立性;长周期性乘同余法友情推荐M=2K=52q+10周期=2K-2232513123010923651312342·101024251712401012**CERN库的随机数产生子PAW用户…gRandom->SetSeed();…Float_trandom=gRandom->Rndm(1);……Realrandom(1)CallRmarin(ISEED,0,0)…CallRanmar(random,1)…注意:用于产生子的随机数种子还可以用来保证后续进程的随机数不重复.
Root用户粒子物理与核物理研究中,大都采用CERN程序库提供的随机数产生子.
**随机数均匀性与相关性检验subroutinemcdoubleprecisionlamda,M,x,x0,ycallhbook1(10,'r',100,0.
,1.
,0.
)callhbook2(20,'r(i+1)vs.
r(i)',&100,0.
,1.
,100,0.
,1.
,0.
)x0=1.
lamda=1220703125!
5**13M=4294967296.
2**32doi=1,10000x=Mod(lamda*x0,M)y=x/Mcallhfill(10,real(y),0.
,1.
0)if(i.
gt.
1)call&hfill(20,real(y_old),real(y),1.
0)x0=xy_old=yenddoreturnend随机变量第I个随机变量第I+1个随机变量频数均匀性相关性**随机数均匀性与相关性检验随机变量第I个随机变量第I+1个随机变量频数均匀性相关性voidrandom(){UInt_tlambda,M,x0;TH1F*h1=newTH1F("h1","",100,0,1);TH2F*h2=newTH2F("h2","",100,0,1,100,0,1);lambda=1220703125;//5^13M=4294967296;//2^32x0=1;doubley,y_old;for(inti=0;iFill(y);if(i>1)h2->Fill(y_old,y);x0=x;y_old=y;}}**用蒙特卡罗法计算积分对于计算积分值解析解:数值解:蒙特卡罗方法:ABOx函数必须解析可积自变量不能太多对函数是否解析可积和是否太多自变量无要求在AB区间均匀投总数为N个点.
**蒙特卡罗方法中的精度问题采用蒙特卡罗方法(MC)计算积分与传统的梯形法相比有如下特点一维积分:多维积分:对于维数大于4的积分,用蒙特卡罗方计算积分总是最好.
**从均匀分布到任意分布的随机数函数变换法舍选法寻找某个函数,当函数的自变量取均匀分布值时,对应的函数值自动满足给定分布.

均匀分布给定分布从一个随机变量与对应概率密度函数最大值构成的二维均匀分布中,按概率密度函数与自变量关系曲线切割得到.

**函数变换法均匀分布任意分布**例子:指数分布抽样抽样效率为100%.
可采用函数变换法抽样的分布指数分布三维各向同性分布二维随机角度的正、余弦分布高斯分布n个自由度的2分布伽马分布二项式分布泊松分布Student分布(http://pdg.
lbl.
gov/2008/reviews/monterpp.
pdf)****舍选法问题:如何找到函数的最大值**舍选法举例subroutineacc_rejrealrvec(1)callhbook1(10,'x(r)',100,0.
,10.
,0.
)callhbook1(20,'x(r)',100,0.
,10.
,0.
)callhbook2(30,'f(x)vs.
x(r)',100,0.
,10.
,100,0.
,1.
1,0.
)fmax=-999.
doi=1,100callranmar(rvec,1)r=0+rvec(1)*(10.
-0.
)f=0.
5*exp(-r/2.
)if(fmax.
lt.
f)fmax=fenddofmax=1.
2*fmaxntot=0doi=1,10000callranmar(rvec,1)r=0+rvec(1)*(10.
-0.
)z=0.
5*exp(-r/2.
)if(z.
gt.
fmax)thenfmax=z*1.
2write(6,*)'zgreaterthanfmax'endifcallhfill(10,r,0.
,1.
0)callranmar(rvec,1)u=rvec(1)*fmaxif(u.
lt.
z)thencallhfill(20,r,0.
,1.
0)callhfill(30,r,u,1.
0)ntot=ntot+1endifenddowrite(6,*)'ntot=',ntotreturnend**舍选法举例voidacc_rej(){TH1F*h11=newTH1F("h11","",100,0,10);TH1F*h12=newTH1F("h12","",100,0,10);TH2F*h2=newTH2F("h2","",100,0,10,100,0,1);doublefmax=-999.
;for(inti=0;iUniform(0,10);doublef=0.
5*exp(-r/2.
if(fmaxUniform(0,10);doublez=0.
5*exp(-r/2.
if(z>fmax)fmax=1.
2*z;h11->Fill(r);doubleu=gRandom->Uniform(0,fmax);if(uFill(r)h2->Fill(r,u)ntot+=1;coutReset();hx=newTH1F("hx","xdis.
100,-10,10);gRandom->SetSeed();Double_tx;constDouble_tsigma=2.
0;constDouble_tmean=1.
0;constInt_tkUPDATE=1000;for(Int_ti=0;iGaus(mean,sigma);hx->Fill(x);}}产生平均值为mean标准偏差为sigma的高斯分布.
可以换为x=gRandom->Rndm(i);x=gRandom->Uniform(xup);x=gRandom->Integer(Imax);x=gRandom->Landau(mean,sigma);x=gRandom->Binomial(ntot,prob);x=gRandom->Poisson(mean);x=gRandom->PoissonD(mean);x=gRandom->Exp(tau);x=gRandom->BreitWigner(me,sig);在ROOT环境下采用已有的分布,可以容易完成布置的练习.
**蒙特卡罗统计检验例如常用来检验理论与实验符合好坏的2分布.
四个服从N(0,1)正态分布的且相互独立的随机变量平方和一定符合自由度为4的2分布思考:如果出现不符合的情况,该如何解释**Toy蒙特卡罗方法粒子物理与核物理在实验的早期设计阶段,通常利用Toy蒙特卡罗来估计可达到的测量精度(也称黑盒子方法).
AB+CCBD+EF+GHIJK末态有D,F,H,I,J,K.
研究测量E质量时实验可以达到的分辨率.
在不做探测器模拟的情况下,可以对稳定的末态粒子动量各分量进行含高斯分辨率的抽样,能损大小进行朗道分布抽样,寿命进行指数分布抽样,等等,然后在所有末态中寻找中间不稳定态E,根据能动量关系计算其对应的质量,得到的质量分布称为Toy蒙特卡罗结果.
**蒙特卡罗物理产生子目的:将理论用于某种物理过程的事例产生输出量:为对应某一物理过程的事例.
对于每个事例,给出过程产生的末态粒子和对应的动量在粒子物理与核物理实验数据分析中,为了验证某一理论或模型,常常需要理论家提供蒙特卡罗物理产生子.

**蒙特卡罗物理产生子(续)简单情形产生与粒子物理与核物理中常用的产生子程序包JETSET(PYTHIA)HERWIGARIADNEISAJETPYTHIAHERWIGKORALWEXCALIBURERATO**蒙特卡罗探测器模拟从产生子中输入粒子种类与动量,然后模拟粒子的输运过程模拟探测器响应多重散射(产生散射角)粒子衰变(产生寿命)电离能损(产生能损)电磁与强子簇射产生信号,电子学响应…输出量=模拟的数据输入重建分析软件用途:预测"物理产生子层面上的"给定假设在"探测器层面上"应该观测到的响应.

通用软件包:GEANT3(FORTRAN),GEANT4(C++)粒子与核物理中模拟的应用用于实验初期的设计阶段建模分析用于了解实验可能遇到物理过程的基本特征用于了解实验仪器自身所受到的各种影响因素与所影响的大小用于数据分析阶段的系统分析…****带电粒子在水中的输运过程模拟给定带电粒子的四动量单位厘米产生多少光子从均匀分布中产生满足一定波长分布的光子沿期伦科夫光锥方向均匀给所有光子动量每个光子开始在水中传播按光与水分子发生作用的概率抽样该光子是否被吸收或散射…**2MeV电子在水中的输运过程模拟结果显示了电子在水中发出期伦科夫光,损失能量直至被停止在水中的过程.

入射电子期伦科夫光子期伦科夫光子被水吸收2米长2米宽2米高水立方空气水**200MeV电子在水中的输运过程入射电子2米长2米宽2米高水立方空气水图中只显示能量大于1MeV的粒子原初电子在水中的轨迹电子韧致辐射产生的光子光子在水中散射发生了康普顿效应打出了电子探测器模拟(几何设置)**探测器模拟(物理过程)**这种模拟可以提供对探测器效率与预期性能的很好估计.
**CERN的蒙特卡罗模拟程序包GEANT4是模拟粒子经过物质时所发生的相互作用的一个软件包.
它的应用范围包括:空间科学医学物理粒子物理,核物理和加速器物理http://geant4.
web.
cern.
ch/geant4**蒙特卡罗方法应用举例如何确定在实验条件下,理论的概率密度函数例如:一质量为m共振宽度为的共振态在实验上观察到的概率密度函数是什么形式布莱特-魏格纳分布探测器分辨率探测效率贝叶斯定理:**应用举例(续一)也就是说,对应于真实的M,实际的M'应该是怎样一个分布如果假设**应用举例(续二)真实物理的图像在实验观测中会发生变化.
如果探测器的影响可以用函数来表达,有时积分可积.
但大多数数情况下,不能用函数表示时,蒙特卡罗方法可以给出最好的近似.