8 * Created by Pat Schloss on 7/6/11.
9 * Copyright 2011 Patrick D. Schloss. All rights reserved.
13 #include "randomnumber.h"
16 /**************************************************************************************************/
18 RandomNumberGenerator::RandomNumberGenerator(){
19 srand( (unsigned)time( NULL ) );
22 /**************************************************************************************************/
24 float RandomNumberGenerator::randomUniform(){
26 float randUnif = 0.0000;
28 while(randUnif == 0.0000){
30 randUnif = rand() / (float)RAND_MAX;
37 /**************************************************************************************************/
39 //Code shamelessly swiped and modified from Numerical Recipes in C++
41 /**************************************************************************************************/
43 float RandomNumberGenerator::randomExp(){
45 float randExp = 0.0000;
47 while(randExp == 0.0000){
49 randExp = -log(randomUniform());
56 /**************************************************************************************************/
58 //Code shamelessly swiped and modified from Numerical Recipes in C++
60 /**************************************************************************************************/
62 float RandomNumberGenerator::randomNorm(){
67 x = 2.0 * randomUniform() - 1.0;
68 y = 2.0 * randomUniform() - 1.0;
70 rsquare = x * x + y * y;
71 } while(rsquare >= 1.0 || rsquare == 0.0);
73 float fac = sqrt(-2.0 * log(rsquare)/rsquare);
79 /**************************************************************************************************/
81 * Slightly modified version of:
83 * Mathlib : A C Library of Special Functions
84 * Copyright (C) 1998 Ross Ihaka
85 * Copyright (C) 2000-2005 The R Development Core Team
87 * This program is free software; you can redistribute it and/or modify
88 * it under the terms of the GNU General Public License as published by
89 * the Free Software Foundation; either version 2 of the License, or
90 * (at your option) any later version.
92 * This program is distributed in the hope that it will be useful,
93 * but WITHOUT ANY WARRANTY; without even the implied warranty of
94 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
95 * GNU General Public License for more details.
97 * You should have received a copy of the GNU General Public License
98 * along with this program; if not, a copy is available at
99 * http://www.r-project.org/Licenses/
104 * float rgamma(float a, float scale);
108 * Random variates from the gamma distribution.
112 * [1] Shape parameter a >= 1. Algorithm GD in:
114 * Ahrens, J.H. and Dieter, U. (1982).
115 * Generating gamma variates by a modified
116 * rejection technique.
117 * Comm. ACM, 25, 47-54.
120 * [2] Shape parameter 0 < a < 1. Algorithm GS in:
122 * Ahrens, J.H. and Dieter, U. (1974).
123 * Computer methods for sampling from gamma, beta,
124 * poisson and binomial distributions.
125 * Computing, 12, 223-246.
127 * Input: a = parameter (mean) of the standard gamma distribution.
128 * Output: a variate from the gamma(a)-distribution
131 float RandomNumberGenerator::randomGamma(float a)
134 const static float sqrt32 = 5.656854;
135 const static float exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
138 /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k))
139 * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k)
140 * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k)
142 const static float q1 = 0.04166669;
143 const static float q2 = 0.02083148;
144 const static float q3 = 0.00801191;
145 const static float q4 = 0.00144121;
146 const static float q5 = -7.388e-5;
147 const static float q6 = 2.4511e-4;
148 const static float q7 = 2.424e-4;
150 const static float a1 = 0.3333333;
151 const static float a2 = -0.250003;
152 const static float a3 = 0.2000062;
153 const static float a4 = -0.1662921;
154 const static float a5 = 0.1423657;
155 const static float a6 = -0.1367177;
156 const static float a7 = 0.1233795;
158 /* State variables [FIXME for threading!] :*/
159 static float aa = 0.;
160 static float aaa = 0.;
161 static float s, s2, d; /* no. 1 (step 1) */
162 static float q0, b, si, c;/* no. 2 (step 4) */
164 float e, p, q, r, t, u, v, w, x, ret_val;
166 if (a <= 0.0 || scale <= 0.0){ cout << "error alpha or scale parameter are less than zero." << endl; exit(1); }
168 if (a < 1.) { /* GS algorithm for parameters a < 1 */
169 e = 1.0 + exp_m1 * a;
171 p = e * randomUniform();
173 x = -log((e - p) / a);
174 if (randomExp() >= (1.0 - a) * log(x))
178 if (randomExp() >= x)
185 /* --- a >= 1 : GD algorithm --- */
187 /* Step 1: Recalculations of s2, s, d if a has changed */
192 d = sqrt32 - s * 12.0;
194 /* Step 2: t = standard normal deviate,
195 x = (s,1/2) -normal deviate. */
197 /* immediate acceptance (i) */
202 return scale * ret_val;
204 /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
206 if (d * u <= t * t * t)
207 return scale * ret_val;
209 /* Step 4: recalculations of q0, b, si, c if necessary */
214 q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
217 /* Approximation depending on size of parameter a */
218 /* The constants in the expressions for b, si and c */
219 /* were established by numerical experiments */
222 b = 0.463 + s + 0.178 * s2;
224 c = 0.195 / s - 0.079 + 0.16 * s;
225 } else if (a <= 13.022) {
226 b = 1.654 + 0.0076 * s2;
227 si = 1.68 / s + 0.275;
228 c = 0.062 / s + 0.024;
235 /* Step 5: no quotient test if x not positive */
238 /* Step 6: calculation of v and quotient q */
241 q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
242 + a3) * v + a2) * v + a1) * v;
244 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
247 /* Step 7: quotient acceptance (q) */
248 if (log(1.0 - u) <= q)
249 return scale * ret_val;
253 /* Step 8: e = standard exponential deviate
254 * u = 0,1 -uniform deviate
255 * t = (b,si)-float exponential (laplace) sample */
263 /* Step 9: rejection if t < tau(1) = -0.71874483771719 */
264 if (t >= -0.71874483771719) {
265 /* Step 10: calculation of v and quotient q */
268 q = q0 + 0.5 * t * t *
269 ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
272 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
273 /* Step 11: hat acceptance (h) */
274 /* (if q not positive go to step 8) */
277 /* ^^^^^ original code had approximation with rel.err < 2e-7 */
278 /* if t is rejected sample again at step 8 */
279 if (c * fabs(u) <= w * exp(e - 0.5 * t * t))
283 } /* for(;;) .. until `t' is accepted */
285 return scale * x * x;
288 /**************************************************************************************************/
290 // essentially swiped from http://en.wikipedia.org/wiki/Dirichlet_distribution#Random_number_generation
292 /**************************************************************************************************/
294 vector<float> RandomNumberGenerator::randomDirichlet(vector<float> alphas){
296 int nAlphas = (int)alphas.size();
297 vector<float> dirs(nAlphas, 0.0000);
300 for(int i=0;i<nAlphas;i++){
301 dirs[i] = randomGamma(alphas[i]);
305 for(int i=0;i<nAlphas;i++){
313 /**************************************************************************************************/