Logo ROOT  
Reference Guide
 
Loading...
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
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("gamma",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{
92 ctx.addResult(this, ctx.buildCall("TMath::GammaDist", x, gamma, mu, beta));
93}
94
95////////////////////////////////////////////////////////////////////////////////
96/// Compute multiple values of Gamma PDF.
98{
100 {ctx.at(x), ctx.at(gamma), ctx.at(beta), ctx.at(mu)});
101}
102
103////////////////////////////////////////////////////////////////////////////////
104
105Int_t RooGamma::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
106{
107 if (matchArgs(allVars,analVars,x)) return 1 ;
108 return 0 ;
109}
110
111////////////////////////////////////////////////////////////////////////////////
112
113double RooGamma::analyticalIntegral(Int_t /*code*/, const char *rangeName) const
114{
115 // integral of the Gamma distribution via ROOT::Math
116 return ROOT::Math::gamma_cdf(x.max(rangeName), gamma, beta, mu) -
117 ROOT::Math::gamma_cdf(x.min(rangeName), gamma, beta, mu);
118}
119
120////////////////////////////////////////////////////////////////////////////////
121
122std::string RooGamma::buildCallToAnalyticIntegral(Int_t /*code*/, const char *rangeName,
124{
125 const std::string a = ctx.buildCall("ROOT::Math::gamma_cdf", x.max(rangeName), gamma, beta, mu);
126 const std::string b = ctx.buildCall("ROOT::Math::gamma_cdf", x.min(rangeName), gamma, beta, mu);
127 return a + " - " + b;
128}
129
130namespace {
131
132inline double randomGamma(double gamma, double beta, double mu, double xmin, double xmax)
133{
134 while (true) {
135
136 double d = gamma - 1. / 3.;
137 double c = 1. / std::sqrt(9. * d);
138 double xgen = 0;
139 double v = 0;
140
141 while (v <= 0.) {
143 v = 1. + c * xgen;
144 }
145 v = v * v * v;
146 double u = RooRandom::randomGenerator()->Uniform();
147 if (u < 1. - .0331 * (xgen * xgen) * (xgen * xgen)) {
148 double x = ((d * v) * beta + mu);
149 if ((x < xmax) && (x > xmin)) {
150 return x;
151 }
152 }
153 if (std::log(u) < 0.5 * xgen * xgen + d * (1. - v + TMath::Log(v))) {
154 double x = ((d * v) * beta + mu);
155 if ((x < xmax) && (x > xmin)) {
156 return x;
157 }
158 }
159 }
160}
161
162} // namespace
163
164////////////////////////////////////////////////////////////////////////////////
165
166Int_t RooGamma::getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
167{
168 if (matchArgs(directVars, generateVars, x))
169 return 1;
170 return 0;
171}
172
173////////////////////////////////////////////////////////////////////////////////
174/// algorithm adapted from code example in:
175/// Marsaglia, G. and Tsang, W. W.
176/// A Simple Method for Generating Gamma Variables
177/// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
178///
179/// The speed of this algorithm depends on the speed of generating normal variates.
180/// The algorithm is limited to \f$ \gamma \geq 0 \f$ !
181
183{
184 if (gamma >= 1) {
185 x = randomGamma(gamma, beta, mu, x.min(), x.max());
186 return;
187 }
188
189 double xVal = 0.0;
190 bool accepted = false;
191
192 while (!accepted) {
193 double u = RooRandom::randomGenerator()->Uniform();
194 double tmp = randomGamma(1 + gamma, beta, mu, 0, std::numeric_limits<double>::infinity());
195 xVal = tmp * std::pow(u, 1.0 / gamma);
196 accepted = xVal < x.max() && xVal > x.min();
197 }
198
199 x = xVal;
200}
#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
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
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
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.
std::span< double > output()
RooBatchCompute::Config config(RooAbsArg const *arg) const
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:122
RooRealProxy beta
Definition RooGamma.h:42
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:166
void doEval(RooFit::EvalContext &) 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:113
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:40
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:182
RooRealProxy gamma
Definition RooGamma.h:41
RooRealProxy mu
Definition RooGamma.h:43
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:105
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:48
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:275
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition TRandom.cxx:682
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
void compute(Config cfg, Computer comp, std::span< double > output, VarSpan vars, ArgSpan extraArgs={})
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 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:2347