/*****************************************************************************
* Project: RooFit                                                           *
* Package: RooFitModels                                                     *
* @(#)root/roofit:\$Id: RooNonCentralChiSquare                               *
* Authors:                                                                  *
*   Kyle Cranmer
*                                                                           *
*****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// The PDF of the Non-Central Chi Square distribution for n degrees of freedom.
// It is the asymptotic distribution of the profile likeihood ratio test q_mu
// when a different mu' is true.  It is Wald's generalization of Wilks' Theorem.
//
// See:
//  Asymptotic formulae for likelihood-based tests of new physics
//     By Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells
//     http://arXiv.org/abs/arXiv:1007.1727
//
//  Wikipedia:
//    http://en.wikipedia.org/wiki/Noncentral_chi-square_distribution#Approximation
//
// It requires MathMore to evaluate for non-integer degrees of freedom, k.
//
// When the Mathmore library is available we can use the MathMore libraries impelmented using GSL.
// It makes use of the modified Bessel function of the first kind (for k > 2). For k < 2 it uses
// the hypergeometric function 0F1.
// When is not available we use explicit summation of normal chi-squared distributions
// The usage of the sum can be forced by calling SetForceSum(true);
//
// This implementation could be improved.  BOOST has a nice implementation:
// 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
// http://wesnoth.repositoryhosting.com/trac/wesnoth_wesnoth/browser/trunk/include/boost/math/distributions/non_central_chi_squared.hpp?rev=6
// END_HTML
//

#include "Riostream.h"

#include "RooNonCentralChiSquare.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include <math.h>
#include "TMath.h"
//#include "RooNumber.h"
#include "Math/DistFunc.h"

#include "RooMsgService.h"

#include "TError.h"

using namespace std;

ClassImp(RooNonCentralChiSquare)

RooNonCentralChiSquare::RooNonCentralChiSquare(const char *name, const char *title,
RooAbsReal& _x,
RooAbsReal& _k,
RooAbsReal& _lambda) :
RooAbsPdf(name,title),
x("x","x",this,_x),
k("k","k",this,_k),
lambda("lambda","lambda",this,_lambda),
fErrorTol(1E-3),
fMaxIters(10),
fHasIssuedConvWarning(false),
fHasIssuedSumWarning(false)
{
#ifdef R__HAS_MATHMORE
ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
"MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
fForceSum = false;
#else
fForceSum = true;
#endif
}

RooNonCentralChiSquare::RooNonCentralChiSquare(const RooNonCentralChiSquare& other, const char* name) :
RooAbsPdf(other,name),
x("x",this,other.x),
k("k",this,other.k),
lambda("lambda",this,other.lambda),
fErrorTol(other.fErrorTol),
fMaxIters(other.fMaxIters),
fHasIssuedConvWarning(false),
fHasIssuedSumWarning(false)
{
#ifdef R__HAS_MATHMORE
ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
"MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
fForceSum = other.fForceSum;
#else
fForceSum = true;
#endif
}

void RooNonCentralChiSquare::SetForceSum(Bool_t flag) {
fForceSum = flag;
#ifndef R__HAS_MATHMORE
if (!fForceSum) {
ccoutD(InputArguments) << "RooNonCentralChiSquare::SetForceSum" << GetName() <<
"MathMore is not available- ForceSum must be on "<< endl ;
fForceSum = true;
}
#endif

}

Double_t RooNonCentralChiSquare::evaluate() const
{
// ENTER EXPRESSION IN TERMS OF VARIABLE ARGUMENTS HERE

// chi^2(0,k) gives inf and causes various problems
// truncate
Double_t xmin = x.min();
Double_t xmax = x.max();
double _x = x;
if(_x<=0){
// options for dealing with this
//     return 0; // gives a funny dip
//     _x = 1./RooNumber::infinity(); // too tall
_x = xmin + 1e-3*(xmax-xmin); // very small fraction of range
}

// special case (form below doesn't work when lambda==0)
if(lambda==0){
return ROOT::Math::chisquared_pdf(_x,k);
}

// three forms
// FIRST FORM
// \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
// could truncate sum

if ( fForceSum  ){
if(!fHasIssuedSumWarning){
coutI(InputArguments) << "RooNonCentralChiSquare sum being forced" << endl ;
fHasIssuedSumWarning=true;
}
double sum = 0;
double ithTerm = 0;
double errorTol = fErrorTol;
int MaxIters = fMaxIters;
int iDominant = (int) TMath::Floor(lambda/2);
//     cout <<"iDominant: " << iDominant << endl;

// do 0th term last
//     if(iDominant==0) iDominant = 1;
for(int i = iDominant; ; ++i){
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
sum+=ithTerm;
//       cout <<"progress: " << i << " " << ithTerm/sum << endl;
if(ithTerm/sum < errorTol)
break;

if( i>iDominant+MaxIters){
if(!fHasIssuedConvWarning){
fHasIssuedConvWarning=true;
coutW(Eval) << "RooNonCentralChiSquare did not converge: for x=" << x <<" k="<<k
<< ", lambda="<<lambda << " fractional error = " << ithTerm/sum
<< "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
<< endl;
}
break;
}
}

for(int i = iDominant - 1; i >= 0; --i){
//       cout <<"Progress: " << i << " " << ithTerm/sum << endl;
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
sum+=ithTerm;
}

return sum;
}

// SECOND FORM (use MathMore function based on Bessel function (if k>2) or
// or  regularized confluent hypergeometric limit function.
#ifdef R__HAS_MATHMORE
return  ROOT::Math::noncentral_chisquared_pdf(_x,k,lambda);
#else
coutF(Eval) << "RooNonCentralChisquare: ForceSum must be set" << endl;
return 0;
#endif

}

//_____________________________________________________________________________
Int_t RooNonCentralChiSquare::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
{
if (matchArgs(allVars,analVars,x)) return 1 ;
return 0 ;
}

//_____________________________________________________________________________
Double_t RooNonCentralChiSquare::analyticalIntegral(Int_t code, const char* rangeName) const
{
R__ASSERT(code==1 );
//  cout << "evaluating analytic integral" << endl;
Double_t xmin = x.min(rangeName);
Double_t xmax = x.max(rangeName);

// if xmin~0 and xmax big, then can return 1. b/c evaluate is normalized.

// special case (form below doesn't work when lambda==0)
if(lambda==0){
return (ROOT::Math::chisquared_cdf(xmax,k) - ROOT::Math::chisquared_cdf(xmin,k));
}

// three forms
// FIRST FORM
// \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
// could truncate sum

double sum = 0;
double ithTerm = 0;
double errorTol = fErrorTol; // for nomralization allow slightly larger error
int MaxIters = fMaxIters; // for normalization use more terms

int iDominant = (int) TMath::Floor(lambda/2);
//     cout <<"iDominant: " << iDominant << endl;
//   iDominant=0;
for(int i = iDominant; ; ++i){
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
*( ROOT::Math::chisquared_cdf(xmax,k+2*i)/TMath::Gamma(i+1)
- ROOT::Math::chisquared_cdf(xmin,k+2*i)/TMath::Gamma(i+1) );
sum+=ithTerm;
//     cout <<"progress: " << i << " " << ithTerm << " " << sum << endl;
if(ithTerm/sum < errorTol)
break;

if( i>iDominant+MaxIters){
if(!fHasIssuedConvWarning){
fHasIssuedConvWarning=true;
coutW(Eval) << "RooNonCentralChiSquare Normalization did not converge: for k="<<k
<< ", lambda="<<lambda << " fractional error = " << ithTerm/sum
<< "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
<< endl;
}
break;
}
}

for(int i = iDominant - 1; i >= 0; --i){
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
*( ROOT::Math::chisquared_cdf(xmax,k+2*i)/TMath::Gamma(i+1)
-ROOT::Math::chisquared_cdf(xmin,k+2*i)/TMath::Gamma(i+1) );
sum+=ithTerm;
}
return sum;
}


RooNonCentralChiSquare.cxx:1
RooNonCentralChiSquare.cxx:2
RooNonCentralChiSquare.cxx:3
RooNonCentralChiSquare.cxx:4
RooNonCentralChiSquare.cxx:5
RooNonCentralChiSquare.cxx:6
RooNonCentralChiSquare.cxx:7
RooNonCentralChiSquare.cxx:8
RooNonCentralChiSquare.cxx:9
RooNonCentralChiSquare.cxx:10
RooNonCentralChiSquare.cxx:11
RooNonCentralChiSquare.cxx:12
RooNonCentralChiSquare.cxx:13
RooNonCentralChiSquare.cxx:14
RooNonCentralChiSquare.cxx:15
RooNonCentralChiSquare.cxx:16
RooNonCentralChiSquare.cxx:17
RooNonCentralChiSquare.cxx:18
RooNonCentralChiSquare.cxx:19
RooNonCentralChiSquare.cxx:20
RooNonCentralChiSquare.cxx:21
RooNonCentralChiSquare.cxx:22
RooNonCentralChiSquare.cxx:23
RooNonCentralChiSquare.cxx:24
RooNonCentralChiSquare.cxx:25
RooNonCentralChiSquare.cxx:26
RooNonCentralChiSquare.cxx:27
RooNonCentralChiSquare.cxx:28
RooNonCentralChiSquare.cxx:29
RooNonCentralChiSquare.cxx:30
RooNonCentralChiSquare.cxx:31
RooNonCentralChiSquare.cxx:32
RooNonCentralChiSquare.cxx:33
RooNonCentralChiSquare.cxx:34
RooNonCentralChiSquare.cxx:35
RooNonCentralChiSquare.cxx:36
RooNonCentralChiSquare.cxx:37
RooNonCentralChiSquare.cxx:38
RooNonCentralChiSquare.cxx:39
RooNonCentralChiSquare.cxx:40
RooNonCentralChiSquare.cxx:41
RooNonCentralChiSquare.cxx:42
RooNonCentralChiSquare.cxx:43
RooNonCentralChiSquare.cxx:44
RooNonCentralChiSquare.cxx:45
RooNonCentralChiSquare.cxx:46
RooNonCentralChiSquare.cxx:47
RooNonCentralChiSquare.cxx:48
RooNonCentralChiSquare.cxx:49
RooNonCentralChiSquare.cxx:50
RooNonCentralChiSquare.cxx:51
RooNonCentralChiSquare.cxx:52
RooNonCentralChiSquare.cxx:53
RooNonCentralChiSquare.cxx:54
RooNonCentralChiSquare.cxx:55
RooNonCentralChiSquare.cxx:56
RooNonCentralChiSquare.cxx:57
RooNonCentralChiSquare.cxx:58
RooNonCentralChiSquare.cxx:59
RooNonCentralChiSquare.cxx:60
RooNonCentralChiSquare.cxx:61
RooNonCentralChiSquare.cxx:62
RooNonCentralChiSquare.cxx:63
RooNonCentralChiSquare.cxx:64
RooNonCentralChiSquare.cxx:65
RooNonCentralChiSquare.cxx:66
RooNonCentralChiSquare.cxx:67
RooNonCentralChiSquare.cxx:68
RooNonCentralChiSquare.cxx:69
RooNonCentralChiSquare.cxx:70
RooNonCentralChiSquare.cxx:71
RooNonCentralChiSquare.cxx:72
RooNonCentralChiSquare.cxx:73
RooNonCentralChiSquare.cxx:74
RooNonCentralChiSquare.cxx:75
RooNonCentralChiSquare.cxx:76
RooNonCentralChiSquare.cxx:77
RooNonCentralChiSquare.cxx:78
RooNonCentralChiSquare.cxx:79
RooNonCentralChiSquare.cxx:80
RooNonCentralChiSquare.cxx:81
RooNonCentralChiSquare.cxx:82
RooNonCentralChiSquare.cxx:83
RooNonCentralChiSquare.cxx:84
RooNonCentralChiSquare.cxx:85
RooNonCentralChiSquare.cxx:86
RooNonCentralChiSquare.cxx:87
RooNonCentralChiSquare.cxx:88
RooNonCentralChiSquare.cxx:89
RooNonCentralChiSquare.cxx:90
RooNonCentralChiSquare.cxx:91
RooNonCentralChiSquare.cxx:92
RooNonCentralChiSquare.cxx:93
RooNonCentralChiSquare.cxx:94
RooNonCentralChiSquare.cxx:95
RooNonCentralChiSquare.cxx:96
RooNonCentralChiSquare.cxx:97
RooNonCentralChiSquare.cxx:98
RooNonCentralChiSquare.cxx:99
RooNonCentralChiSquare.cxx:100
RooNonCentralChiSquare.cxx:101
RooNonCentralChiSquare.cxx:102
RooNonCentralChiSquare.cxx:103
RooNonCentralChiSquare.cxx:104
RooNonCentralChiSquare.cxx:105
RooNonCentralChiSquare.cxx:106
RooNonCentralChiSquare.cxx:107
RooNonCentralChiSquare.cxx:108
RooNonCentralChiSquare.cxx:109
RooNonCentralChiSquare.cxx:110
RooNonCentralChiSquare.cxx:111
RooNonCentralChiSquare.cxx:112
RooNonCentralChiSquare.cxx:113
RooNonCentralChiSquare.cxx:114
RooNonCentralChiSquare.cxx:115
RooNonCentralChiSquare.cxx:116
RooNonCentralChiSquare.cxx:117
RooNonCentralChiSquare.cxx:118
RooNonCentralChiSquare.cxx:119
RooNonCentralChiSquare.cxx:120
RooNonCentralChiSquare.cxx:121
RooNonCentralChiSquare.cxx:122
RooNonCentralChiSquare.cxx:123
RooNonCentralChiSquare.cxx:124
RooNonCentralChiSquare.cxx:125
RooNonCentralChiSquare.cxx:126
RooNonCentralChiSquare.cxx:127
RooNonCentralChiSquare.cxx:128
RooNonCentralChiSquare.cxx:129
RooNonCentralChiSquare.cxx:130
RooNonCentralChiSquare.cxx:131
RooNonCentralChiSquare.cxx:132
RooNonCentralChiSquare.cxx:133
RooNonCentralChiSquare.cxx:134
RooNonCentralChiSquare.cxx:135
RooNonCentralChiSquare.cxx:136
RooNonCentralChiSquare.cxx:137
RooNonCentralChiSquare.cxx:138
RooNonCentralChiSquare.cxx:139
RooNonCentralChiSquare.cxx:140
RooNonCentralChiSquare.cxx:141
RooNonCentralChiSquare.cxx:142
RooNonCentralChiSquare.cxx:143
RooNonCentralChiSquare.cxx:144
RooNonCentralChiSquare.cxx:145
RooNonCentralChiSquare.cxx:146
RooNonCentralChiSquare.cxx:147
RooNonCentralChiSquare.cxx:148
RooNonCentralChiSquare.cxx:149
RooNonCentralChiSquare.cxx:150
RooNonCentralChiSquare.cxx:151
RooNonCentralChiSquare.cxx:152
RooNonCentralChiSquare.cxx:153
RooNonCentralChiSquare.cxx:154
RooNonCentralChiSquare.cxx:155
RooNonCentralChiSquare.cxx:156
RooNonCentralChiSquare.cxx:157
RooNonCentralChiSquare.cxx:158
RooNonCentralChiSquare.cxx:159
RooNonCentralChiSquare.cxx:160
RooNonCentralChiSquare.cxx:161
RooNonCentralChiSquare.cxx:162
RooNonCentralChiSquare.cxx:163
RooNonCentralChiSquare.cxx:164
RooNonCentralChiSquare.cxx:165
RooNonCentralChiSquare.cxx:166
RooNonCentralChiSquare.cxx:167
RooNonCentralChiSquare.cxx:168
RooNonCentralChiSquare.cxx:169
RooNonCentralChiSquare.cxx:170
RooNonCentralChiSquare.cxx:171
RooNonCentralChiSquare.cxx:172
RooNonCentralChiSquare.cxx:173
RooNonCentralChiSquare.cxx:174
RooNonCentralChiSquare.cxx:175
RooNonCentralChiSquare.cxx:176
RooNonCentralChiSquare.cxx:177
RooNonCentralChiSquare.cxx:178
RooNonCentralChiSquare.cxx:179
RooNonCentralChiSquare.cxx:180
RooNonCentralChiSquare.cxx:181
RooNonCentralChiSquare.cxx:182
RooNonCentralChiSquare.cxx:183
RooNonCentralChiSquare.cxx:184
RooNonCentralChiSquare.cxx:185
RooNonCentralChiSquare.cxx:186
RooNonCentralChiSquare.cxx:187
RooNonCentralChiSquare.cxx:188
RooNonCentralChiSquare.cxx:189
RooNonCentralChiSquare.cxx:190
RooNonCentralChiSquare.cxx:191
RooNonCentralChiSquare.cxx:192
RooNonCentralChiSquare.cxx:193
RooNonCentralChiSquare.cxx:194
RooNonCentralChiSquare.cxx:195
RooNonCentralChiSquare.cxx:196
RooNonCentralChiSquare.cxx:197
RooNonCentralChiSquare.cxx:198
RooNonCentralChiSquare.cxx:199
RooNonCentralChiSquare.cxx:200
RooNonCentralChiSquare.cxx:201
RooNonCentralChiSquare.cxx:202
RooNonCentralChiSquare.cxx:203
RooNonCentralChiSquare.cxx:204
RooNonCentralChiSquare.cxx:205
RooNonCentralChiSquare.cxx:206
RooNonCentralChiSquare.cxx:207
RooNonCentralChiSquare.cxx:208
RooNonCentralChiSquare.cxx:209
RooNonCentralChiSquare.cxx:210
RooNonCentralChiSquare.cxx:211
RooNonCentralChiSquare.cxx:212
RooNonCentralChiSquare.cxx:213
RooNonCentralChiSquare.cxx:214
RooNonCentralChiSquare.cxx:215
RooNonCentralChiSquare.cxx:216
RooNonCentralChiSquare.cxx:217
RooNonCentralChiSquare.cxx:218
RooNonCentralChiSquare.cxx:219
RooNonCentralChiSquare.cxx:220
RooNonCentralChiSquare.cxx:221
RooNonCentralChiSquare.cxx:222
RooNonCentralChiSquare.cxx:223
RooNonCentralChiSquare.cxx:224
RooNonCentralChiSquare.cxx:225
RooNonCentralChiSquare.cxx:226
RooNonCentralChiSquare.cxx:227
RooNonCentralChiSquare.cxx:228
RooNonCentralChiSquare.cxx:229
RooNonCentralChiSquare.cxx:230
RooNonCentralChiSquare.cxx:231
RooNonCentralChiSquare.cxx:232
RooNonCentralChiSquare.cxx:233
RooNonCentralChiSquare.cxx:234
RooNonCentralChiSquare.cxx:235
RooNonCentralChiSquare.cxx:236
RooNonCentralChiSquare.cxx:237
RooNonCentralChiSquare.cxx:238
RooNonCentralChiSquare.cxx:239
RooNonCentralChiSquare.cxx:240
RooNonCentralChiSquare.cxx:241
RooNonCentralChiSquare.cxx:242
RooNonCentralChiSquare.cxx:243
RooNonCentralChiSquare.cxx:244
RooNonCentralChiSquare.cxx:245
RooNonCentralChiSquare.cxx:246
RooNonCentralChiSquare.cxx:247
RooNonCentralChiSquare.cxx:248
RooNonCentralChiSquare.cxx:249
RooNonCentralChiSquare.cxx:250
RooNonCentralChiSquare.cxx:251
RooNonCentralChiSquare.cxx:252
RooNonCentralChiSquare.cxx:253
RooNonCentralChiSquare.cxx:254
RooNonCentralChiSquare.cxx:255
RooNonCentralChiSquare.cxx:256
RooNonCentralChiSquare.cxx:257
RooNonCentralChiSquare.cxx:258
RooNonCentralChiSquare.cxx:259
RooNonCentralChiSquare.cxx:260
RooNonCentralChiSquare.cxx:261
RooNonCentralChiSquare.cxx:262
RooNonCentralChiSquare.cxx:263
RooNonCentralChiSquare.cxx:264
RooNonCentralChiSquare.cxx:265