ROOT   6.10/09 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
12 Implementation of the Gamma PDF for RooFit/RooStats.
13 \f[
14 f(x) = \frac{(x-\mu)^{\gamma-1} \cdot \exp^{(-(x-mu) / \beta}}{\Gamma(\gamma) \cdot \beta^{\gamma}}
15 \f]
16 defined for \f$x \geq 0 \f$ if \f$\mu = 0 \f$
17
18 Notes from Kyle Cranmer:
19
20 Wikipedia and several sources refer to the Gamma distribution as
21
22 \f[
23 G(\mu,\alpha,\beta) = \beta^\alpha \mu^{(\alpha-1)} \frac{e^{(-\beta \mu)}}{\Gamma(\alpha)}
24 \f]
25
26 Below is the correspondance:
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
36 Note, that for a model Pois(N|mu), a uniform prior on mu, and a measurement N
37 the posterior is in the Wikipedia parameterization Gamma(mu, alpha=N+1, beta=1)
38 thus for this implementation it is:
39
40 RooGamma(_x=mu,_gamma=N+1,_beta=1,_mu=0)
41
42 Also note, that in this case it is equivalent to
43 RooPoison(N,mu) and treating the function as a PDF in mu.
44 **/
45
46 #include "RooFit.h"
47
48 #include "Riostream.h"
49 #include "Riostream.h"
50 #include <math.h>
51
52 #include "RooGamma.h"
53 #include "RooAbsReal.h"
54 #include "RooRealVar.h"
55 #include "RooRandom.h"
56 #include "RooMath.h"
57
58 #include <iostream>
59 #include "TMath.h"
60
61 #include <Math/SpecFuncMathCore.h>
62 #include <Math/PdfFuncMathCore.h>
63 #include <Math/ProbFuncMathCore.h>
64
65 #include "TError.h"
66
67 using namespace std;
68
70
71 ////////////////////////////////////////////////////////////////////////////////
72
73 RooGamma::RooGamma(const char *name, const char *title,
74  RooAbsReal& _x, RooAbsReal& _gamma,
75  RooAbsReal& _beta, RooAbsReal& _mu) :
76  RooAbsPdf(name,title),
77  x("x","Observable",this,_x),
78  gamma("gamma","Mean",this,_gamma),
79  beta("beta","Width",this,_beta),
80  mu("mu","Para",this,_mu)
81 {
82 }
83
84 ////////////////////////////////////////////////////////////////////////////////
85
86 RooGamma::RooGamma(const RooGamma& other, const char* name) :
87  RooAbsPdf(other,name), x("x",this,other.x), gamma("mean",this,other.gamma),
88  beta("beta",this,other.beta), mu("mu",this,other.mu)
89 {
90 }
91
92 ////////////////////////////////////////////////////////////////////////////////
93
95 {
96  Double_t arg= x ;
97  Double_t ret = TMath::GammaDist(arg, gamma, mu, beta) ;
98  return ret ;
99 }
100
101 ////////////////////////////////////////////////////////////////////////////////
102
103 Int_t RooGamma::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
104 {
105  if (matchArgs(allVars,analVars,x)) return 1 ;
106  return 0 ;
107 }
108
109 ////////////////////////////////////////////////////////////////////////////////
110
111 Double_t RooGamma::analyticalIntegral(Int_t code, const char* rangeName) const
112 {
113  R__ASSERT(code==1) ;
114
115  //integral of the Gamma distribution via ROOT::Math
116  Double_t integral = ROOT::Math::gamma_cdf(x.max(rangeName), gamma, beta, mu) - ROOT::Math::gamma_cdf(x.min(rangeName), gamma, beta, mu);
117  return integral ;
118 }
119
120 ////////////////////////////////////////////////////////////////////////////////
121
122 Int_t RooGamma::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
123 {
124  if (matchArgs(directVars,generateVars,x)) return 1 ;
125  return 0 ;
126 }
127
128 ////////////////////////////////////////////////////////////////////////////////
129 /// algorithm adapted from code example in:
130 /// Marsaglia, G. and Tsang, W. W.
131 /// A Simple Method for Generating Gamma Variables
132 /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
133 ///
134 /// The speed of this algorithm depends on the speed of generating normal variates.
135 /// The algorithm is limited to \f$\gamma \geq 0 \f$ !
136
138 {
139  R__ASSERT(code==1) ;
140
141
142  while(1) {
143
144  double d = 0;
145  double c = 0;
146  double xgen = 0;
147  double v = 0;
148  double u = 0;
149  d = gamma -1./3.; c = 1./TMath::Sqrt(9.*d);
150
151  while(v <= 0.){
152  xgen = RooRandom::randomGenerator()->Gaus(); v = 1. + c*xgen;
153  }
154  v = v*v*v; u = RooRandom::randomGenerator()->Uniform();
155  if( u < 1.-.0331*(xgen*xgen)*(xgen*xgen) ) {
156  if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
157  x = ((d*v)* beta + mu) ;
158  break;
159  }
160  }
161  if( TMath::Log(u) < 0.5*xgen*xgen + d*(1.-v + TMath::Log(v)) ) {
162  if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
163  x = ((d*v)* beta + mu) ;
164  break;
165  }
166  }
167
168  }
169
170
171  return;
172 }
RooRealProxy x
Definition: RooGamma.h:39
Double_t Log(Double_t x)
Definition: TMath.h:649
Double_t evaluate() const
Definition: RooGamma.cxx:94
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:235
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooRealProxy beta
Definition: RooGamma.h:41
#define R__ASSERT(e)
Definition: TError.h:96
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
RooRealProxy gamma
Definition: RooGamma.h:40
double beta(double x, double y)
Calculates the beta function.
Double_t x[n]
Definition: legend1.C:17
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
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:103
double gamma(double x)
void generateEvent(Int_t code)
algorithm adapted from code example in: Marsaglia, G.
Definition: RooGamma.cxx:137
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
Definition: RooGamma.cxx:111
SVector< double, 2 > v
Definition: Dict.h:5
RooRealProxy mu
Definition: RooGamma.h:42
Implementation of the Gamma PDF for RooFit/RooStats.
Definition: RooGamma.h:22
Double_t Mean(Long64_t n, const T *a, const Double_t *w=0)
Definition: TMath.h:973
RooGamma()
Definition: RooGamma.h:24
#define ClassImp(name)
Definition: Rtypes.h:336
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:122
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
double gamma_cdf(double x, double alpha, double theta, double x0=0)
Cumulative distribution function of the gamma distribution (lower tail).
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
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:2296
Double_t Sqrt(Double_t x)
Definition: TMath.h:591