计算材料学

02/04/20201粒子物理与核物理实验中的数据分析杨振伟清华大学第四讲:蒙特卡罗方法02/04/20202上一讲回顾概率的基本概念随机变量与概率密度函数随机变量的平均值与方差能不通过实验对随机变量进行研究吗？02/04/20203本讲要点蒙特卡罗方法随机数产生子任意分布抽样之函数变换法与舍选法蒙特卡罗方法中的精度问题在粒子物理与核物理中的应用02/04/20204蒙特卡罗方法简介蒙特卡罗方法就是利用一系列随机数来计算各种概率大小和随机变量均值等等的数值分析技术。通常的步骤为：1)产生一系列在[0，1]之间均匀分布的随机数。2)利用这些随机数按某些概率密度函数抽样生成我们感兴趣的另一随机序列。3)利用这些值来估计的一些特性，例如：通过找到在区间的比例，给出积分值。mrrr,...,,21)(xfnxxx,...,,21xx)(xfbxabadxxf)(第一层面上的应用：蒙特卡罗计算=积分第二层面上的应用：蒙特卡罗变量=“模拟的数据”02/04/20205随机数的产生用物理方法产生真正的随机数•不可重复•产生速度慢用数学方法产生伪随机数•可以重复•产生的速度快02/04/20206真随机数与伪随机数美国兰德（RAND）公司在1950年代，利用真空管中产生的噪音制作了一个含十万个真正的随机数表，并运用于其开展的所有模拟研究中。真正的随机数与伪随机数之间的区别在于：数据串是否具有可压缩性，即能否用更短的形式来表示。真正的随机数是不可压缩的，非常不规则，以至于无法用更短的形式来表示它。在粒子物理与核物理研究中，随机数的可重复性经常也是非常有用的，尤其是程序的调试(debugging)。02/04/20207随机数产生子目的是使在[0，1]范围内产生的伪随机数满足：均匀性；相互独立性；长周期性乘同余法110(),/iiiiMrMM余数随机数与为选定常数为随机数种子。友情推荐M=2K=52q+10周期=2K-2232513123010923651312342·10102425171240101202/04/20208CERN库的随机数产生子PAW用户…gRandom-SetSeed();…Float_trandom=gRandom-Rndm(1);……Realrandom(1)CallRmarin(ISEED,0,0)…CallRanmar(random,1)…注意:用于产生子的随机数种子还可以用来保证后续进程的随机数不重复。Root用户粒子物理与核物理研究中，大都采用CERN程序库提供的随机数产生子。02/04/20209随机数均匀性与相关性检验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个随机变量频数均匀性相关性20h/pl10;h/pl2;1zonePAWmc.fcallPAW02/04/202010随机数均匀性与相关性检验随机变量第I个随机变量第I+1个随机变量频数均匀性相关性Draw()h2-root[2]Draw()h1-root[1]random.C.xroot[0]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;i10000;i++){x=(lambda*x0)%M;y=(double)x/M;h1-Fill(y);if(i1)h2-Fill(y_old,y);x0=x;y_old=y;}}02/04/202011用蒙特卡罗法计算积分对于计算积分值()bafxdx解析解：()()()bxbxaafxdxFxFx数值解：1()()nbiaibafxdxfxn蒙特卡罗方法：abABOx()fx函数必须解析可积自变量不能太多对函数是否解析可积和是否太多自变量无要求in()()()()baNfxfxdxAOBON在AB区间均匀投总数为N个点。02/04/202012蒙特卡罗方法中的精度问题采用蒙特卡罗方法(MC)计算积分与传统的梯形法相比有如下特点一维积分:)(,1/:)(,1/:MC2为子区间的数目梯形法精度为产生的随机数精度nnnn多维积分:)(,1/:)(,1/:MC/2为子区间的数目梯形法精度为产生的随机数且与维数无关精度维数nnnn对于维数大于4的积分,用蒙特卡罗方计算积分总是最好。02/04/202013从均匀分布到任意分布的随机数函数变换法舍选法寻找某个函数，当函数的自变量取均匀分布值时，对应的函数值自动满足给定分布。均匀分布给定分布从一个随机变量与对应概率密度函数最大值构成的二维均匀分布中，按概率密度函数与自变量关系曲线切割得到。02/04/202014函数变换法nxxxxfrx,...,,)()(,1][0,21分布的随机数找出服从通过适当的变换均匀分布的随机数从在')'())'((')'(')(，))'(()'(:rrxrxFdxxfrdrrgrxxPrrP即要求)()(1rFxrxF解出令均匀分布任意分布02/04/202015例子:指数分布抽样)1log()(),('1)0(1),(:)(0/'/rrxrxrdxexexfrxxx得到并解出令指数概率密度函数2呈一一对应关系与随机变量)(rxr抽样效率为100%。可采用函数变换法抽样的分布指数分布三维各向同性分布二维随机角度的正、余弦分布高斯分布n个自由度的2分布伽马分布二项式分布泊松分布Student分布（）02/04/20201602/04/202017舍选法)(:,,)(minmax1minxxrxxxxxf即变量产生均匀分布的由第一个随机取值范围的自变量根据概率密度函数max2max,0,fruf即范围内与均匀分布在变量产生第二个独立的随机。，，,),(新进行抽样从值拒绝该否则值则接受该如果xxxfu问题:如何找到函数的最大值?x的最大值)(xf)(xf02/04/202018舍选法举例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的最大值找)(xf30h/pl20;h/pl10;h/pl3;1zonePAWacc_rej.fcallPAW02/04/202019舍选法举例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;i100;i++){doubler=gRandom-Uniform(0,10);doublef=0.5*exp(-r/2.);if(fmaxf)fmax=f;}fmax*=1.2;coutfmax=fmaxendl;intntot=0;for(inti=0;i10000;i++){doubler=gRandom-Uniform(0,10);doublez=0.5*exp(-r/2.);if(zfmax)fmax=1.2*z;h11-Fill(r);doubleu=gRandom-Uniform(0,fmax);if(uz){h12-Fill(r);h2-Fill(r,u);ntot+=1;}}coutntot=ntotendl;}的最大值找)(xfDraw();h13-root[3]Draw();h12-root[2]Draw();h11-root[1]acc_rej.C.xroot[0]ROOT脚本02/04/202020舍选法举例(续))(rx随机变量)(rx随机变量)(rx随机变量频数)(xf频数边缘分布舍选法存在效率问题。二维均匀分布02/04/202021函数变换法与舍选法函数变换法优点：100%的抽样效率缺点：函数须解析可积舍选法优点：方法简单，可用于非常复杂的函数缺点：需要估计函数最大值，而且抽样效率低粒子物理与核物理中，对常用的概率密度函数有各种建议采用的方法(见)。除此之外，舍选法最为常用。初学者常犯的错误：对同一个过程做计算机模拟批处理，没有考虑在批处理结果中存在随机数的重复性。02/04/202022常用概率密度分布函数的抽样高斯(正态)分布{gROOT-Reset();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;ikUPDATE;i++){x=gRandom-Gaus(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环境下采用已有的分布，可以容