randomnumber.cpp

   1 #ifndef RANDOMNUMBER_H
   2 #define RANDOMNUMBER_H
   3
   4 /*
   5  *  randomnumber.cpp
   6  *
   7  *
   8  *  Created by Pat Schloss on 7/6/11.
   9  *  Copyright 2011 Patrick D. Schloss. All rights reserved.
  10  *
  11  */
  12
  13 #include "randomnumber.h"
  14 #include <cmath>
  15
  16 /**************************************************************************************************/
  17
  18 RandomNumberGenerator::RandomNumberGenerator(){
  19     srand( (unsigned)time( NULL ) );
  20 }
  21
  22 /**************************************************************************************************/
  23
  24 float RandomNumberGenerator::randomUniform(){
  25
  26         float randUnif = 0.0000;
  27
  28         while(randUnif == 0.0000){
  29
  30                 randUnif = rand() / (float)RAND_MAX;
  31
  32         }
  33
  34         return randUnif;
  35 }
  36
  37 /**************************************************************************************************/
  38 //
  39 //Code shamelessly swiped and modified from Numerical Recipes in C++
  40 //
  41 /**************************************************************************************************/
  42
  43 float RandomNumberGenerator::randomExp(){
  44
  45         float randExp = 0.0000;
  46
  47         while(randExp == 0.0000){
  48
  49                 randExp = -log(randomUniform());
  50
  51         }
  52
  53         return randExp;
  54 }
  55
  56 /**************************************************************************************************/
  57 //
  58 //Code shamelessly swiped and modified from Numerical Recipes in C++
  59 //
  60 /**************************************************************************************************/
  61
  62 float RandomNumberGenerator::randomNorm(){
  63
  64         float x, y, rsquare;
  65
  66         do{
  67                 x = 2.0 * randomUniform() - 1.0;
  68                 y = 2.0 * randomUniform() - 1.0;
  69
  70                 rsquare = x * x + y * y;
  71         } while(rsquare >= 1.0 || rsquare == 0.0);
  72
  73         float fac = sqrt(-2.0 * log(rsquare)/rsquare);
  74
  75         return x * fac;
  76 }
  77
  78
  79 /**************************************************************************************************/
  80 /*
  81  *      Slightly modified version of:
  82  *
  83  *  Mathlib : A C Library of Special Functions
  84  *  Copyright (C) 1998 Ross Ihaka
  85  *  Copyright (C) 2000-2005 The R Development Core Team
  86  *
  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.
  91  *
  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.
  96  *
  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/
 100  *
 101  *  SYNOPSIS
 102  *
 103  *    #include <Rmath.h>
 104  *    float rgamma(float a, float scale);
 105  *
 106  *  DESCRIPTION
 107  *
 108  *    Random variates from the gamma distribution.
 109  *
 110  *  REFERENCES
 111  *
 112  *    [1] Shape parameter a >= 1.  Algorithm GD in:
 113  *
 114  *        Ahrens, J.H. and Dieter, U. (1982).
 115  *        Generating gamma variates by a modified
 116  *        rejection technique.
 117  *        Comm. ACM, 25, 47-54.
 118  *
 119  *
 120  *    [2] Shape parameter 0 < a < 1. Algorithm GS in:
 121  *
 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.
 126  *
 127  *    Input: a = parameter (mean) of the standard gamma distribution.
 128  *    Output: a variate from the gamma(a)-distribution
 129  */
 130
 131 float RandomNumberGenerator::randomGamma(float a)
 132 {
 133         /* Constants : */
 134     const static float sqrt32 = 5.656854;
 135     const static float exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
 136         float scale = 1.0;
 137
 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)
 141      */
 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;
 149
 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;
 157
 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) */
 163
 164     float e, p, q, r, t, u, v, w, x, ret_val;
 165
 166     if (a <= 0.0 || scale <= 0.0){      cout << "error alpha or scale parameter are less than zero." << endl; exit(1);  }
 167
 168     if (a < 1.) { /* GS algorithm for parameters a < 1 */
 169         e = 1.0 + exp_m1 * a;
 170         for(;;) {
 171             p = e * randomUniform();
 172             if (p >= 1.0) {
 173                 x = -log((e - p) / a);
 174                 if (randomExp() >= (1.0 - a) * log(x))
 175                     break;
 176             } else {
 177                 x = exp(log(p) / a);
 178                 if (randomExp() >= x)
 179                     break;
 180             }
 181         }
 182         return scale * x;
 183     }
 184
 185     /* --- a >= 1 : GD algorithm --- */
 186
 187     /* Step 1: Recalculations of s2, s, d if a has changed */
 188     if (a != aa) {
 189         aa = a;
 190         s2 = a - 0.5;
 191         s = sqrt(s2);
 192         d = sqrt32 - s * 12.0;
 193     }
 194     /* Step 2: t = standard normal deviate,
 195          x = (s,1/2) -normal deviate. */
 196
 197     /* immediate acceptance (i) */
 198     t = randomNorm();
 199     x = s + 0.5 * t;
 200     ret_val = x * x;
 201     if (t >= 0.0)
 202         return scale * ret_val;
 203
 204     /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
 205     u = randomUniform();
 206     if (d * u <= t * t * t)
 207         return scale * ret_val;
 208
 209     /* Step 4: recalculations of q0, b, si, c if necessary */
 210
 211     if (a != aaa) {
 212         aaa = a;
 213         r = 1.0 / a;
 214         q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
 215                + q2) * r + q1) * r;
 216
 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 */
 220
 221         if (a <= 3.686) {
 222             b = 0.463 + s + 0.178 * s2;
 223             si = 1.235;
 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;
 229         } else {
 230             b = 1.77;
 231             si = 0.75;
 232             c = 0.1515 / s;
 233         }
 234     }
 235     /* Step 5: no quotient test if x not positive */
 236
 237     if (x > 0.0) {
 238         /* Step 6: calculation of v and quotient q */
 239         v = t / (s + s);
 240         if (fabs(v) <= 0.25)
 241             q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
 242                                       + a3) * v + a2) * v + a1) * v;
 243         else
 244             q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
 245
 246
 247         /* Step 7: quotient acceptance (q) */
 248         if (log(1.0 - u) <= q)
 249             return scale * ret_val;
 250     }
 251
 252     for(;;) {
 253         /* Step 8: e = standard exponential deviate
 254          *      u =  0,1 -uniform deviate
 255          *      t = (b,si)-float exponential (laplace) sample */
 256         e = randomExp();
 257         u = randomUniform();
 258         u = u + u - 1.0;
 259         if (u < 0.0)
 260             t = b - si * e;
 261         else
 262             t = b + si * e;
 263         /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */
 264         if (t >= -0.71874483771719) {
 265             /* Step 10:  calculation of v and quotient q */
 266             v = t / (s + s);
 267             if (fabs(v) <= 0.25)
 268                 q = q0 + 0.5 * t * t *
 269                                 ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
 270                                   + a2) * v + a1) * v;
 271             else
 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) */
 275             if (q > 0.0) {
 276                 w = expm1(q);
 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))
 280                     break;
 281             }
 282         }
 283     } /* for(;;) .. until  `t' is accepted */
 284     x = s + 0.5 * t;
 285     return scale * x * x;
 286 }
 287
 288 /**************************************************************************************************/
 289 //
 290 //      essentially swiped from http://en.wikipedia.org/wiki/Dirichlet_distribution#Random_number_generation
 291 //
 292 /**************************************************************************************************/
 293
 294 vector<float> RandomNumberGenerator::randomDirichlet(vector<float> alphas){
 295
 296         int nAlphas = (int)alphas.size();
 297         vector<float> dirs(nAlphas, 0.0000);
 298
 299         float sum = 0.0000;
 300         for(int i=0;i<nAlphas;i++){
 301                 dirs[i] = randomGamma(alphas[i]);
 302                 sum += dirs[i];
 303         }
 304
 305         for(int i=0;i<nAlphas;i++){
 306                 dirs[i] /= sum;
 307         }
 308
 309         return dirs;
 310
 311 }
 312
 313 /**************************************************************************************************/
 314
 315 #endif