52 lines
1.6 KiB
Plaintext
52 lines
1.6 KiB
Plaintext
define sample_gamma(alpha){
|
|
/*
|
|
A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001
|
|
https://dl.acm.org/doi/pdf/10.1145/358407.358414
|
|
see also the references/ folder
|
|
Note that the Wikipedia page for the gamma distribution includes a scaling parameter
|
|
k or beta
|
|
https://en.wikipedia.org/wiki/Gamma_distribution
|
|
such that gamma_k(alpha, k) = k * gamma(alpha)
|
|
or gamma_beta(alpha, beta) = gamma(alpha) / beta
|
|
So far I have not needed to use this, and thus the second parameter is by default 1.
|
|
*/
|
|
|
|
if (alpha >= 1) {
|
|
d = alpha - (1/3);
|
|
c = 1.0 / sqrt(9.0 * d);
|
|
while (1) {
|
|
v=-1
|
|
while(v<=0) {
|
|
x = sample_unit_normal();
|
|
v = 1 + c * x;
|
|
}
|
|
v = v * v * v;
|
|
u = sample_unit_uniform();
|
|
if (u < (1 - (0.0331 * (x * x * x * x)))) { /* Condition 1 */
|
|
/*
|
|
the 0.0331 doesn't inspire much confidence
|
|
however, this isn't the whole story
|
|
by knowing that Condition 1 implies condition 2
|
|
we realize that this is just a way of making the algorithm faster
|
|
i.e., of not using the logarithms
|
|
*/
|
|
return d * v;
|
|
}
|
|
if (log(u, 2) < ((0.5 * (x * x)) + (d * (1 - v + log(v, 2))))) { /* Condition 2 */
|
|
return d * v;
|
|
}
|
|
}
|
|
} else {
|
|
return sample_gamma(1 + alpha) * p(sample_unit_uniform(), 1 / alpha);
|
|
/* see note in p. 371 of https://dl.acm.org/doi/pdf/10.1145/358407.358414 */
|
|
}
|
|
}
|
|
|
|
define sample_beta(a, b)
|
|
{
|
|
/* See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions */
|
|
gamma_a = sample_gamma(a);
|
|
gamma_b = sample_gamma(b);
|
|
return gamma_a / (gamma_a + gamma_b);
|
|
}
|