清华大学：《粒子物理与核物理实验中的数据分析》课程教学资源（课件讲稿）第八讲 最大似然法（二）

双参数拟合：logL等高线 在前面蒙特卡罗样本中，其中的 0.7 次拟合结果在，B平面表示为 0.6 0.5 true value ML fit result 0.4 对于大的样本容量n,logL具有下列 0.3 形式： 0.3 0.4 0.5 0.6 0.7 C log L(a,B)=log Lmax =川84 6

6 双参数拟合 :log L 等高线 在前面蒙特卡罗样本中，其中的 一次拟合结果在 α， β 平面表示为 对于大的样本容量 n，log L 具有下列 形式： ma x 2 2 2 ˆ ˆ ˆ ˆ lo g ( , ) lo g ˆ ˆ 1 ˆ ˆ 2 2( 1 ) L L α α β β α β α α β β α α β β ρ ρ σ σ σ σ = ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ − − ⎛ ⎞ − − ⎢ ⎥ − + ⎜ ⎟ − ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ − ⎢ ⎜ ⎟ ⎜ ⎟⎥ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦

双参数拟合logL等高线（续） 等高线logL(a,β)-log L max-么是根据下式定义得到的 j 其中，等高线的切线给出了标准偏差 注意：参数之间相关的效应体 椭圆的倾角与相关系数有关 现在对估计量的误差或方差方 面的影响（或偏大或偏小）。 2poa tan2φ= 7

7 双参数拟合 :log L 等高线 ( 续 ) 等高线 log L ( α,β)=log L ma x – ½ 是根据下式定义得到的 2 2 2 ˆ ˆ ˆ ˆ ˆ ˆ 1 ˆ ˆ 2 1 1 α α β β α α β β α α β β ρ ρ σ σ σ σ ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ − − ⎛ ⎞ − − ⎢ ⎥ + ⎜ ⎟ − = ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ − ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ 其中，等高线的切线给出了标准偏差 椭圆的倾角与相关系数有关 α β σ ˆ σ ˆ , 2 ˆ 2 ˆ 2 ˆ ˆ tan 2 α β α β σ σ ρσ σ φ − = 注意：参数之间相关的效应体 现在对估计量的误差或方差方 面的影响 (或偏大或偏小 )。 注意：参数之间相关的效应体 现在对估计量的误差或方差方 面的影响 (或偏大或偏小 )

一个实用的求函数最小值程序 在核与粒子物理研究中，大家普遍使用的求函数最小值程序是 MINUIT软件包 http://hep.tsinghua.edu.cn/~chensm/lectures/minuit.pdf 求logL最大值等效于求-logL的最小值。 它可以对目标函数为似然函数，最小二乘函数或用户自定义函数求极值。 它提供了几种最小化过程和相应的误差分析。原初版本用Fortran语言。 写于25年前，用于当时西欧核子中心(CERN)的实验数据分析。现也被 应用于粒子物理以外的领域。目前，在核与粒子物理研究中，有三种使 用方式 >在MINUIT框架下单独使用（适于data-driven模式）： >在物理分析工作站(PAW)环境下互动调用MINUIT软件包（基于Fortran); >在PAW升级软件包ROOT环境下互动调用MINUIT软件包（基于C++)。 8

8 一个实用的求函数最小值程序 在核与粒子物理研究中，大家普遍使用的求函数最小值程序是 MINUIT 软件包 它可以对目标函数为似然函数，最小二乘函数或用户自定义函数求极值。 它提供了几种最小化过程和相应的误差分析。原初版本用Fortran语言。 写于25年前，用于当时西欧核子中心 (CERN )的实验数据分析。现也被 应用于粒子物理以外的领域。目前，在核与粒子物理研究中，有三种使 用方式 ¾ 在MINUIT框架下单独使用 (适于data-driven模式 )； ¾在物理分析工作站(PAW)环境下互动调用MINUIT软件包 (基于Fortran)； ¾ 在PAW升级软件包ROOT环境下互动调用MINUIT软件包 (基于C++ )。 求 log L 最大值等效于求 –log L 的最小值。 http://hep.tsinghua.edu.cn/~chensm/lectures/minuit.pdf

MINUIT数值求解极小值(1) 实际应用中，通常采用-ogL()的形式来求最小值。例如，对前面的例子 program MINUIT FIT c Minimize using SIMPLEX and MIGRAD,get covariance implicit NONE c matrix with HESSE integer npar call MNEXCM(FCN,'SIMPLEX',arglis,0,ierr,0) parameter(npar=2) call MNEXCM(FCN,'MIGRAD,arglis,0,ierr,0) character*10 chnam(npar) call MNEXCM(FCN,'HESSE,arglis,0,ierr,0) integer ierr,ird,isav,istat,ivarbl,iwr c Get result of fit(for least squares,fmin is chi2) integer npari,nparx call MNSTAT(fmin,fedm,errdef,npari,nparx,istat) double precision arglis(10),bndl,bnd2,deriv(npar),dpar(npar) call MNPOUT(1,chnam(1),par(1),dpar(1),bnd1,bnd2,ivarbl) double precision fmin,fedm,errdef,covmat(npar,npar),log I call MNPOUT(2,chnam(2),par(2),dpar(2),bnd1,bnd2,ivarbl) external FCN call MNEMAT(covmat,npar) double precision par(npar) log 1=-0.5*fmin c Initialize MINUIT,set print level to-1 write(6,*)'Fit results:' ird=5!unit number for input to Minuit(keyboard) write(6,*) iwr=6!unit number for output from Minuit(screen) write(6,*)'alpha =,par(1),'+/-',dpar(1) isav=7 unit number for use of the SAVE command write(6,*)'beta =,par(2),'+-',dpar(2) call MNINIT(ird,iwr,isav) write(6,*)'cov[alpha,beta]=',covmat(1,2) arglis(1)=-1.d0 write(6,*)'rho[alpha,beta]=',covmat(1,2)/(dpar(1)*dpar(2)) call MNEXCM(FCN,'SET PRIN,arglis,1,ierr,0) write(6,*)'log 1 =,log 1 c Define parameters alpha and beta,give initial values and step sizes stop call MNPARM(1,alpha,0.5d0,0.1d0,0.d0,0.d0,ierr) end call MNPARM(2,beta,0.5d0,0.1d0,0.d0,0.d0,ierr) c Get input data by calling FCN with iflag=l arglis(1)=1.do 编译：f77 minuit fit,f-o minuit fit`cernlib<return> call MNEXCM(FCN,'CALL,arglis,1,ierr,0) 运行：minuit fit<return>

9 MINUIT数值求解极小值(1) − log L(θ ) G 实际应用中，通常采用 的形式来求最小值。例如，对前面的例子 program MINUIT_FIT implicit NONE integer npar parameter (npar=2) character*10 chnam(npar) integer ierr, ird, isav, istat, ivarbl, iwr integer npari, nparx double precision arglis(10), bnd1, bnd2, deriv(npar), dpar(npar) double precision fmin, fedm, errdef, covmat(npar,npar), log_l external FCN double precision par(npar) c Initialize MINUIT, set print level to -1 ird = 5 ! unit number for input to Minuit (keyboard) iwr = 6 ! unit number for output from Minuit (screen) isav= 7 ! unit number for use of the SAVE command call MNINIT(ird,iwr,isav) arglis(1)=-1.d0 call MNEXCM(FCN,'SET PRIN',arglis,1,ierr,0) c Define parameters alpha and beta, give initial values and step sizes call MNPARM(1,'alpha',0.5d0,0.1d0,0.d0,0.d0,ierr) call MNPARM(2,'beta ',0.5d0,0.1d0,0.d0,0.d0,ierr) c Get input data by calling FCN with iflag=1 arglis(1)=1.d0 call MNEXCM(FCN,'CALL',arglis,1,ierr,0) c Minimize using SIMPLEX and MIGRAD, get covariance c matrix with HESSE call MNEXCM(FCN,'SIMPLEX',arglis,0,ierr,0) call MNEXCM(FCN,'MIGRAD',arglis,0,ierr,0) call MNEXCM(FCN,'HESSE',arglis,0,ierr,0) c Get result of fit (for least squares, fmin is chi2) call MNSTAT(fmin,fedm,errdef,npari,nparx,istat) call MNPOUT(1,chnam(1),par(1),dpar(1),bnd1,bnd2,ivarbl) call MNPOUT(2,chnam(2),par(2),dpar(2),bnd1,bnd2,ivarbl) call MNEMAT(covmat,npar) log_l=-0.5*fmin write(6,*)'Fit results:' 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_l stop end 编译:f77 minuit_fit.f -o minuit_fit `cernlib` <return> 运行:minuit_fit <return>

MINUIT数值求解极小值(2) 用户应提供似然函数FCN 注意：概率密度函数需要在有 subroutine FCN(npar,grad,chi2,par,iflag,futil) 效测量范围进行归一化处理。 c Input:integer npar number of parameters to fit 如果计算有困难，可以将归一 c double precision par(npar)parameter vector C integer iflag select what to do 化因子作为参数来拟合。这种 c double precision futil optional external function 折衷方法虽简便，但会给其它 c Output:double precision grad(npar)gradient vector(not filled) double precision chi2 function to be minimized 待定参数的确定带来误差。 implicit NONE c Calculate log-likelihood integer npar,iflag alpha=par(1) double precision futil,chi2,par(*),grad(*) beta =par(2) integer n max parameter(n_max=10000) 1og1=0. c Normalization factor for [-angle_cut,angle_cut] integer i,n f nor =2*angle cut+2*beta/3.*angle cut**3 double precision alpha,beta,f,log 1,x(n_max),angle_cut,f_nor do i=1,n c Begin n=2000 f-(1.+alpha*x(i)+beta*x(i)**2)/f_nor angle_cut=0.95 log_1=log_1+DLOG(f) end do if(iflag.eq.1)then get n,array x call GET INPUT_DATA(x,n,n_max,angle_cut) chi2=-2.*log_1 2 gets errors right return end if end 10

10 MINUIT数值求解极小值(2) subroutine FCN(npar,grad,chi2,par,i flag,futil) c Input: integ er npar number of parameters to fit c double precision par(npar) parameter vector c integer iflag s elect what to do c double precision futil optional external function c Output: double precision grad(npar) gradient vector (not filled) c double precision chi2 function to be minimized implicit NONE integ er npar,iflag double precision futil,chi2,par(*),grad(*) integ e r n_max parameter (n_max=10000) integ e r i,n double precision alpha,beta,f,log_l,x(n_max),angle_cut,f_nor c Begin n=2000 angle_cut=0.95 if(iflag.eq.1)then ! get n, array x call GET_INPUT_DATA(x,n,n_max,angle_cut) end if c Calculate log-likelihood alpha = par(1) beta = par(2) log_l = 0. c Normalization factor for [-angle_cut,angle_cut] f_nor = 2*angle_cut+2*beta/3.*angle_cut**3 do i=1,n f=(1. +alpha*x(i)+beta*x(i)**2)/ f_nor log_l=log_l+DLOG(f) end do chi2=-2.*log_l ! 2 gets errors right return end 用户应提供似然函数 FCN 注意：概率密度函数需要在有 效测量范围进行归一化处理。 如果计算有困难，可以将归一 化因子作为参数来拟合。这种 折衷方法虽简便，但会给其它 待定参数的确定带来误差