Logo 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
40`RooGamma(_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
84{
85 return TMath::GammaDist(x, gamma, mu, beta) ;
86}
87
88////////////////////////////////////////////////////////////////////////////////
89/// Compute multiple values of Gamma PDF.
91 return RooBatchCompute::dispatch->computeGamma(this, evalData, x->getValues(evalData, normSet), gamma->getValues(evalData, normSet), beta->getValues(evalData, normSet), mu->getValues(evalData, normSet));
92}
93
94////////////////////////////////////////////////////////////////////////////////
95
96Int_t RooGamma::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
97{
98 if (matchArgs(allVars,analVars,x)) return 1 ;
99 return 0 ;
100}
101
102////////////////////////////////////////////////////////////////////////////////
103
104Double_t RooGamma::analyticalIntegral(Int_t code, const char* rangeName) const
105{
106 R__ASSERT(code==1) ;
107
108 //integral of the Gamma distribution via ROOT::Math
109 Double_t integral = ROOT::Math::gamma_cdf(x.max(rangeName), gamma, beta, mu) - ROOT::Math::gamma_cdf(x.min(rangeName), gamma, beta, mu);
110 return integral ;
111}
112
113////////////////////////////////////////////////////////////////////////////////
114
115Int_t RooGamma::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
116{
117 if (matchArgs(directVars,generateVars,x)) return 1 ;
118 return 0 ;
119}
120
121////////////////////////////////////////////////////////////////////////////////
122/// algorithm adapted from code example in:
123/// Marsaglia, G. and Tsang, W. W.
124/// A Simple Method for Generating Gamma Variables
125/// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
126///
127/// The speed of this algorithm depends on the speed of generating normal variates.
128/// The algorithm is limited to \f$ \gamma \geq 0 \f$ !
129
131{
132 R__ASSERT(code==1) ;
133
134
135 while(1) {
136
137 double d = 0;
138 double c = 0;
139 double xgen = 0;
140 double v = 0;
141 double u = 0;
142 d = gamma -1./3.; c = 1./TMath::Sqrt(9.*d);
143
144 while(v <= 0.){
145 xgen = RooRandom::randomGenerator()->Gaus(); v = 1. + c*xgen;
146 }
148 if( u < 1.-.0331*(xgen*xgen)*(xgen*xgen) ) {
149 if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
150 x = ((d*v)* beta + mu) ;
151 break;
152 }
153 }
154 if( TMath::Log(u) < 0.5*xgen*xgen + d*(1.-v + TMath::Log(v)) ) {
155 if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
156 x = ((d*v)* beta + mu) ;
157 break;
158 }
159 }
160
161 }
162
163
164 return;
165}
166
167
#define d(i)
Definition: RSha256.hxx:102
#define c(i)
Definition: RSha256.hxx:101
int Int_t
Definition: RtypesCore.h:45
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:364
#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:61
virtual RooSpan< const double > getValues(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet=nullptr) const
by this change, please consult the release notes for ROOT 6.24 for guidance on how to make this trans...
Definition: RooAbsReal.cxx:311
Bool_t 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:35
virtual RooSpan< double > computeGamma(const RooAbsReal *, RunContext &, RooSpan< const double > x, RooSpan< const double > gamma, RooSpan< const double > beta, RooSpan< const double > mu)=0
Implementation of the Gamma PDF for RooFit/RooStats.
Definition: RooGamma.h:20
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooGamma.cxx:83
void generateEvent(Int_t code)
algorithm adapted from code example in: Marsaglia, G.
Definition: RooGamma.cxx:130
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooGamma.cxx:115
RooRealProxy beta
Definition: RooGamma.h:39
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooGamma.cxx:96
RooSpan< double > evaluateSpan(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet) const
Compute multiple values of Gamma PDF.
Definition: RooGamma.cxx:90
RooGamma()
Definition: RooGamma.h:22
RooRealProxy x
Definition: RooGamma.h:37
RooRealProxy gamma
Definition: RooGamma.h:38
RooRealProxy mu
Definition: RooGamma.h:40
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooGamma.cxx:104
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:53
A simple container to hold a batch of data values.
Definition: RooSpan.h:34
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 * dispatch
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:118
Double_t Log(Double_t x)
Definition: TMath.h:710
Double_t Sqrt(Double_t x)
Definition: TMath.h:641
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:2308
This struct enables passing computation data around between elements of a computation graph.
Definition: RunContext.h:31