Logo ROOT   6.10/09
Reference Guide
RooNonCentralChiSquare.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitModels *
4  * @(#)root/roofit:$Id: RooNonCentralChiSquare *
5  * Authors: *
6  * Kyle Cranmer
7  * *
8  *****************************************************************************/
9 
10 /** \class RooNonCentralChiSquare
11  \ingroup Roofit
12 
13 The PDF of the Non-Central Chi Square distribution for n degrees of freedom.
14 It is the asymptotic distribution of the profile likelihood ratio test q_mu
15 when a different mu' is true. It is Wald's generalization of Wilks' Theorem.
16 
17 See:
18 
19  Asymptotic formulae for likelihood-based tests of new physics
20 
21  By Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells
22  http://arXiv.org/abs/arXiv:1007.1727
23 
24  [Wikipedia](http://en.wikipedia.org/wiki/Noncentral_chi-square_distribution#Approximation)
25 
26 It requires MathMore to evaluate for non-integer degrees of freedom, k.
27 
28 When the Mathmore library is available we can use the MathMore libraries implemented using GSL.
29 It makes use of the modified Bessel function of the first kind (for k > 2). For k < 2 it uses
30 the hypergeometric function 0F1.
31 When is not available we use explicit summation of normal chi-squared distributions
32 The usage of the sum can be forced by calling SetForceSum(true);
33 
34 This implementation could be improved. BOOST has a nice implementation:
35 
36 http://live.boost.org/doc/libs/1_42_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/nc_chi_squared_dist.html
37 
38 http://wesnoth.repositoryhosting.com/trac/wesnoth_wesnoth/browser/trunk/include/boost/math/distributions/non_central_chi_squared.hpp?rev=6
39 **/
40 
41 #include "Riostream.h"
42 
43 #include "RooNonCentralChiSquare.h"
44 #include "RooAbsReal.h"
45 #include "RooAbsCategory.h"
46 #include <math.h>
47 #include "TMath.h"
48 //#include "RooNumber.h"
49 #include "Math/DistFunc.h"
50 
51 
52 #include "RooMsgService.h"
53 
54 #include "TError.h"
55 
56 using namespace std;
57 
59 
60 ////////////////////////////////////////////////////////////////////////////////
61 
62 RooNonCentralChiSquare::RooNonCentralChiSquare(const char *name, const char *title,
63  RooAbsReal& _x,
64  RooAbsReal& _k,
65  RooAbsReal& _lambda) :
66  RooAbsPdf(name,title),
67  x("x","x",this,_x),
68  k("k","k",this,_k),
69  lambda("lambda","lambda",this,_lambda),
70  fErrorTol(1E-3),
71  fMaxIters(10),
72  fHasIssuedConvWarning(false),
73  fHasIssuedSumWarning(false)
74 {
75 #ifdef R__HAS_MATHMORE
76  ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
77  "MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
78  fForceSum = false;
79 #else
80  fForceSum = true;
81 #endif
82 }
83 
84 ////////////////////////////////////////////////////////////////////////////////
85 
87  RooAbsPdf(other,name),
88  x("x",this,other.x),
89  k("k",this,other.k),
90  lambda("lambda",this,other.lambda),
91  fErrorTol(other.fErrorTol),
92  fMaxIters(other.fMaxIters),
93  fHasIssuedConvWarning(false),
94  fHasIssuedSumWarning(false)
95 {
96 #ifdef R__HAS_MATHMORE
97  ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
98  "MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
99  fForceSum = other.fForceSum;
100 #else
101  fForceSum = true;
102 #endif
103 }
104 
105 ////////////////////////////////////////////////////////////////////////////////
106 
108  fForceSum = flag;
109 #ifndef R__HAS_MATHMORE
110  if (!fForceSum) {
111  ccoutD(InputArguments) << "RooNonCentralChiSquare::SetForceSum" << GetName() <<
112  "MathMore is not available- ForceSum must be on "<< endl ;
113  fForceSum = true;
114  }
115 #endif
116 
117 }
118 
119 ////////////////////////////////////////////////////////////////////////////////
120 
122 {
123  // ENTER EXPRESSION IN TERMS OF VARIABLE ARGUMENTS HERE
124 
125 
126  // chi^2(0,k) gives inf and causes various problems
127  // truncate
128  Double_t xmin = x.min();
129  Double_t xmax = x.max();
130  double _x = x;
131  if(_x<=0){
132  // options for dealing with this
133  // return 0; // gives a funny dip
134  // _x = 1./RooNumber::infinity(); // too tall
135  _x = xmin + 1e-3*(xmax-xmin); // very small fraction of range
136  }
137 
138  // special case (form below doesn't work when lambda==0)
139  if(lambda==0){
140  return ROOT::Math::chisquared_pdf(_x,k);
141  }
142 
143  // three forms
144  // FIRST FORM
145  // \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
146  // could truncate sum
147 
148  if ( fForceSum ){
150  coutI(InputArguments) << "RooNonCentralChiSquare sum being forced" << endl ;
152  }
153  double sum = 0;
154  double ithTerm = 0;
155  double errorTol = fErrorTol;
156  int MaxIters = fMaxIters;
157  int iDominant = (int) TMath::Floor(lambda/2);
158  // cout <<"iDominant: " << iDominant << endl;
159 
160  // do 0th term last
161  // if(iDominant==0) iDominant = 1;
162  for(int i = iDominant; ; ++i){
163  ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
164  sum+=ithTerm;
165  // cout <<"progress: " << i << " " << ithTerm/sum << endl;
166  if(ithTerm/sum < errorTol)
167  break;
168 
169  if( i>iDominant+MaxIters){
172  coutW(Eval) << "RooNonCentralChiSquare did not converge: for x=" << x <<" k="<<k
173  << ", lambda="<<lambda << " fractional error = " << ithTerm/sum
174  << "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
175  << endl;
176  }
177  break;
178  }
179  }
180 
181  for(int i = iDominant - 1; i >= 0; --i){
182  // cout <<"Progress: " << i << " " << ithTerm/sum << endl;
183  ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
184  sum+=ithTerm;
185  }
186 
187 
188  return sum;
189  }
190 
191  // SECOND FORM (use MathMore function based on Bessel function (if k>2) or
192  // or regularized confluent hypergeometric limit function.
193 #ifdef R__HAS_MATHMORE
195 #else
196  coutF(Eval) << "RooNonCentralChisquare: ForceSum must be set" << endl;
197  return 0;
198 #endif
199 
200 }
201 
202 ////////////////////////////////////////////////////////////////////////////////
203 
204 Int_t RooNonCentralChiSquare::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
205 {
206  if (matchArgs(allVars,analVars,x)) return 1 ;
207  return 0 ;
208 }
209 
210 ////////////////////////////////////////////////////////////////////////////////
211 
213 {
214  R__ASSERT(code==1 );
215  // cout << "evaluating analytic integral" << endl;
216  Double_t xmin = x.min(rangeName);
217  Double_t xmax = x.max(rangeName);
218 
219  // if xmin~0 and xmax big, then can return 1. b/c evaluate is normalized.
220 
221  // special case (form below doesn't work when lambda==0)
222  if(lambda==0){
224  }
225 
226  // three forms
227  // FIRST FORM
228  // \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
229  // could truncate sum
230 
231 
232  double sum = 0;
233  double ithTerm = 0;
234  double errorTol = fErrorTol; // for nomralization allow slightly larger error
235  int MaxIters = fMaxIters; // for normalization use more terms
236 
237  int iDominant = (int) TMath::Floor(lambda/2);
238  // cout <<"iDominant: " << iDominant << endl;
239  // iDominant=0;
240  for(int i = iDominant; ; ++i){
241  ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
242  *( ROOT::Math::chisquared_cdf(xmax,k+2*i)/TMath::Gamma(i+1)
243  - ROOT::Math::chisquared_cdf(xmin,k+2*i)/TMath::Gamma(i+1) );
244  sum+=ithTerm;
245  // cout <<"progress: " << i << " " << ithTerm << " " << sum << endl;
246  if(ithTerm/sum < errorTol)
247  break;
248 
249  if( i>iDominant+MaxIters){
252  coutW(Eval) << "RooNonCentralChiSquare Normalization did not converge: for k="<<k
253  << ", lambda="<<lambda << " fractional error = " << ithTerm/sum
254  << "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
255  << endl;
256  }
257  break;
258  }
259  }
260 
261  for(int i = iDominant - 1; i >= 0; --i){
262  ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
263  *( ROOT::Math::chisquared_cdf(xmax,k+2*i)/TMath::Gamma(i+1)
264  -ROOT::Math::chisquared_cdf(xmin,k+2*i)/TMath::Gamma(i+1) );
265  sum+=ithTerm;
266  }
267  return sum;
268 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:145
static long int sum(long int i)
Definition: Factory.cxx:2162
float xmin
Definition: THbookFile.cxx:93
Double_t Floor(Double_t x)
Definition: TMath.h:600
double chisquared_pdf(double x, double r, double x0=0)
Probability density function of the distribution with degrees of freedom.
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
#define coutI(a)
Definition: RooMsgService.h:31
#define R__ASSERT(e)
Definition: TError.h:96
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:352
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:33
The PDF of the Non-Central Chi Square distribution for n degrees of freedom.
Double_t x[n]
Definition: legend1.C:17
double pow(double, double)
float xmax
Definition: THbookFile.cxx:93
#define coutF(a)
Definition: RooMsgService.h:35
constexpr Double_t E()
Definition: TMath.h:74
double noncentral_chisquared_pdf(double x, double r, double lambda)
Probability density function of the non central distribution with degrees of freedom and the noon-c...
#define ClassImp(name)
Definition: Rtypes.h:336
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
double chisquared_cdf(double x, double r, double x0=0)
Cumulative distribution function of the distribution with degrees of freedom (lower tail)...
double exp(double)
#define ccoutD(a)
Definition: RooMsgService.h:37