c BAYES - program to check equivalency of Bayes and frequentist c upper limit calculations c c John Conway Nov. 2004 implicit none integer nev ! number of observed events real*8 sac ! relative uncertainty in acceptance real*8 xbg ! expected background real*8 sbg ! uncertainty in background real*8 cl ! desired CL real*8 smax ! upper limit of integration real*8 prec ! integration step size real*8 plim ! upper limit on Poisson process real*8 perr ! error on upper limit on Poisson process logical done character*1 yes print '('' ------------------------------------------------- '')' print '('' Bayes version 2 (Nov. 2004) '')' print '('' ------------------------------------------------- '')' print '('' This program, bayes.f, performs a calculation of '')' print '('' the upper limit on a Poisson process with '')' print '('' background, incorporating the effects of '')' print '('' uncertainty on the background and signal '')' print '('' acceptance. The program assumes a uniform prior '')' print '('' pdf in the signal rate, and Gaussian priors for '')' print '('' the acceptance and expected background. The '')' print '('' Gaussian priors are truncated at zero. Formally, '')' print '('' therefore, the integral of the posterior density '')' print '('' with respect to the signal rate diverges '')' print '('' logarithmically. This is a small effect '')' print '('' (requiring an integral out to very large signal '')' print '('' rates to be observed) if the acceptance '')' print '('' uncertainty is small, less than about 15-20%. If '')' print '('' your acceptance uncertainty is larger than this '')' print '('' your resulting limit may be sensitive to the upper'')' print '('' limit of integration. In this case, you have a '')' print '('' problem, and you should consider an alternative '')' print '('' approach. '')' print '('' '')' print '('' If you have any questions about these issues '')' print '('' please contact John Conway at conway@fnal.gov. '')' print '('' ------------------------------------------------- '')' print '('' Be careful to enter a real number (with a decimal '')' print '('' point) for the values below requesting it. '')' print '('' ------------------------------------------------- '')' print '(\$,'' enter number of observed events (integer): '')' read(*,'(i4)') nev print '(\$,'' enter relative error on acceptance (real): '')' read(*,'(f10.4)') sac if(sac.gt.0.15) then print '('' Your error is quite large - please contact '')' print '('' John Conway at conway@fnal.gov ! '')' endif print '(\$,'' enter expected number of background events '', . ''(real): '')' read(*,'(f10.4)') xbg print '(\$,'' enter error on number of background events '', . ''(real): '')' read(*,'(f10.4)') sbg print '(\$,'' enter desired confidence level (real): '')' read(*,'(f10.4)') cl smax = 10*float(nev) if(smax.lt.10.0) smax = 10.0 done = .false. do while(.not.done) print '(\$,'' enter integration upper limit (real): '')' read(*,'(f10.4)') smax if(smax.lt.10.0) then print '(/,'' Value of upper limit too small! '')' elseif(smax.lt.float(nev)+5.0*sqrt(float(nev))) then print '(/,'' Value of upper limit too small! '')' else done = .true. endif enddo print '(\$,'' enter integration step size '', . ''(recommend 0.02 or less):'')' read(*,'(f10.4)') prec print '(/)' print '('' ----------------------------------------- '')' print '('' observed events '',i10)',nev print '('' relative error on acceptance '',f10.3)',sac print '('' expected background '',f10.3)',xbg print '('' absolute error on background '',f10.3)',sbg print '('' desired confidence level '',f10.2)',cl print '('' integration upper limit '',f10.2)',smax print '('' integration step size '',f10.4)',prec print '('' ----------------------------------------- '')' print '(/)' print '(\$,'' Are the above correct? '')' read(*,'(a)') yes if(yes(1:1).ne.'y') stop call baylim(nev,sac,xbg,sbg,cl,prec,smax,plim,perr) print '(//,'' limit: less than'', . f10.3,'' signal events '',//)',plim end subroutine baylim(nev,sac,xbg,sbg,cl,prec,smax,plim,perr) c calculate Bayesian upper limit on Poisson process implicit none c arguments integer nev ! number of events real*8 sac ! error on acceptance real*8 xbg ! expected background real*8 sbg ! error on background real*8 cl ! desired CL real*8 smax ! upper limit of integration real*8 prec ! integration step size real*8 plim ! Poisson upper limit (returned) real*8 perr ! error on upper limit (returned) real*8 xev,xevmax,dxev logical done real*8 blike,xlike integer i,nlist,icl real*8 blist(10000),xlist(10000),bint,bsum,bcl c scan likelihood L(nev|xev,...) xevmax = smax dxev = prec nlist = 0 bsum = 0. done = .false. xev = dxev/2. do while(.not.done) xlike = blike(nev,sac,xbg,sbg,xev) nlist = nlist + 1 xlist(nlist) = xev blist(nlist) = xlike bsum = bsum + xlike print '(\$,i5,2g12.5,a)',nev,xev,xlike,char(13) xev = xev + dxev if(blist(nlist)/blist(1).lt.1.0e-6) done = .true. if(xev.gt.xevmax) done = .true. enddo c find place just below threshold bint = 0. icl = 0 do i=1,nlist if(bint .le.cl*bsum.and. . bint+blist(i).gt.cl*bsum) then icl = i bcl = bint endif bint = bint + blist(i) enddo c interpolate linearly plim = xlist(icl) + (xlist(icl+1)-xlist(icl))* . (cl*bsum-bcl)/blist(icl+1) perr = 0. return end real*8 function blike(nev,sac,xbg,sbg,xev) c return likelihood to observe nev events given expected implicit none integer nev ! number of events observed real*8 xev ! number of events expected real*8 sac ! error on acceptance real*8 xbg ! expected background real*8 sbg ! error on background integer i,nmax parameter(nmax=100000) real*4 a,b real*8 yyev,ssev,yyex,xxev,ssac,xxbg,ssbg,yybg,xint,xxx real*8 stirling logical done c perform double gaussian integral by Monte Carlo xxev = nev xxbg = xbg ssbg = sbg ssac = sac xxev = xev ssev = xxev*ssac xint = 0. done = .false. do i=1,nmax c pick expected background and signal from Gaussian yybg = -1. yyev = -1. do while(yybg.lt.0..or. . yyev.lt.0.) call rannor(a,b) yybg = xxbg + a*ssbg yyev = xxev*(1.0+b*ssac) enddo c total expected yyex = yybg + yyev c value of integrand xxx = yyex**nev * dexp(-yyex) / stirling(nev) ! Poisson term xint = xint + xxx enddo blike = xint/float(nmax) return end real*8 function stirling(n) c Stirling's approximation for n! with first correction term real*8 x,fact(0:9) parameter (root2pi=2.506628) data fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./ if(n.ge.0.and.n.le.9) then stirling=fact(n) return endif if(n.gt.0) then x=n stirling=root2pi*sqrt(x)*x**n*exp(-x)*(1.+1./12./x) else stirling=1. endif return end