Go back to Random distributions

<project-file type=“source”/> <content> /* this file implements a gamma distribution generator. It is shamelessly

 copied from the gsl (Gnu Science Library) code, and adapted to my own 
 need. The original copyright message is provided below. */

/* randist/gamma.c * * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */

#include “genrand2.h” #include <math.h>

static float gamma_large ( const float a); static float gamma_frac (const float a); static float gsl_ran_gamma_int (const unsigned int a);

/* The Gamma distribution of order a>0 is defined by:

 p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx

 for x>0.  If X and Y are independent gamma-distributed random
 variables of order a1 and a2 with the same scale parameter b, then
 X+Y has gamma distribution of order a1+a2.

 The algorithms below are from Knuth, vol 2, 2nd ed, p. 129. */

float gammaRandom (const float a, const float b) {

/* assume a > 0 */
unsigned int na = (int)floor (a);

if (a == na)
  {
    return b * gsl_ran_gamma_int (na);
  }
else if (na == 0)
  {
    return b * gamma_frac (a);
  }
else
  {
    return b * (gsl_ran_gamma_int (na) + gamma_frac (a - na)) ;
  }

}

float gsl_ran_gamma_int (const unsigned int a) {

if (a < 12)
  {
    unsigned int i;
    float prod = 1,r;

    for (i = 0; i < a; i++)
{
        do {
    r=genrand2d();
  } while (r==0);
  prod *= r;
}

    /* Note: for 12 iterations we are safe against underflow, since
 the smallest positive random number is O(2^-32). This means
 the smallest possible product is 2^(-12*32) = 10^-116 which
 is within the range of float precision. */

    return -log (prod);
  }
else
  {
    return gamma_large ((float) a);
  }

}

static float gamma_large (const float a) {

/* Works only if a > 1, and is most efficient if a is large

   This algorithm, reported in Knuth, is attributed to Ahrens.  A
   faster one, we are told, can be found in: J. H. Ahrens and
   U. Dieter, Computing 12 (1974) 223-246.  */

float sqa, x, y, v;
sqa = sqrt (2 * a - 1);
do
  {
    do
{
  y = tan (M_PI * genrand2d());
  x = sqa * y + a - 1;
}
    while (x <= 0);
    v = genrand2d();
  }
while (v > (1 + y * y) * exp ((a - 1) * log (x / (a - 1)) - sqa * y));

return x;

}

static float gamma_frac (const float a) {

/* This is exercise 16 from Knuth; see page 135, and the solution is
   on page 551.  */

float p, q, x, u, v;
p = M_E / (a + M_E);
do
  {
    u = genrand2d();
    do {
       v = genrand2d();
    } while (v==0);

    if (u < p)
{
  x = exp ((1 / a) * log (v));
  q = exp (-x);
}
    else
{
  x = 1 - log (v);
  q = exp ((a - 1) * log (x));
}
  }
while (genrand2d() >= q);

return x;

} </content> <use name=“genrand2.h”/>