Logo ROOT  
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
13The PDF of the Non-Central Chi Square distribution for n degrees of freedom.
14It is the asymptotic distribution of the profile likelihood ratio test q_mu
15when a different mu' is true. It is Wald's generalization of Wilks' Theorem.
16
17See:
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
26It requires MathMore to evaluate for non-integer degrees of freedom, k.
27
28When the Mathmore library is available we can use the MathMore libraries implemented using GSL.
29It makes use of the modified Bessel function of the first kind (for k > 2). For k < 2 it uses
30the hypergeometric function 0F1.
31When is not available we use explicit summation of normal chi-squared distributions
32The usage of the sum can be forced by calling SetForceSum(true);
33
34This implementation could be improved. BOOST has a nice implementation:
35
36http://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
38http://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
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
56using namespace std;
57
59
60////////////////////////////////////////////////////////////////////////////////
61
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 xmin = x.min();
129 double 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
204Int_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
212double RooNonCentralChiSquare::analyticalIntegral(Int_t code, const char* rangeName) const
213{
214 R__ASSERT(code==1 );
215 // cout << "evaluating analytic integral" << endl;
216 double xmin = x.min(rangeName);
217 double 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)
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)
265 sum+=ithTerm;
266 }
267 return sum;
268}
#define e(i)
Definition: RSha256.hxx:103
#define coutI(a)
Definition: RooMsgService.h:34
#define coutW(a)
Definition: RooMsgService.h:36
#define coutF(a)
Definition: RooMsgService.h:38
#define ccoutD(a)
Definition: RooMsgService.h:41
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
The PDF of the Non-Central Chi Square distribution for n degrees of freedom.
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
double chisquared_pdf(double x, double r, double x0=0)
Probability density function of the distribution with degrees of freedom.
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...
double chisquared_cdf(double x, double r, double x0=0)
Cumulative distribution function of the distribution with degrees of freedom (lower tail).
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1753
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
Double_t x[n]
Definition: legend1.C:17
@ InputArguments
Definition: RooGlobalFunc.h:64
Double_t Floor(Double_t x)
Rounds x downward, returning the largest integral value that is not greater than x.
Definition: TMath.h:677
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:93
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:353
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345