粒子物理与核物理实验中的数据分析lecture-8-最大似然法(II)

1粒子物理与核物理实验中的数据分析陈少敏清华大学第八讲：昀大似然法(II)2本讲要点„双参数情况下的昀大似然法举例„利用MINUIT数值求解昀小值„推广的昀大似然法„利用昀大似然法处理分区数据„检验昀大似然法的拟合优度„与贝叶斯参数估计之间的关系3例子:昀大似然法处理双参数考虑散射角分布例如：对于参数为α=0.5，β=0.5，探测器有效范围xmin=-0.95,xmax=0.95,产生n=2000个蒙特卡罗事例，研究用昀大似然法处理双参数的问题。θcos=x3/221),;(2ββαβα+++=xxxfminmaxxxx在实际情况下，探测器也许不可能做全空间探测，即。这个时候需要考虑效率影响，并对函数进行归一化处理，例如11(;,)()1fxxdxαβε−⋅=∫θminmax1()0xxxxε⎧=⎨⎩效率否则4双参数问题(续)用数值解找到logL的昀大值mˆ0.50ˆ0.51ˆˆcov,0.0024ˆ0.40.050.121rαβαβρ=±=±⎡⎤=⎣⎦==协方差m21ˆlog()ijijLVθθθθ−=∂=−∂∂GG“误差”：ˆˆˆˆ,αβσσ由联合概率密度得到样本的似然函数20001(;,)(;,)iiLxfxαβαβ==∏G(MINUITHESSE)注意：拟合中采用的是点拟合，与直方图的区间大小划分无关。注意：拟合中采用的是点拟合，与直方图的区间大小划分无关。5双参数拟合：蒙特卡罗研究做500次模拟实验，重复昀大似然拟合，每次都是模拟n=2000个事例：ˆˆˆˆ0.4990.4980.0510.111ˆˆcov[,]0.0024ˆ0.42ssrαβαβαβρ=======ˆˆ,αβ正相关边缘概率密度函数近似为高斯分布。ˆαˆβˆαˆβˆˆαβ表明与并不完全独立。6双参数拟合:logL等高线在前面蒙特卡罗样本中，其中的一次拟合结果在α，β平面表示为对于大的样本容量n，logL具有下列形式：max222ˆˆˆˆlog(,)logˆˆˆˆ122(1)LLααββαβααββααββρρσσσσ=⎡⎤⎛⎞⎛⎞⎛⎞⎛⎞−−−−⎢⎥⎜⎟⎜⎟−+−⎜⎟⎜⎟⎜⎟⎜⎟⎢⎥−⎝⎠⎝⎠⎝⎠⎝⎠⎣⎦7双参数拟合:logL等高线(续)等高线logL(α,β)=logLmax–½是根据下式定义得到的222ˆˆˆˆˆˆˆˆ1211ααββααββααββρρσσσσ⎡⎤⎛⎞⎛⎞⎛⎞⎛⎞−−−−⎢⎥⎜⎟⎜⎟+−=⎜⎟⎜⎟⎜⎟⎜⎟⎢⎥−⎝⎠⎝⎠⎝⎠⎝⎠⎣⎦其中，等高线的切线给出了标准偏差椭圆的倾角与相关系数有关βασσˆˆ,2ˆ2ˆˆˆ22tanβαβασσσρσφ−=注意：参数之间相关的效应体现在对估计量的误差或方差方面的影响(或偏大或偏小)。注意：参数之间相关的效应体现在对估计量的误差或方差方面的影响(或偏大或偏小)。8一个实用的求函数昀小值程序在核与粒子物理研究中，大家普遍使用的求函数昀小值程序是MINUIT软件包它可以对目标函数为似然函数，昀小二乘函数或用户自定义函数求极值。它提供了几种昀小化过程和相应的误差分析。原初版本用Fortran语言。写于25年前，用于当时西欧核子中心(CERN)的实验数据分析。现也被应用于粒子物理以外的领域。目前，在核与粒子物理研究中，有三种使用方式¾在MINUIT框架下单独使用(适于data-driven模式)；¾在物理分析工作站(PAW)环境下互动调用MINUIT软件包(基于Fortran)；¾在PAW升级软件包ROOT环境下互动调用MINUIT软件包(基于C++)。求logL昀大值等效于求–logL的昀小值。~chensm/lectures/minuit.pdf9MINUIT数值求解极小值(1)log()Lθ−G实际应用中，通常采用的形式来求昀小值。例如，对前面的例子programMINUIT_FITimplicitNONEintegernparparameter(npar=2)character*10chnam(npar)integerierr,ird,isav,istat,ivarbl,iwrintegernpari,nparxdoubleprecisionarglis(10),bnd1,bnd2,deriv(npar),dpar(npar)doubleprecisionfmin,fedm,errdef,covmat(npar,npar),log_lexternalFCNdoubleprecisionpar(npar)cInitializeMINUIT,setprintlevelto-1ird=5!unitnumberforinputtoMinuit(keyboard)iwr=6!unitnumberforoutputfromMinuit(screen)isav=7!unitnumberforuseoftheSAVEcommandcallMNINIT(ird,iwr,isav)arglis(1)=-1.d0callMNEXCM(FCN,'SETPRIN',arglis,1,ierr,0)cDefineparametersalphaandbeta,giveinitialvaluesandstepsizescallMNPARM(1,'alpha',0.5d0,0.1d0,0.d0,0.d0,ierr)callMNPARM(2,'beta',0.5d0,0.1d0,0.d0,0.d0,ierr)cGetinputdatabycallingFCNwithiflag=1arglis(1)=1.d0callMNEXCM(FCN,'CALL',arglis,1,ierr,0)cMinimizeusingSIMPLEXandMIGRAD,getcovariancecmatrixwithHESSEcallMNEXCM(FCN,'SIMPLEX',arglis,0,ierr,0)callMNEXCM(FCN,'MIGRAD',arglis,0,ierr,0)callMNEXCM(FCN,'HESSE',arglis,0,ierr,0)cGetresultoffit(forleastsquares,fminischi2)callMNSTAT(fmin,fedm,errdef,npari,nparx,istat)callMNPOUT(1,chnam(1),par(1),dpar(1),bnd1,bnd2,ivarbl)callMNPOUT(2,chnam(2),par(2),dpar(2),bnd1,bnd2,ivarbl)callMNEMAT(covmat,npar)log_l=-0.5*fminwrite(6,*)'Fitresults:'write(6,*)write(6,*)'alpha=',par(1),'+/-',dpar(1)write(6,*)'beta=',par(2),'+/-',dpar(2)write(6,*)'cov[alpha,beta]=',covmat(1,2)write(6,*)'rho[alpha,beta]=',covmat(1,2)/(dpar(1)*dpar(2))write(6,*)'log_l=',log_lstopend编译:f77minuit_fit.f-ominuit_fit`cernlib`return运行:minuit_fitreturn10MINUIT数值求解极小值(2)subroutineFCN(npar,grad,chi2,par,iflag,futil)cInput:integernparnumberofparameterstofitcdoubleprecisionpar(npar)parametervectorcintegeriflagselectwhattodocdoubleprecisionfutiloptionalexternalfunctioncOutput:doubleprecisiongrad(npar)gradientvector(notfilled)cdoubleprecisionchi2functiontobeminimizedimplicitNONEintegernpar,iflagdoubleprecisionfutil,chi2,par(*),grad(*)integern_maxparameter(n_max=10000)integeri,ndoubleprecisionalpha,beta,f,log_l,x(n_max),angle_cut,f_norcBeginn=2000angle_cut=0.95if(iflag.eq.1)then!getn,arrayxcallGET_INPUT_DATA(x,n,n_max,angle_cut)endifcCalculatelog-likelihoodalpha=par(1)beta=par(2)log_l=0.cNormalizationfactorfor[-angle_cut,angle_cut]f_nor=2*angle_cut+2*beta/3.*angle_cut**3doi=1,nf=(1.+alpha*x(i)+beta*x(i)**2)/f_norlog_l=log_l+DLOG(f)enddochi2=-2.*log_l!2getserrorsrightreturnend用户应提供似然函数FCN注意：概率密度函数需要在有效测量范围进行归一化处理。如果计算有困难，可以将归一化因子作为参数来拟合。这种折衷方法虽简便，但会给其它待定参数的确定带来误差。11MINUIT数值求解极小值(3)subroutineget_input_data(x,n,n_max,angle_cut)integern,n_max,ntotdoubleprecisionx(n_max),angle_cutrealrvec(1),alpha,beta,r,fmax,z,ualpha=0.5beta=0.5fmax=-999.doi=1,100callranmar(rvec,1)r=rvec(1)*2.0-1.0!from(0,1)to(-1,1)f=(1.+alpha*r+beta*r**2)/(2.+2.*beta/3.)if(fmax.lt.f)fmax=fenddofmax=1.05*fmaxntot=010callranmar(rvec,1)r=rvec(1)*2.0-1.0!from(0,1)to(-1,1)z=(1.+alpha*r+beta*r**2)/(2.+2.*beta/3.)if(z.gt.fmax)thenfmax=zwrite(6,*)'zgreaterthanfmax'endifcallranmar(rvec,1)u=rvec(1)*fmaxif(u.le.z)thenif(abs(r).lt.angle_cut)thenntot=ntot+1x(ntot)=rif(ntot.ge.n)goto20endifendifgoto1020continuereturnend对于采用蒙特卡罗研究参数估计值问题，可用下面的舍选法例子对任意概率密度函数进行抽样，得到模拟的数据样本。12MINUIT数值求解极小值(4)做500次模拟实验，重复昀大似然拟合，每次都是模拟n=2000个事例：programMINUIT_FIT…doubleprecisionpar(npar)realh(80000)common/pawc/hintegericallhlimit(80000)callhbook2(200,'alphavsbeta',100,0.,1.,100,0.,1.,0.)doi=1,500cInitializeMINUIT,setprintlevelto–1…callhfill(200,real(par(1)),real(par(2)),1.0)enddocallhrput(0,'minuit_fit.hbook'