Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 <cmath>
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 std::endl;
57
59
60////////////////////////////////////////////////////////////////////////////////
61
63 RooAbsReal &_lambda)
64 : RooAbsPdf(name, title),
65 x("x", "x", this, _x),
66 k("k", "k", this, _k),
67 lambda("lambda", "lambda", this, _lambda),
68 fErrorTol(1E-3),
69 fMaxIters(10),
70 fForceSum(false),
71 fHasIssuedConvWarning(false),
72 fHasIssuedSumWarning(false)
73{
74 ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
75 "MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
76}
77
78////////////////////////////////////////////////////////////////////////////////
79
81 : RooAbsPdf(other, name),
82 x("x", this, other.x),
83 k("k", this, other.k),
84 lambda("lambda", this, other.lambda),
85 fErrorTol(other.fErrorTol),
86 fMaxIters(other.fMaxIters),
87 fForceSum(other.fForceSum),
88 fHasIssuedConvWarning(false),
89 fHasIssuedSumWarning(false)
90{
91 ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
92 "MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
93}
94
95////////////////////////////////////////////////////////////////////////////////
96
98 fForceSum = flag;
99}
100
101////////////////////////////////////////////////////////////////////////////////
102
104{
105 // ENTER EXPRESSION IN TERMS OF VARIABLE ARGUMENTS HERE
106
107
108 // chi^2(0,k) gives inf and causes various problems
109 // truncate
110 double xmin = x.min();
111 double xmax = x.max();
112 double _x = x;
113 if(_x<=0){
114 // options for dealing with this
115 // return 0; // gives a funny dip
116 // _x = 1./RooNumber::infinity(); // too tall
117 _x = xmin + 1e-3*(xmax-xmin); // very small fraction of range
118 }
119
120 // special case (form below doesn't work when lambda==0)
121 if(lambda==0){
122 return ROOT::Math::chisquared_pdf(_x,k);
123 }
124
125 // three forms
126 // FIRST FORM
127 // \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
128 // could truncate sum
129
130 if ( fForceSum ){
132 coutI(InputArguments) << "RooNonCentralChiSquare sum being forced" << endl ;
134 }
135 double sum = 0;
136 double ithTerm = 0;
137 double errorTol = fErrorTol;
138 int MaxIters = fMaxIters;
139 int iDominant = (int) std::floor(lambda/2);
140 // cout <<"iDominant: " << iDominant << endl;
141
142 // do 0th term last
143 // if(iDominant==0) iDominant = 1;
144 for(int i = iDominant; ; ++i){
145 ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
146 sum+=ithTerm;
147 // cout <<"progress: " << i << " " << ithTerm/sum << endl;
148 if(ithTerm/sum < errorTol)
149 break;
150
151 if( i>iDominant+MaxIters){
154 coutW(Eval) << "RooNonCentralChiSquare did not converge: for x=" << x <<" k="<<k
155 << ", lambda="<<lambda << " fractional error = " << ithTerm/sum
156 << "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
157 << endl;
158 }
159 break;
160 }
161 }
162
163 for(int i = iDominant - 1; i >= 0; --i){
164 // cout <<"Progress: " << i << " " << ithTerm/sum << endl;
165 ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
166 sum+=ithTerm;
167 }
168
169
170 return sum;
171 }
172
173 // SECOND FORM (use MathMore function based on Bessel function (if k>2) or
174 // or regularized confluent hypergeometric limit function.
176}
177
178////////////////////////////////////////////////////////////////////////////////
179
180Int_t RooNonCentralChiSquare::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
181{
182 if (matchArgs(allVars,analVars,x)) return 1 ;
183 return 0 ;
184}
185
186////////////////////////////////////////////////////////////////////////////////
187
188double RooNonCentralChiSquare::analyticalIntegral(Int_t code, const char* rangeName) const
189{
190 R__ASSERT(code==1 );
191 // cout << "evaluating analytic integral" << endl;
192 double xmin = x.min(rangeName);
193 double xmax = x.max(rangeName);
194
195 // if xmin~0 and xmax big, then can return 1. b/c evaluate is normalized.
196
197 // special case (form below doesn't work when lambda==0)
198 if(lambda==0){
200 }
201
202 // three forms
203 // FIRST FORM
204 // \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
205 // could truncate sum
206
207
208 double sum = 0;
209 double ithTerm = 0;
210 double errorTol = fErrorTol; // for normalization allow slightly larger error
211 int MaxIters = fMaxIters; // for normalization use more terms
212
213 int iDominant = (int) std::floor(lambda/2);
214 // cout <<"iDominant: " << iDominant << endl;
215 // iDominant=0;
216 for(int i = iDominant; ; ++i){
217 ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
220 sum+=ithTerm;
221 // cout <<"progress: " << i << " " << ithTerm << " " << sum << endl;
222 if(ithTerm/sum < errorTol)
223 break;
224
225 if( i>iDominant+MaxIters){
228 coutW(Eval) << "RooNonCentralChiSquare Normalization did not converge: for k="<<k
229 << ", lambda="<<lambda << " fractional error = " << ithTerm/sum
230 << "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
231 << endl;
232 }
233 break;
234 }
235 }
236
237 for(int i = iDominant - 1; i >= 0; --i){
238 ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
241 sum+=ithTerm;
242 }
243 return sum;
244}
#define e(i)
Definition RSha256.hxx:103
#define coutI(a)
#define coutW(a)
#define ccoutD(a)
#define ClassImp(name)
Definition Rtypes.h:382
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
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:24
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=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
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).
Double_t x[n]
Definition legend1.C:17
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