PDF《粒子物理与核物理实验中的》

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` 运行:minuit_fit

10 MINUIT数值求解极小值(2) 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) c double precision chi2 function to be minimized implicit NONE integer npar,iflag double precision futil,chi2,par(*),grad(*) integer n_max parameter (n_max=10000) integer 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 注意:概率密度函数需要在有 效测量范围进行归一化处理. 如果计算有困难,可以将归一 化因子作为参数来拟合.这种 折衷方法虽简便,但会给其它 待定参数的确定带来误差.

11 MINUIT数值求解极小值(3) subroutine get_input_data(x,n,n_max,angle_cut) integer n,n_max,ntot double precision x(n_max),angle_cut real rvec(1),alpha,beta,r,fmax,z,u alpha=0.5 beta =0.5 fmax=-999. do i=1,100 call ranmar(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=f end do fmax=1.05*fmax ntot=0