ROOT   Reference Guide
RooGamma.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2 * Project: RooFit *
3 * *
4 * Simple Gamma distribution
5 * authors: Stefan A. Schmitz, Gregory Schott
6 * *
7 *****************************************************************************/
8
9/** \class RooGamma
10 \ingroup Roofit
11
12Implementation of the Gamma PDF for RooFit/RooStats.
13\f[
14f(x) = \frac{(x-\mu)^{\gamma-1} \cdot \exp^{(-(x-mu) / \beta}}{\Gamma(\gamma) \cdot \beta^{\gamma}}
15\f]
16defined for \f$x \geq 0 \f$ if \f$\mu = 0 \f$
17
18Notes from Kyle Cranmer:
19
20Wikipedia and several sources refer to the Gamma distribution as
21
22\f[
23G(\mu,\alpha,\beta) = \beta^\alpha \mu^{(\alpha-1)} \frac{e^{(-\beta \mu)}}{\Gamma(\alpha)}
24\f]
25
26Below is the correspondence:
27
28| Wikipedia | This Implementation |
29|-----------------|--------------------------|
30| \f$\alpha \f$ | \f$\gamma \f$ |
31| \f$\beta \f$ | \f$\frac{1}{\beta} \f$ |
32| \f$\mu \f$ | x |
33| 0 | \f$\mu \f$ |
34
35
36Note, that for a model Pois(N|mu), a uniform prior on mu, and a measurement N
37the posterior is in the Wikipedia parameterization Gamma(mu, alpha=N+1, beta=1)
38thus for this implementation it is:
39
40RooGamma(_x=mu,_gamma=N+1,_beta=1,_mu=0)
41
42Also note, that in this case it is equivalent to
43RooPoison(N,mu) and treating the function as a PDF in mu.
44**/
45
46#include "RooGamma.h"
47
48#include "RooRandom.h"
49#include "RooHelpers.h"
50#include "RooBatchCompute.h"
51
52#include "TMath.h"
54
55#include <cmath>
56
58
59////////////////////////////////////////////////////////////////////////////////
60
61RooGamma::RooGamma(const char *name, const char *title,
62 RooAbsReal& _x, RooAbsReal& _gamma,
63 RooAbsReal& _beta, RooAbsReal& _mu) :
64 RooAbsPdf(name,title),
65 x("x","Observable",this,_x),
66 gamma("gamma","Mean",this,_gamma),
67 beta("beta","Width",this,_beta),
68 mu("mu","Para",this,_mu)
69{
70 RooHelpers::checkRangeOfParameters(this, {&_gamma, &_beta}, 0.);
71}
72
73////////////////////////////////////////////////////////////////////////////////
74
75RooGamma::RooGamma(const RooGamma& other, const char* name) :
76 RooAbsPdf(other,name), x("x",this,other.x), gamma("mean",this,other.gamma),
77 beta("beta",this,other.beta), mu("mu",this,other.mu)
78{
79}
80
81////////////////////////////////////////////////////////////////////////////////
82
83double RooGamma::evaluate() const
84{
85 return TMath::GammaDist(x, gamma, mu, beta) ;
86}
87
88////////////////////////////////////////////////////////////////////////////////
89/// Compute multiple values of Gamma PDF.
90void RooGamma::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
91{
93 dispatch->compute(stream, RooBatchCompute::Gamma, output, nEvents,
94 {dataMap.at(x), dataMap.at(gamma), dataMap.at(beta), dataMap.at(mu)});
95}
96
97////////////////////////////////////////////////////////////////////////////////
98
99Int_t RooGamma::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
100{
101 if (matchArgs(allVars,analVars,x)) return 1 ;
102 return 0 ;
103}
104
105////////////////////////////////////////////////////////////////////////////////
106
107double RooGamma::analyticalIntegral(Int_t code, const char* rangeName) const
108{
109 R__ASSERT(code==1) ;
110
111 //integral of the Gamma distribution via ROOT::Math
112 double integral = ROOT::Math::gamma_cdf(x.max(rangeName), gamma, beta, mu) - ROOT::Math::gamma_cdf(x.min(rangeName), gamma, beta, mu);
113 return integral ;
114}
115
116////////////////////////////////////////////////////////////////////////////////
117
118Int_t RooGamma::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
119{
120 if (matchArgs(directVars,generateVars,x)) return 1 ;
121 return 0 ;
122}
123
124////////////////////////////////////////////////////////////////////////////////
125/// algorithm adapted from code example in:
126/// Marsaglia, G. and Tsang, W. W.
127/// A Simple Method for Generating Gamma Variables
128/// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
129///
130/// The speed of this algorithm depends on the speed of generating normal variates.
131/// The algorithm is limited to \f$\gamma \geq 0 \f$ !
132
134{
135 R__ASSERT(code==1) ;
136
137
138 while(1) {
139
140 double d = 0;
141 double c = 0;
142 double xgen = 0;
143 double v = 0;
144 double u = 0;
145 d = gamma -1./3.; c = 1./TMath::Sqrt(9.*d);
146
147 while(v <= 0.){
148 xgen = RooRandom::randomGenerator()->Gaus(); v = 1. + c*xgen;
149 }
151 if( u < 1.-.0331*(xgen*xgen)*(xgen*xgen) ) {
152 if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
153 x = ((d*v)* beta + mu) ;
154 break;
155 }
156 }
157 if( TMath::Log(u) < 0.5*xgen*xgen + d*(1.-v + TMath::Log(v)) ) {
158 if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
159 x = ((d*v)* beta + mu) ;
160 break;
161 }
162 }
163
164 }
165
166
167 return;
168}
169
170
#define d(i)
Definition: RSha256.hxx:102
#define c(i)
Definition: RSha256.hxx:101
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
char name[80]
Definition: TGX11.cxx:110
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
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
Implementation of the Gamma PDF for RooFit/RooStats.
Definition: RooGamma.h:20
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooGamma.cxx:99
RooRealProxy beta
Definition: RooGamma.h:39
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooGamma.cxx:118
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of Gamma PDF.
Definition: RooGamma.cxx:90
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooGamma.cxx:107
RooGamma()
Definition: RooGamma.h:22
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooGamma.cxx:83
RooRealProxy x
Definition: RooGamma.h:37
void generateEvent(Int_t code) override
algorithm adapted from code example in: Marsaglia, G.
Definition: RooGamma.cxx:133
RooRealProxy gamma
Definition: RooGamma.h:38
RooRealProxy mu
Definition: RooGamma.h:40
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:51
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.
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:274
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:672
double gamma_cdf(double x, double alpha, double theta, double x0=0)
Cumulative distribution function of the gamma distribution (lower tail).
double beta(double x, double y)
Calculates the beta function.
Double_t x[n]
Definition: legend1.C:17
double gamma(double x)
R__EXTERN RooBatchComputeInterface * dispatchCUDA
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
void checkRangeOfParameters(const RooAbsReal *callingClass, std::initializer_list< const RooAbsReal * > pars, double min=-std::numeric_limits< double >::max(), double max=std::numeric_limits< double >::max(), bool limitsInAllowedRange=false, std::string const &extraMessage="")
Check if the parameters have a range, and warn if the range extends below / above the set limits.
Definition: RooHelpers.cxx:119
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition: TMath.h:753
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition: TMath.h:659
Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1)
Computes the density function of Gamma distribution at point x.
Definition: TMath.cxx:2324
static void output()