ROOT   Reference Guide
Searching...
No Matches
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
91{
93}
94
95////////////////////////////////////////////////////////////////////////////////
96/// Compute multiple values of Gamma PDF.
97void RooGamma::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
98{
100 dispatch->compute(stream, RooBatchCompute::Gamma, output, nEvents,
101 {dataMap.at(x), dataMap.at(gamma), dataMap.at(beta), dataMap.at(mu)});
102}
103
104////////////////////////////////////////////////////////////////////////////////
105
106Int_t RooGamma::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
107{
108 if (matchArgs(allVars,analVars,x)) return 1 ;
109 return 0 ;
110}
111
112////////////////////////////////////////////////////////////////////////////////
113
114double RooGamma::analyticalIntegral(Int_t /*code*/, const char *rangeName) const
115{
116 // integral of the Gamma distribution via ROOT::Math
117 return ROOT::Math::gamma_cdf(x.max(rangeName), gamma, beta, mu) -
118 ROOT::Math::gamma_cdf(x.min(rangeName), gamma, beta, mu);
119}
120
121////////////////////////////////////////////////////////////////////////////////
122
123std::string RooGamma::buildCallToAnalyticIntegral(Int_t /*code*/, const char *rangeName,
125{
126 const std::string a = ctx.buildCall("ROOT::Math::gamma_cdf", x.max(rangeName), gamma, beta, mu);
127 const std::string b = ctx.buildCall("ROOT::Math::gamma_cdf", x.min(rangeName), gamma, beta, mu);
128 return a + " - " + b;
129}
130
131namespace {
132
133inline double randomGamma(double gamma, double beta, double mu, double xmin, double xmax)
134{
135 while (true) {
136
137 double d = gamma - 1. / 3.;
138 double c = 1. / TMath::Sqrt(9. * d);
139 double xgen = 0;
140 double v = 0;
141
142 while (v <= 0.) {
144 v = 1. + c * xgen;
145 }
146 v = v * v * v;
147 double u = RooRandom::randomGenerator()->Uniform();
148 if (u < 1. - .0331 * (xgen * xgen) * (xgen * xgen)) {
149 double x = ((d * v) * beta + mu);
150 if ((x < xmax) && (x > xmin)) {
151 return x;
152 }
153 }
154 if (TMath::Log(u) < 0.5 * xgen * xgen + d * (1. - v + TMath::Log(v))) {
155 double x = ((d * v) * beta + mu);
156 if ((x < xmax) && (x > xmin)) {
157 return x;
158 }
159 }
160 }
161}
162
163} // namespace
164
165////////////////////////////////////////////////////////////////////////////////
166
167Int_t RooGamma::getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
168{
169 if (matchArgs(directVars, generateVars, x))
170 return 1;
171 return 0;
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// algorithm adapted from code example in:
176/// Marsaglia, G. and Tsang, W. W.
177/// A Simple Method for Generating Gamma Variables
178/// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
179///
180/// The speed of this algorithm depends on the speed of generating normal variates.
181/// The algorithm is limited to \f$\gamma \geq 0 \f$ !
182
184{
185 if (gamma >= 1) {
186 x = randomGamma(gamma, beta, mu, x.min(), x.max());
187 return;
188 }
189
190 double xVal = 0.0;
191 bool accepted = false;
192
193 while (!accepted) {
194 double u = RooRandom::randomGenerator()->Uniform();
195 double tmp = randomGamma(1 + gamma, beta, mu, 0, std::numeric_limits<double>::infinity());
196 xVal = tmp * std::pow(u, 1.0 / gamma);
197 accepted = xVal < x.max() && xVal > x.min();
198 }
199
200 x = xVal;
201}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define ClassImp(name)
Definition Rtypes.h:377
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:61
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:55
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, ArgVector &)=0
A class to maintain the context for squashing of RooFit models into code.
std::string buildCall(std::string const &funcname, Args_t const &...args)
Build the code to call the function with name funcname, passing some arguments.
void addResult(RooAbsArg const *key, std::string const &value)
A function to save an expression that includes/depends on the result of the input node.
RooSpan< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
Definition DataMap.cxx:21
Implementation of the Gamma PDF for RooFit/RooStats.
Definition RooGamma.h:20
std::string buildCallToAnalyticIntegral(Int_t code, const char *rangeName, RooFit::Detail::CodeSquashContext &ctx) const override
This function defines the analytical integral translation for the class.
Definition RooGamma.cxx:123
RooRealProxy beta
Definition RooGamma.h:43
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:167
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of Gamma PDF.
Definition RooGamma.cxx:97
RooGamma()
Definition RooGamma.h:22
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition RooGamma.cxx:114
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:41
void translate(RooFit::Detail::CodeSquashContext &ctx) const override
This function defines a translation for each RooAbsReal based object that can be used to express the ...
Definition RooGamma.cxx:90
void generateEvent(Int_t code) override
algorithm adapted from code example in: Marsaglia, G.
Definition RooGamma.cxx:183
RooRealProxy gamma
Definition RooGamma.h:42
RooRealProxy mu
Definition RooGamma.h:44
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition RooGamma.cxx:106
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:51
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower 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_t x[n]
Definition legend1.C:17
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.
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition TMath.h:756
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
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:2342
static void output()