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 "RooAbsReal.h"
49#include "RooRealVar.h"
50#include "RooRandom.h"
51#include "RooMath.h"
52#include "RooHelpers.h"
53#include "BatchHelpers.h"
55
56#include "TMath.h"
60
61#include <iostream>
62#include <cmath>
63using namespace std;
64
66
67////////////////////////////////////////////////////////////////////////////////
68
69RooGamma::RooGamma(const char *name, const char *title,
70 RooAbsReal& _x, RooAbsReal& _gamma,
71 RooAbsReal& _beta, RooAbsReal& _mu) :
72 RooAbsPdf(name,title),
73 x("x","Observable",this,_x),
74 gamma("gamma","Mean",this,_gamma),
75 beta("beta","Width",this,_beta),
76 mu("mu","Para",this,_mu)
77{
78 RooHelpers::checkRangeOfParameters(this, {&_gamma, &_beta}, 0.);
79}
80
81////////////////////////////////////////////////////////////////////////////////
82
83RooGamma::RooGamma(const RooGamma& other, const char* name) :
84 RooAbsPdf(other,name), x("x",this,other.x), gamma("mean",this,other.gamma),
85 beta("beta",this,other.beta), mu("mu",this,other.mu)
86{
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
92{
93 return TMath::GammaDist(x, gamma, mu, beta) ;
94}
95
96////////////////////////////////////////////////////////////////////////////////
97
98namespace {
99//Author: Emmanouil Michalainas, CERN 22 August 2019
100
101template<class Tx, class Tgamma, class Tbeta, class Tmu>
102void compute( size_t batchSize,
103 double * __restrict output,
104 Tx X, Tgamma G, Tbeta B, Tmu M)
105{
106 constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
107 for (size_t i=0; i<batchSize; i++) {
108 if (X[i]<M[i] || G[i] <= 0.0 || B[i] <= 0.0) {
109 output[i] = NaN;
110 }
111 if (X[i] == M[i]) {
112 output[i] = (G[i]==1.0)/B[i];
113 }
114 else {
115 output[i] = 0.0;
116 }
117 }
118
119 if (G.isBatch()) {
120 for (size_t i=0; i<batchSize; i++) {
121 if (output[i] == 0.0) {
122 output[i] = -std::lgamma(G[i]);
123 }
124 }
125 }
126 else {
127 double gamma = -std::lgamma(G[2019]);
128 for (size_t i=0; i<batchSize; i++) {
129 if (output[i] == 0.0) {
130 output[i] = gamma;
131 }
132 }
133 }
134
135 for (size_t i=0; i<batchSize; i++) {
136 if (X[i] != M[i]) {
137 const double invBeta = 1/B[i];
138 double arg = (X[i]-M[i])*invBeta;
139 output[i] -= arg;
140 arg = _rf_fast_log(arg);
141 output[i] += arg*(G[i]-1);
142 output[i] = _rf_fast_exp(output[i]);
143 output[i] *= invBeta;
144 }
145 }
146}
147};
148
149RooSpan<double> RooGamma::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
150 using namespace BatchHelpers;
151
152 EvaluateInfo info = getInfo( {&x, &gamma, &beta, &mu}, begin, batchSize );
153 if (info.nBatches == 0) {
154 return {};
155 }
156 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
157 auto xData = x.getValBatch(begin, info.size);
158
159 if (info.nBatches==1 && !xData.empty()) {
160 compute(info.size, output.data(), xData.data(),
164 }
165 else {
166 compute(info.size, output.data(),
171 }
172 return output;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176
177Int_t RooGamma::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
178{
179 if (matchArgs(allVars,analVars,x)) return 1 ;
180 return 0 ;
181}
182
183////////////////////////////////////////////////////////////////////////////////
184
185Double_t RooGamma::analyticalIntegral(Int_t code, const char* rangeName) const
186{
187 R__ASSERT(code==1) ;
188
189 //integral of the Gamma distribution via ROOT::Math
190 Double_t integral = ROOT::Math::gamma_cdf(x.max(rangeName), gamma, beta, mu) - ROOT::Math::gamma_cdf(x.min(rangeName), gamma, beta, mu);
191 return integral ;
192}
193
194////////////////////////////////////////////////////////////////////////////////
195
196Int_t RooGamma::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
197{
198 if (matchArgs(directVars,generateVars,x)) return 1 ;
199 return 0 ;
200}
201
202////////////////////////////////////////////////////////////////////////////////
203/// algorithm adapted from code example in:
204/// Marsaglia, G. and Tsang, W. W.
205/// A Simple Method for Generating Gamma Variables
206/// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
207///
208/// The speed of this algorithm depends on the speed of generating normal variates.
209/// The algorithm is limited to \f$\gamma \geq 0 \f$ !
210
212{
213 R__ASSERT(code==1) ;
214
215
216 while(1) {
217
218 double d = 0;
219 double c = 0;
220 double xgen = 0;
221 double v = 0;
222 double u = 0;
223 d = gamma -1./3.; c = 1./TMath::Sqrt(9.*d);
224
225 while(v <= 0.){
226 xgen = RooRandom::randomGenerator()->Gaus(); v = 1. + c*xgen;
227 }
229 if( u < 1.-.0331*(xgen*xgen)*(xgen*xgen) ) {
230 if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
231 x = ((d*v)* beta + mu) ;
232 break;
233 }
234 }
235 if( TMath::Log(u) < 0.5*xgen*xgen + d*(1.-v + TMath::Log(v)) ) {
236 if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
237 x = ((d*v)* beta + mu) ;
238 break;
239 }
240 }
241
242 }
243
244
245 return;
246}
247
248
#define NaN
#define d(i)
Definition: RSha256.hxx:102
#define c(i)
Definition: RSha256.hxx:101
double _rf_fast_exp(double x)
double _rf_fast_log(double x)
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
#define R__ASSERT(e)
Definition: TError.h:96
char name[80]
Definition: TGX11.cxx:109
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
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:59
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Implementation of the Gamma PDF for RooFit/RooStats.
Definition: RooGamma.h:22
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooGamma.cxx:91
void generateEvent(Int_t code)
algorithm adapted from code example in: Marsaglia, G.
Definition: RooGamma.cxx:211
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:196
RooRealProxy beta
Definition: RooGamma.h:41
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:177
RooGamma()
Definition: RooGamma.h:24
RooRealProxy x
Definition: RooGamma.h:39
RooRealProxy gamma
Definition: RooGamma.h:40
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
Definition: RooGamma.cxx:149
RooRealProxy mu
Definition: RooGamma.h:42
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooGamma.cxx:185
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
Double_t min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize) const
Double_t 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:263
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:635
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 lgamma(double x)
Calculates the logarithm of the gamma function.
Double_t x[n]
Definition: legend1.C:17
#define G(x, y, z)
EvaluateInfo getInfo(std::vector< const RooRealProxy * > parameters, size_t begin, size_t batchSize)
static double B[]
double gamma(double x)
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 extraMessage="")
Check if the parameters have a range, and warn if the range extends below / above the set limits.
Definition: RooHelpers.cxx:75
Double_t Log(Double_t x)
Definition: TMath.h:750
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
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:2315
static void output(int code)
Definition: gifencode.c:226