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