RooNonCentralChiSquare.cxx

Go to the documentation of this file.

/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitModels                                                     *
 * @(#)root/roofit:$Id: RooNonCentralChiSquare                               *
 * Authors:                                                                  *
 *   Kyle Cranmer
 *                                                                           *
 *****************************************************************************/
 
/** \class RooNonCentralChiSquare
    \ingroup Roofit
 
The PDF of the Non-Central Chi Square distribution for n degrees of freedom.
It is the asymptotic distribution of the profile likelihood 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 implemented 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
**/
 
#include "Riostream.h"
 
#include "RooNonCentralChiSquare.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include <cmath>
#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)
{
   ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
      "MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
   fForceSum = false;
}
 
////////////////////////////////////////////////////////////////////////////////
 
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)
{
   ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
     "MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
   fForceSum = other.fForceSum;
}
 
////////////////////////////////////////////////////////////////////////////////
 
void RooNonCentralChiSquare::SetForceSum(bool flag) {
   fForceSum = flag;
}
 
////////////////////////////////////////////////////////////////////////////////
 
double RooNonCentralChiSquare::evaluate() const
{
   // ENTER EXPRESSION IN TERMS OF VARIABLE ARGUMENTS HERE
 
 
   // chi^2(0,k) gives inf and causes various problems
   // truncate
   double xmin = x.min();
   double 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.
   return  ROOT::Math::noncentral_chisquared_pdf(_x,k,lambda);
}
 
////////////////////////////////////////////////////////////////////////////////
 
Int_t RooNonCentralChiSquare::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
{
   if (matchArgs(allVars,analVars,x)) return 1 ;
   return 0 ;
}
 
////////////////////////////////////////////////////////////////////////////////
 
double RooNonCentralChiSquare::analyticalIntegral(Int_t code, const char* rangeName) const
{
   R__ASSERT(code==1 );
   //  cout << "evaluating analytic integral" << endl;
   double xmin = x.min(rangeName);
   double 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 normalization 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;
}