Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 8085) +++ trunk/MagicSoft/Mars/Changelog (revision 8087) @@ -18,4 +18,14 @@ -*-*- END OF LINE -*-*- + 2006/10/17 Markus Meyer + + * mtools + - added a new class called MRolke, which is a modification of + TRolke from root_v5.12.00b. There is now a new function, called + LogFactorial() which enables to calculate confidence intervals + even for a large number of events (larger than 170). + + + 2006/10/17 Stefan Ruegamer Index: trunk/MagicSoft/Mars/mtools/MRolke.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MRolke.cc (revision 8087) +++ trunk/MagicSoft/Mars/mtools/MRolke.cc (revision 8087) @@ -0,0 +1,799 @@ +// @(#)root/physics:$Name: not supported by cvs2svn $:$Id: MRolke.cc,v 1.1 2006-10-17 12:51:31 meyer Exp $ +// Author: Jan Conrad 9/2/2004 + +/************************************************************************* + * Copyright (C) 1995-2004, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////////// +// +// MRolke +// +// This class computes confidence intervals for the rate of a Poisson +// in the presence of background and efficiency with a fully frequentist +// treatment of the uncertainties in the efficiency and background estimate +// using the profile likelihood method. +// +// The signal is always assumed to be Poisson. +// +// The method is very similar to the one used in MINUIT (MINOS). +// +// Two options are offered to deal with cases where the maximum likelihood +// estimate (MLE) is not in the physical region. Version "bounded likelihood" +// is the one used by MINOS if bounds for the physical region are chosen. Versi// on "unbounded likelihood (the default) allows the MLE to be in the +// unphysical region. It has however better coverage. +// For more details consult the reference (see below). +// +// +// It allows the following Models: +// +// 1: Background - Poisson, Efficiency - Binomial (cl,x,y,z,tau,m) +// 2: Background - Poisson, Efficiency - Gaussian (cl,xd,y,em,tau,sde) +// 3: Background - Gaussian, Efficiency - Gaussian (cl,x,bm,em,sd) +// 4: Background - Poisson, Efficiency - known (cl,x,y,tau,e) +// 5: Background - Gaussian, Efficiency - known (cl,x,y,z,sdb,e) +// 6: Background - known, Efficiency - Binomial (cl,x,z,m,b) +// 7: Background - known, Efficiency - Gaussian (cl,x,em,sde,b) +// +// Parameter definition: +// +// cl = Confidence level +// +// x = number of observed events +// +// y = number of background events +// +// z = number of simulated signal events +// +// em = measurement of the efficiency. +// +// bm = background estimate +// +// tau = ratio between signal and background region (in case background is +// observed) ratio between observed and simulated livetime in case +// background is determined from MC. +// +// sd(x) = sigma of the Gaussian +// +// e = true efficiency (in case known) +// +// b = expected background (in case known) +// +// m = number of MC runs +// +// mid = ID number of the model ... +// +// For a description of the method and its properties: +// +// W.Rolke, A. Lopez, J. Conrad and Fred James +// "Limits and Confidence Intervals in presence of nuisance parameters" +// http://lanl.arxiv.org/abs/physics/0403059 +// Nucl.Instrum.Meth.A551:493-503,2005 +// +// Should I use MRolke, TFeldmanCousins, TLimit? +// ============================================ +// 1. I guess MRolke makes TFeldmanCousins obsolete? +// +// Certainly not. TFeldmanCousins is the fully frequentist construction and +// should be used in case of no (or negligible uncertainties). It is however +// not capable of treating uncertainties in nuisance parameters. +// MRolke is desined for this case and it is shown in the reference above +// that it has good coverage properties for most cases, ie it might be +// used where FeldmannCousins can't. +// +// 2. What are the advantages of MRolke over TLimit? +// +// MRolke is fully frequentist. TLimit treats nuisance parameters Bayesian. +// For a coverage study of a Bayesian method refer to +// physics/0408039 (Tegenfeldt & J.C). However, this note studies +// the coverage of Feldman&Cousins with Bayesian treatment of nuisance +// parameters. To make a long story short: using the Bayesian method you +// might introduce a small amount of over-coverage (though I haven't shown it +// for TLimit). On the other hand, coverage of course is a not so interesting +// when you consider yourself a Bayesian. +// +// Author: Jan Conrad (CERN) +// +// see example in tutorial Rolke.C +// +// Copyright CERN 2004 Jan.Conrad@cern.ch +// +/////////////////////////////////////////////////////////////////////////// + + +#include "MRolke.h" +#include "TMath.h" +#include "Riostream.h" + +ClassImp(MRolke) + +//__________________________________________________________________________ +MRolke::MRolke(Double_t CL, Option_t * /*option*/) +{ + //constructor + fUpperLimit = 0.0; + fLowerLimit = 0.0; + fCL = CL; + fSwitch = 0; // 0: unbounded likelihood + // 1: bounded likelihood +} + +//___________________________________________________________________________ +MRolke::~MRolke() +{ +} + + +//___________________________________________________________________________ +Double_t MRolke::CalculateInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em,Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m) +{ + //calculate interval + Int_t done = 0; + Double_t limit[2]; + + limit[1] = Interval(x,y,z,bm,em,e,mid, sde,sdb,tau,b,m); + + if (limit[1] > 0) { + done = 1; + } + + if (fSwitch == 0) { + + Int_t trial_x = x; + + while (done == 0) { + trial_x++; + limit[1] = Interval(trial_x,y,z,bm,em,e,mid, sde,sdb,tau,b,m); + if (limit[1] > 0) done = 1; + } + } + + return limit[1]; +} + + + + +//_____________________________________________________________________ +Double_t MRolke::Interval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em,Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m) +{ + // Calculates the Confidence Interval + + //Double_t dchi2 = Chi2Percentile(1,1-fCL); + Double_t dchi2 = TMath::ChisquareQuantile(fCL, 1); + + Double_t tempxy[2],limits[2] = {0,0}; + Double_t slope,fmid,low,flow,high,fhigh,test,ftest,mu0,maximum,target,l,f0; + Double_t med = 0; + Double_t maxiter=1000, acc = 0.00001; + Int_t i; + Int_t bp = 0; + + if ((mid != 3) && (mid != 5)) bm = (Double_t)y; + + if ((mid == 3) || (mid == 5)) { + if (bm == 0) bm = 0.00001; + } + + if ((mid <= 2) || (mid == 4)) bp = 1; + + + if (bp == 1 && x == 0 && bm > 0 ){ + + for(Int_t i = 0; i < 2; i++) { + x++; + tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m); + } + slope = tempxy[1] - tempxy[0]; + limits[1] = tempxy[0] - slope; + limits[0] = 0.0; + if (limits[1] < 0) limits[1] = 0.0; + goto done; + } + + if (bp != 1 && x == 0){ + + for(Int_t i = 0; i < 2; i++) { + x++; + tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m); + } + slope = tempxy[1] - tempxy[0]; + limits[1] = tempxy[0] - slope; + limits[0] = 0.0; + if (limits[1] < 0) limits[1] = 0.0; + goto done; + } + + if (bp != 1 && bm == 0){ + for(Int_t i = 0; i < 2; i++) { + bm++; + limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m); + tempxy[i] = limits[1]; + } + slope = tempxy[1] - tempxy[0]; + limits[1] = tempxy[0] - slope; + if (limits[1] < 0) limits[1] = 0; + goto done; + } + + + if (x == 0 && bm == 0){ + x++; + bm++; + + limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m); + tempxy[0] = limits[1]; + x = 1; + bm = 2; + limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m); + tempxy[1] = limits[1]; + x = 2; + bm = 1; + limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m); + limits[1] = 3*tempxy[0] -tempxy[1] - limits[1]; + if (limits[1] < 0) limits[1] = 0; + goto done; + } + + mu0 = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,1); + + maximum = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,2); + + test = 0; + + f0 = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3); + + if ( fSwitch == 1 ) { // do this only for the unbounded likelihood case + if ( mu0 < 0 ) maximum = f0; + } + + target = maximum - dchi2; + + if (f0 > target) { + limits[0] = 0; + } else { + if (mu0 < 0){ + limits[0] = 0; + limits[1] = 0; + } + + low = 0; + flow = f0; + high = mu0; + fhigh = maximum; + + for(Int_t i = 0; i < maxiter; i++) { + l = (target-fhigh)/(flow-fhigh); + if (l < 0.2) l = 0.2; + if (l > 0.8) l = 0.8; + + med = l*low + (1-l)*high; + if(med < 0.01){ + limits[1]=0.0; + goto done; + } + + fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3); + + if (fmid > target) { + high = med; + fhigh = fmid; + } else { + low = med; + flow = fmid; + } + if ((high-low) < acc*high) break; + } + limits[0] = med; + } + + + if(mu0 > 0) { + low = mu0; + flow = maximum; + } else { + low = 0; + flow = f0; + } + + test = low +1 ; + + ftest = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3); + + if (ftest < target) { + high = test; + fhigh = ftest; + } else { + slope = (ftest - flow)/(test - low); + high = test + (target -ftest)/slope; + fhigh = Likelihood(high,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3); + } + + for(i = 0; i < maxiter; i++) { + l = (target-fhigh)/(flow-fhigh); + if (l < 0.2) l = 0.2; + if (l > 0.8) l = 0.8; + med = l * low + (1.-l)*high; + fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3); + + if (fmid < target) { + high = med; + fhigh = fmid; + } else { + low = med; + flow = fmid; + } + if (high-low < acc*high) break; + } + limits[1] = med; + +done: + + if ( (mid == 4) || (mid==5) ) { + limits[0] /= e; + limits[1] /= e; + } + + + fUpperLimit = limits[1]; + fLowerLimit = TMath::Max(limits[0],0.0); + + + return limits[1]; +} + + +//___________________________________________________________________________ +Double_t MRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm,Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m, Int_t what) +{ + // Chooses between the different profile likelihood functions to use for the + // different models. + // Returns evaluation of the profile likelihood functions. + + switch (mid) { + case 1: return EvalLikeMod1(mu,x,y,z,e,tau,b,m,what); + case 2: return EvalLikeMod2(mu,x,y,em,e,sde,tau,b,what); + case 3: return EvalLikeMod3(mu,x,bm,em,e,sde,sdb,b,what); + case 4: return EvalLikeMod4(mu,x,y,tau,b,what); + case 5: return EvalLikeMod5(mu,x,bm,sdb,b,what); + case 6: return EvalLikeMod6(mu,x,z,e,b,m,what); + case 7: return EvalLikeMod7(mu,x,em,e,sde,b,what); + } + + return 0; +} + +//_________________________________________________________________________ +Double_t MRolke::EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t e, Double_t tau, Double_t b, Int_t m, Int_t what) +{ + // Calculates the Profile Likelihood for MODEL 1: + // Poisson background/ Binomial Efficiency + // what = 1: Maximum likelihood estimate is returned + // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned. + // what = 3: Profile Likelihood of Test hypothesis is returned + // otherwise parameters as described in the beginning of the class) + + Double_t f = 0; + Double_t zm = Double_t(z)/m; + + if (what == 1) { + f = (x-y/tau)/zm; + } + + if (what == 2) { + mu = (x-y/tau)/zm; + b = y/tau; + Double_t e = zm; + f = LikeMod1(mu,b,e,x,y,z,tau,m); + } + + if (what == 3) { + if (mu == 0){ + b = (x+y)/(1.0+tau); + e = zm; + f = LikeMod1(mu,b,e,x,y,z,tau,m); + } else { + MRolke g; + g.ProfLikeMod1(mu,b,e,x,y,z,tau,m); + f = LikeMod1(mu,b,e,x,y,z,tau,m); + } + } + + return f; +} + +//________________________________________________________________________ +Double_t MRolke::LikeMod1(Double_t mu,Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m) +{ + // Profile Likelihood function for MODEL 1: + // Poisson background/ Binomial Efficiency + + return 2*(x*TMath::Log(e*mu+b)-(e*mu +b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) + TMath::Log(TMath::Factorial(m)) - TMath::Log(TMath::Factorial(m-z)) - TMath::Log(TMath::Factorial(z))+ z * TMath::Log(e) + (m-z)*TMath::Log(1-e)); +} + +//________________________________________________________________________ +void MRolke::ProfLikeMod1(Double_t mu,Double_t &b,Double_t &e,Int_t x,Int_t y, Int_t z,Double_t tau,Int_t m) +{ + // Void needed to calculate estimates of efficiency and background for model 1 + + Double_t med = 0.0,fmid; + Int_t maxiter =1000; + Double_t acc = 0.00001; + Double_t emin = ((m+mu*tau)-TMath::Sqrt((m+mu*tau)*(m+mu*tau)-4 * mu* tau * z))/2/mu/tau; + + Double_t low = TMath::Max(1e-10,emin+1e-10); + Double_t high = 1 - 1e-10; + + for(Int_t i = 0; i < maxiter; i++) { + med = (low+high)/2.; + + fmid = LikeGradMod1(med,mu,x,y,z,tau,m); + + if(high < 0.5) acc = 0.00001*high; + else acc = 0.00001*(1-high); + + if ((high - low) < acc*high) break; + + if(fmid > 0) low = med; + else high = med; + } + + e = med; + Double_t eta = Double_t(z)/e -Double_t(m-z)/(1-e); + + b = Double_t(y)/(tau -eta/mu); +} + +//___________________________________________________________________________ +Double_t MRolke::LikeGradMod1(Double_t e, Double_t mu, Int_t x,Int_t y,Int_t z,Double_t tau,Int_t m) +{ + //gradient model + Double_t eta, etaprime, bprime,f; + eta = static_cast(z)/e - static_cast(m-z)/(1.0 - e); + etaprime = (-1) * (static_cast(m-z)/((1.0 - e)*(1.0 - e)) + static_cast(z)/(e*e)); + Double_t b = y/(tau - eta/mu); + bprime = (b*b * etaprime)/mu/y; + f = (mu + bprime) * (x/(e * mu + b) - 1)+(y/b - tau) * bprime + eta; + return f; +} + +//___________________________________________________________________________ +Double_t MRolke::EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t e,Double_t sde, Double_t tau, Double_t b, Int_t what) +{ + // Calculates the Profile Likelihood for MODEL 2: + // Poisson background/ Gauss Efficiency + // what = 1: Maximum likelihood estimate is returned + // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned. + // what = 3: Profile Likelihood of Test hypothesis is returned + // otherwise parameters as described in the beginning of the class) + + Double_t v = sde*sde; + Double_t coef[4],roots[3]; + Double_t f = 0; + + if (what == 1) { + f = (x-y/tau)/em; + } + + if (what == 2) { + mu = (x-y/tau)/em; + b = y/tau; + e = em; + f = LikeMod2(mu,b,e,x,y,em,tau,v); + } + + if (what == 3) { + if (mu == 0 ) { + b = (x+y)/(1+tau); + f = LikeMod2(mu,b,e,x,y,em,tau,v); + } else { + coef[3] = mu; + coef[2] = mu*mu*v-2*em*mu-mu*mu*v*tau; + coef[1] = ( - x)*mu*v - mu*mu*mu*v*v*tau - mu*mu*v*em + em*mu*mu*v*tau + em*em*mu - y*mu*v; + coef[0] = x*mu*mu*v*v*tau + x*em*mu*v - y*mu*mu*v*v + y*em*mu*v; + + TMath::RootsCubic(coef,roots[0],roots[1],roots[2]); + + e = roots[1]; + b = y/(tau + (em - e)/mu/v); + f = LikeMod2(mu,b,e,x,y,em,tau,v); + } + } + + return f; +} + +//_________________________________________________________________________ +Double_t MRolke::LikeMod2(Double_t mu, Double_t b, Double_t e,Int_t x,Int_t y,Double_t em,Double_t tau, Double_t v) +{ + // Profile Likelihood function for MODEL 2: + // Poisson background/Gauss Efficiency + + return 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2); +} + +//_____________________________________________________________________ + +Double_t MRolke::EvalLikeMod3(Double_t mu, Int_t x, Double_t bm, Double_t em, Double_t e, Double_t sde, Double_t sdb, Double_t b, Int_t what) +{ + // Calculates the Profile Likelihood for MODEL 3: + // Gauss background/ Gauss Efficiency + // what = 1: Maximum likelihood estimate is returned + // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned. + // what = 3: Profile Likelihood of Test hypothesis is returned + // otherwise parameters as described in the beginning of the class) + + Double_t f = 0.; + Double_t v = sde*sde; + Double_t u = sdb*sdb; + + if (what == 1) { + f = (x-bm)/em; + } + + + if (what == 2) { + mu = (x-bm)/em; + b = bm; + e = em; + f = LikeMod3(mu,b,e,x,bm,em,u,v); + } + + + if(what == 3) { + if(mu == 0.0){ + b = ((bm-u)+TMath::Sqrt((bm-u)*(bm-u)+4*x*u))/2.; + e = em; + f = LikeMod3(mu,b,e,x,bm,em,u,v); + } else { + Double_t temp[3]; + temp[0] = mu*mu*v+u; + temp[1] = mu*mu*mu*v*v+mu*v*u-mu*mu*v*em+mu*v*bm-2*u*em; + temp[2] = mu*mu*v*v*bm-mu*v*u*em-mu*v*bm*em+u*em*em-mu*mu*v*v*x; + e = (-temp[1]+TMath::Sqrt(temp[1]*temp[1]-4*temp[0]*temp[2]))/2/temp[0]; + b = bm-(u*(em-e))/v/mu; + f = LikeMod3(mu,b,e,x,bm,em,u,v); + } + } + + return f; +} + +//____________________________________________________________________ +Double_t MRolke::LikeMod3(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t bm,Double_t em,Double_t u,Double_t v) +{ + // Profile Likelihood function for MODEL 3: + // Gauss background/Gauss Efficiency + + return 2*(x * TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)-1.837877-TMath::Log(u)/2-(bm-b)*(bm-b)/u/2-TMath::Log(v)/2-(em-e)*(em-e)/v/2); +} + +//____________________________________________________________________ +Double_t MRolke::EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Double_t b, Int_t what) +{ + // Calculates the Profile Likelihood for MODEL 4: + // Poiss background/Efficiency known + // what = 1: Maximum likelihood estimate is returned + // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned. + // what = 3: Profile Likelihood of Test hypothesis is returned + // otherwise parameters as described in the beginning of the class) + + Double_t f = 0.0; + + if (what == 1) f = x-y/tau; + if (what == 2) { + mu = x-y/tau; + b = Double_t(y)/tau; + f = LikeMod4(mu,b,x,y,tau); + } + if (what == 3) { + if (mu == 0.0) { + b = Double_t(x+y)/(1+tau); + f = LikeMod4(mu,b,x,y,tau); + } else { + b = (x+y-(1+tau)*mu+sqrt((x+y-(1+tau)*mu)*(x+y-(1+tau)*mu)+4*(1+tau)*y*mu))/2/(1+tau); + f = LikeMod4(mu,b,x,y,tau); + } + } + return f; +} + +//___________________________________________________________________ +Double_t MRolke::LikeMod4(Double_t mu,Double_t b,Int_t x,Int_t y,Double_t tau) +{ + // Profile Likelihood function for MODEL 4: + // Poiss background/Efficiency known + + return 2*(x*TMath::Log(mu+b)-(mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) ); +} + +//___________________________________________________________________ +Double_t MRolke::EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Double_t b, Int_t what) +{ + // Calculates the Profile Likelihood for MODEL 5: + // Gauss background/Efficiency known + // what = 1: Maximum likelihood estimate is returned + // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned. + // what = 3: Profile Likelihood of Test hypothesis is returned + // otherwise parameters as described in the beginning of the class) + + Double_t u=sdb*sdb; + Double_t f = 0; + + if(what == 1) { + f = x - bm; + } + if(what == 2) { + mu = x-bm; + b = bm; + f = LikeMod5(mu,b,x,bm,u); + } + + if (what == 3) { + b = ((bm-u-mu)+TMath::Sqrt((bm-u-mu)*(bm-u-mu)-4*(mu*u-mu*bm-u*x)))/2; + f = LikeMod5(mu,b,x,bm,u); + } + return f; +} + +//_______________________________________________________________________ +Double_t MRolke::LikeMod5(Double_t mu,Double_t b,Int_t x,Double_t bm,Double_t u) +{ + // Profile Likelihood function for MODEL 5: + // Gauss background/Efficiency known + + return 2*(x*TMath::Log(mu+b)-(mu + b)-LogFactorial(x)-0.9189385-TMath::Log(u)/2-((bm-b)*(bm-b))/u/2); +} + +//_______________________________________________________________________ +Double_t MRolke::EvalLikeMod6(Double_t mu, Int_t x, Int_t z, Double_t e, Double_t b, Int_t m, Int_t what) +{ + // Calculates the Profile Likelihood for MODEL 6: + // Gauss known/Efficiency binomial + // what = 1: Maximum likelihood estimate is returned + // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned. + // what = 3: Profile Likelihood of Test hypothesis is returned + // otherwise parameters as described in the beginning of the class) + + Double_t coef[4],roots[3]; + Double_t f = 0.; + Double_t zm = Double_t(z)/m; + + if(what==1){ + f = (x-b)/zm; + } + + if(what==2){ + mu = (x-b)/zm; + e = zm; + f = LikeMod6(mu,b,e,x,z,m); + } + if(what == 3){ + if(mu==0){ + e = zm; + } else { + coef[3] = mu*mu; + coef[2] = mu * b - mu * x - mu*mu - mu * m; + coef[1] = mu * x - mu * b + mu * z - m * b; + coef[0] = b * z; + TMath::RootsCubic(coef,roots[0],roots[1],roots[2]); + e = roots[1]; + } + f =LikeMod6(mu,b,e,x,z,m); + } + return f; +} + +//_______________________________________________________________________ +Double_t MRolke::LikeMod6(Double_t mu,Double_t b,Double_t e,Int_t x,Int_t z,Int_t m) +{ + // Profile Likelihood function for MODEL 6: + // background known/ Efficiency binomial + + Double_t f = 0.0; + + if (z > 100 || m > 100) { + f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+(m*TMath::Log(m) - m)-(z*TMath::Log(z) - z) - ((m-z)*TMath::Log(m-z) - m + z)+z*TMath::Log(e)+(m-z)*TMath::Log(1-e)); + } else { + f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+TMath::Log(TMath::Factorial(m))-TMath::Log(TMath::Factorial(z))-TMath::Log(TMath::Factorial(m-z))+z*TMath::Log(e)+(m-z)*TMath::Log(1-e)); + } + return f; +} + + +//___________________________________________________________________________ +Double_t MRolke::EvalLikeMod7(Double_t mu, Int_t x, Double_t em, Double_t e, Double_t sde, Double_t b, Int_t what) +{ + // Calculates the Profile Likelihood for MODEL 7: + // background known/Efficiency Gauss + // what = 1: Maximum likelihood estimate is returned + // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned. + // what = 3: Profile Likelihood of Test hypothesis is returned + // otherwise parameters as described in the beginning of the class) + + Double_t v=sde*sde; + Double_t f = 0.; + + if(what == 1) { + f = (x-b)/em; + } + + if(what == 2) { + mu = (x-b)/em; + e = em; + f = LikeMod7(mu, b, e, x, em, v); + } + + if(what == 3) { + if(mu==0) { + e = em; + } else { + e = ( -(mu*em-b-mu*mu*v)-TMath::Sqrt((mu*em-b-mu*mu*v)*(mu*em-b-mu*mu*v)+4*mu*(x*mu*v-mu*b*v + b * em)))/( - mu)/2; + } + f = LikeMod7(mu, b, e, x, em, v); + } + + return f; +} + +//___________________________________________________________________________ +Double_t MRolke::LikeMod7(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t em,Double_t v) +{ + // Profile Likelihood function for MODEL 6: + // background known/ Efficiency binomial + + return 2*(x*TMath::Log(e*mu+b)-(e*mu + b)-LogFactorial(x)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2); +} + +//______________________________________________________________________ +Double_t MRolke::EvalPolynomial(Double_t x, const Int_t coef[], Int_t N) +{ + // evaluate polynomial + + const Int_t *p; + p = coef; + Double_t ans = *p++; + Int_t i = N; + + do + ans = ans * x + *p++; + while( --i ); + + return ans; +} + +//______________________________________________________________________ +Double_t MRolke::EvalMonomial(Double_t x, const Int_t coef[], Int_t N) +{ + // evaluate mononomial + + Double_t ans; + const Int_t *p; + + p = coef; + ans = x + *p++; + Int_t i = N-1; + + do + ans = ans * x + *p++; + while( --i ); + + return ans; +} +//-------------------------------------------------------------------------- +// +// Uses Stirling-Wells formula: ln(N!) ~ N*ln(N) - N + 0.5*ln(2piN) +// if N exceeds 70, otherwise use the TMath functions +// +// +Double_t MRolke::LogFactorial(Int_t n) +{ + if (n>69) + return n*TMath::Log(n)-n + 0.5*TMath::Log(TMath::TwoPi()*n); + else + return TMath::Log(TMath::Factorial(n)); +} Index: trunk/MagicSoft/Mars/mtools/MRolke.h =================================================================== --- trunk/MagicSoft/Mars/mtools/MRolke.h (revision 8087) +++ trunk/MagicSoft/Mars/mtools/MRolke.h (revision 8087) @@ -0,0 +1,91 @@ +// @(#)root/physics:$Name: not supported by cvs2svn $:$Id: MRolke.h,v 1.1 2006-10-17 12:52:15 meyer Exp $ +// Author: Jan Conrad 9/2/2004 + +/************************************************************************* + * Copyright (C) 1995-2004, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +#ifndef MARS_MRolke +#define MARS_MRolke + +#ifndef ROOT_TObject +#include "TObject.h" +#endif + +class MRolke : public TObject { + +protected: + Double_t fCL; // confidence level as a fraction [e.g. 90% = 0.9] + Double_t fUpperLimit; // the calculated upper limit + Double_t fLowerLimit; // the calculated lower limit + Int_t fSwitch; // 0: for unbounded likelihood + // 1: for bounded likelihood + + // The Calculator + + Double_t Interval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b,Int_t m); + + // LIKELIHOOD ROUTINE + + Double_t Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m, Int_t what); + + //MODEL 1 + Double_t EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t e, Double_t tau, Double_t b, Int_t m, Int_t what); + Double_t LikeMod1(Double_t mu,Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m); + void ProfLikeMod1(Double_t mu,Double_t &b, Double_t &e,Int_t x,Int_t y, Int_t z,Double_t tau,Int_t m); + Double_t LikeGradMod1(Double_t e, Double_t mu, Int_t x,Int_t y,Int_t z,Double_t tau,Int_t m); + + //MODEL 2 + Double_t EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t e,Double_t sde, Double_t tau, Double_t b, Int_t what); + + Double_t LikeMod2(Double_t mu, Double_t b, Double_t e,Int_t x,Int_t y,Double_t em,Double_t tau, Double_t v); + + //MODEL 3 + Double_t EvalLikeMod3(Double_t mu, Int_t x, Double_t bm, Double_t em, Double_t e, Double_t sde, Double_t sdb, Double_t b, Int_t what); + Double_t LikeMod3(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t bm,Double_t em,Double_t u,Double_t v); + + //MODEL 4 + Double_t EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Double_t b, Int_t what); + Double_t LikeMod4(Double_t mu,Double_t b,Int_t x,Int_t y,Double_t tau); + + //MODEL 5 + Double_t EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Double_t b, Int_t what); + Double_t LikeMod5(Double_t mu,Double_t b,Int_t x,Double_t bm,Double_t u); + + //MODEL 6 + Double_t EvalLikeMod6(Double_t mu, Int_t x, Int_t z, Double_t e, Double_t b, Int_t m, Int_t what); + Double_t LikeMod6(Double_t mu,Double_t b,Double_t e,Int_t x,Int_t z,Int_t m); + + //MODEL 7 + Double_t EvalLikeMod7(Double_t mu, Int_t x, Double_t em, Double_t e, Double_t sde, Double_t b, Int_t what); + Double_t LikeMod7(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t em,Double_t v); + + //MISC + static Double_t EvalPolynomial(Double_t x, const Int_t coef[], Int_t N); + static Double_t EvalMonomial (Double_t x, const Int_t coef[], Int_t N); + + Double_t LogFactorial(Int_t n); + + +public: + + MRolke(Double_t CL=0.9, Option_t *option = ""); + + virtual ~MRolke(); + + Double_t CalculateInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b,Int_t m); + Double_t GetUpperLimit(void) const { return fUpperLimit;} + Double_t GetLowerLimit(void) const { return fLowerLimit;} + Int_t GetSwitch(void) const { return fSwitch;} + void SetSwitch(Int_t sw) { fSwitch = sw; } + Double_t GetCL(void) const { return fCL;} + void SetCL(Double_t CL) { fCL = CL; } + + ClassDef(MRolke,1) //calculate confidence limits using the Rolke method +}; +#endif + Index: trunk/MagicSoft/Mars/mtools/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtools/Makefile (revision 8085) +++ trunk/MagicSoft/Mars/mtools/Makefile (revision 8087) @@ -37,5 +37,6 @@ MagicDomino.cc \ MagicCivilization.cc \ - MineSweeper.cc + MineSweeper.cc \ + MRolke.cc ############################################################ Index: trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h (revision 8085) +++ trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h (revision 8087) @@ -24,3 +24,5 @@ #pragma link C++ class MagicCivilization+; +#pragma link C++ class MRolke+; + #endif