|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
//FAMOS headers
#include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
GammaFunctionGenerator::GammaFunctionGenerator() {
xmax = 30.;
for (unsigned i = 1; i <= 12; ++i) {
// For all let's put the limit at 2.*(alpha-1) (alpha-1 is the max of the dist)
approxLimit.push_back(2 * ((double)i));
myIncompleteGamma.a().setValue((double)i);
integralToApproxLimit.push_back(myIncompleteGamma(approxLimit[i - 1]));
theGammas.push_back(GammaNumericalGenerator((double)i, 1., 0, approxLimit[i - 1] + 1.));
}
coreCoeff.push_back(0.); // alpha=1 not used
coreCoeff.push_back(1. / 8.24659e-01);
coreCoeff.push_back(1. / 7.55976e-01);
coreCoeff.push_back(1. / 7.12570e-01);
coreCoeff.push_back(1. / 6.79062e-01);
coreCoeff.push_back(1. / 6.65496e-01);
coreCoeff.push_back(1. / 6.48736e-01);
coreCoeff.push_back(1. / 6.25185e-01);
coreCoeff.push_back(1. / 6.09188e-01);
coreCoeff.push_back(1. / 6.06221e-01);
coreCoeff.push_back(1. / 6.05057e-01);
}
GammaFunctionGenerator::~GammaFunctionGenerator() {}
double GammaFunctionGenerator::shoot(RandomEngineAndDistribution const* random) const {
if (alpha < 0.)
return -1.;
if (badRange)
return xmin / beta;
if (alpha < 12) {
if (alpha == na) {
return gammaInt(random) / beta;
} else if (na == 0) {
return gammaFrac(random) / beta;
} else {
double gi = gammaInt(random);
double gf = gammaFrac(random);
return (gi + gf) / beta;
}
} else {
// an other generator has to be used in such a case
return -1.;
}
}
double GammaFunctionGenerator::gammaFrac(RandomEngineAndDistribution const* random) const {
/* This is exercise 16 from Knuth; see page 135, and the solution is
on page 551. */
double p, q, x, u, v;
p = M_E / (frac + M_E);
do {
u = random->flatShoot();
v = random->flatShoot();
if (u < p) {
x = exp((1 / frac) * log(v));
q = exp(-x);
} else {
x = 1 - log(v);
q = exp((frac - 1) * log(x));
}
} while (random->flatShoot() >= q);
return x;
}
double GammaFunctionGenerator::gammaInt(RandomEngineAndDistribution const* random) const {
// Exponential distribution : no approximation
if (na == 1) {
return xmin - log(random->flatShoot());
}
unsigned gn = na - 1;
// are we sure to be in the tail
if (coreProba == 0.)
return xmin - coreCoeff[gn] * log(random->flatShoot());
// core-tail interval
if (random->flatShoot() < coreProba) {
return theGammas[gn].gamma_lin(random);
}
// std::cout << " Tail called " << std::endl;
return approxLimit[gn] - coreCoeff[gn] * log(random->flatShoot());
}
void GammaFunctionGenerator::setParameters(double a, double b, double xm) {
// std::cout << "Setting parameters " << std::endl;
alpha = a;
beta = b;
xmin = xm * beta;
if (xm > xmax) {
badRange = true;
return;
}
badRange = false;
na = 0;
if (alpha > 0. && alpha < 12)
na = (unsigned)floor(alpha);
frac = alpha - na;
// Now calculate the probability to shoot between approxLimit and xmax
// The Incomplete gamma is normalized to 1
if (na <= 1)
return;
myIncompleteGamma.a().setValue((double)na);
unsigned gn = na - 1;
// std::cout << " na " << na << " xm " << xm << " beta " << beta << " xmin " << xmin << " approxLimit " << approxLimit[gn] << std::endl;
if (xmin > approxLimit[gn]) {
coreProba = 0.;
} else {
double tmp = (xmin != 0.) ? myIncompleteGamma(xmin) : 0.;
coreProba = (integralToApproxLimit[gn] - tmp) / (1. - tmp);
theGammas[gn].setSubInterval(xmin, approxLimit[gn]);
}
// std::cout << " Proba " << coreProba << std::endl;
}
|