ROOT   Reference Guide
RooChiSquarePdf.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * Kyle Cranmer
7 * *
8 *****************************************************************************/
9
10/** \class RooChiSquarePdf
11 \ingroup Roofit
12
13The PDF of the Chi Square distribution for n degrees of freedom.
14Oddly, this is hard to find in ROOT (except via relation to GammaDist).
15Here we also implement the analytic integral.
16**/
17
18#include "RooChiSquarePdf.h"
19#include "RooFit.h"
20#include "RooAbsReal.h"
21#include "RooRealVar.h"
22#include "BatchHelpers.h"
24
25#include "TMath.h"
26
27#include <cmath>
28using namespace std;
29
31
32////////////////////////////////////////////////////////////////////////////////
33
35{
36}
37
38////////////////////////////////////////////////////////////////////////////////
39
40RooChiSquarePdf::RooChiSquarePdf(const char* name, const char* title,
41 RooAbsReal& x, RooAbsReal& ndof):
42 RooAbsPdf(name, title),
43 _x("x", "Dependent", this, x),
44 _ndof("ndof","ndof", this, ndof)
45{
46}
47
48////////////////////////////////////////////////////////////////////////////////
49
51 RooAbsPdf(other, name),
52 _x("x", this, other._x),
53 _ndof("ndof",this,other._ndof)
54{
55}
56
57////////////////////////////////////////////////////////////////////////////////
58
60{
61 if(_x <= 0) return 0;
62
63 return pow(_x,(_ndof/2.)-1.) * exp(-_x/2.) / TMath::Gamma(_ndof/2.) / pow(2.,_ndof/2.);
64}
65
66////////////////////////////////////////////////////////////////////////////////
67
68namespace {
69//Author: Emmanouil Michalainas, CERN 28 Aug 2019
70
71template<class T_x, class T_ndof>
72void compute( size_t batchSize,
73 double * __restrict output,
74 T_x X, T_ndof N)
75{
76 if ( N.isBatch() ) {
77 for (size_t i=0; i<batchSize; i++) {
78 if (X[i] > 0) {
79 output[i] = 1/std::tgamma(N[i]/2.0);
80 }
81 }
82 }
83 else {
84 // N is just a scalar so bracket adapter ignores index.
85 const double gamma = 1/std::tgamma(N[2019]/2.0);
86 for (size_t i=0; i<batchSize; i++) {
87 output[i] = gamma;
88 }
89 }
90
91 constexpr double ln2 = 0.693147180559945309417232121458;
92 const double lnx0 = std::log(X[0]);
93 for (size_t i=0; i<batchSize; i++) {
94 double lnx;
95 if ( X.isBatch() ) lnx = _rf_fast_log(X[i]);
96 else lnx = lnx0;
97
98 double arg = (N[i]-2)*lnx -X[i] -N[i]*ln2;
99 output[i] *= _rf_fast_exp(0.5*arg);
100 }
101}
102};
103
104RooSpan<double> RooChiSquarePdf::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
105 using namespace BatchHelpers;
106 auto _xData = _x.getValBatch(begin, batchSize);
107 auto _ndofData = _ndof.getValBatch(begin, batchSize);
108 const bool batch_x = !_xData.empty();
109 const bool batch_ndof = !_ndofData.empty();
110
111 if (!batch_x && !batch_ndof) {
112 return {};
113 }
114 batchSize = findSize({ _xData, _ndofData });
115 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
116
117 if (batch_x && !batch_ndof ) {
119 }
120 else if (!batch_x && batch_ndof ) {
122 }
123 else if (batch_x && batch_ndof ) {
124 compute(batchSize, output.data(), _xData, _ndofData);
125 }
126 return output;
127}
128
129
130////////////////////////////////////////////////////////////////////////////////
131/// No analytical calculation available (yet) of integrals over subranges
132
133Int_t RooChiSquarePdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
134{
135 if (rangeName && strlen(rangeName)) {
136 return 0 ;
137 }
138
139 if (matchArgs(allVars, analVars, _x)) return 1;
140 return 0;
141}
142
143////////////////////////////////////////////////////////////////////////////////
144
145Double_t RooChiSquarePdf::analyticalIntegral(Int_t code, const char* rangeName) const
146{
147 assert(1 == code); (void)code;
148 Double_t xmin = _x.min(rangeName); Double_t xmax = _x.max(rangeName);
149
150 // TMath::Prob needs ndof to be an integer, or it returns 0.
151 // return TMath::Prob(xmin, _ndof) - TMath::Prob(xmax,_ndof);
152
153 // cumulative is known based on lower incomplete gamma function, or regularized gamma function
154 // Wikipedia defines lower incomplete gamma function without the normalization 1/Gamma(ndof),
155 // but it is included in the ROOT implementation.
156 Double_t pmin = TMath::Gamma(_ndof/2,xmin/2);
157 Double_t pmax = TMath::Gamma(_ndof/2,xmax/2);
158
159 // only use this if range is appropriate
160 return pmax-pmin;
161}
double _rf_fast_exp(double x)
double _rf_fast_log(double x)
#define ClassImp(name)
Definition: Rtypes.h:361
#define N
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
double pow(double, double)
double exp(double)
double log(double)
typedef void((*Func_t)())
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize, const RooArgSet *const normSet=nullptr, Tag_t ownerTag=kUnspecified)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:118
Little adapter that gives a bracket operator to types that don't have one.
Definition: BatchHelpers.h:58
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:450
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
The PDF of the Chi Square distribution for n degrees of freedom.
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
RooRealProxy _ndof
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
No analytical calculation available (yet) of integrals over subranges.
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
RooRealProxy _x
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
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.
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
double tgamma(double x)
The gamma function is defined to be the extension of the factorial to real numbers.
Double_t x[n]
Definition: legend1.C:17
size_t findSize(std::vector< RooSpan< const double > > parameters)
This function returns the minimum size of the non-zero-sized batches.
double gamma(double x)
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:348
static void output(int code)
Definition: gifencode.c:226