beta.bc (1684B)
1 define sample_gamma(alpha){ 2 /* 3 A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001 4 https://dl.acm.org/doi/pdf/10.1145/358407.358414 5 see also the references/ folder 6 Note that the Wikipedia page for the gamma distribution includes a scaling parameter 7 k or beta 8 https://en.wikipedia.org/wiki/Gamma_distribution 9 such that gamma_k(alpha, k) = k * gamma(alpha) 10 or gamma_beta(alpha, beta) = gamma(alpha) / beta 11 So far I have not needed to use this, and thus the second parameter is by default 1. 12 */ 13 14 if (alpha >= 1) { 15 d = alpha - (1/3); 16 c = 1.0 / sqrt(9.0 * d); 17 while (1) { 18 v=-1 19 while(v<=0) { 20 x = sample_unit_normal(); 21 v = 1 + c * x; 22 } 23 v = v * v * v; 24 u = sample_unit_uniform(); 25 if (u < (1 - (0.0331 * (x * x * x * x)))) { /* Condition 1 */ 26 /* 27 the 0.0331 doesn't inspire much confidence 28 however, this isn't the whole story 29 by knowing that Condition 1 implies condition 2 30 we realize that this is just a way of making the algorithm faster 31 i.e., of not using the logarithms 32 */ 33 return d * v; 34 } 35 if (log(u, 2) < ((0.5 * (x * x)) + (d * (1 - v + log(v, 2))))) { /* Condition 2 */ 36 return d * v; 37 } 38 } 39 } else { 40 return sample_gamma(1 + alpha) * p(sample_unit_uniform(), 1 / alpha); 41 /* see note in p. 371 of https://dl.acm.org/doi/pdf/10.1145/358407.358414 */ 42 } 43 } 44 45 define sample_beta(a, b) 46 { 47 /* See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions */ 48 gamma_a = sample_gamma(a); 49 gamma_b = sample_gamma(b); 50 return gamma_a / (gamma_a + gamma_b); 51 }