diff options
Diffstat (limited to 'gsl-1.9/randist/gamma.c')
| -rw-r--r-- | gsl-1.9/randist/gamma.c | 220 |
1 files changed, 220 insertions, 0 deletions
diff --git a/gsl-1.9/randist/gamma.c b/gsl-1.9/randist/gamma.c new file mode 100644 index 0000000..37d167c --- /dev/null +++ b/gsl-1.9/randist/gamma.c @@ -0,0 +1,220 @@ +/* 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_sf_gamma.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_randist.h> + +static double gamma_large (const gsl_rng * r, const double a); +static double gamma_frac (const gsl_rng * r, const double 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. */ + +double +gsl_ran_gamma_knuth (const gsl_rng * r, const double a, const double b) +{ + /* assume a > 0 */ + unsigned int na = floor (a); + + if (a == na) + { + return b * gsl_ran_gamma_int (r, na); + } + else if (na == 0) + { + return b * gamma_frac (r, a); + } + else + { + return b * (gsl_ran_gamma_int (r, na) + gamma_frac (r, a - na)) ; + } +} + +double +gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a) +{ + if (a < 12) + { + unsigned int i; + double prod = 1; + + for (i = 0; i < a; i++) + { + prod *= gsl_rng_uniform_pos (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 double precision. */ + + return -log (prod); + } + else + { + return gamma_large (r, (double) a); + } +} + +static double +gamma_large (const gsl_rng * r, const double 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. */ + + double sqa, x, y, v; + sqa = sqrt (2 * a - 1); + do + { + do + { + y = tan (M_PI * gsl_rng_uniform (r)); + x = sqa * y + a - 1; + } + while (x <= 0); + v = gsl_rng_uniform (r); + } + while (v > (1 + y * y) * exp ((a - 1) * log (x / (a - 1)) - sqa * y)); + + return x; +} + +static double +gamma_frac (const gsl_rng * r, const double a) +{ + /* This is exercise 16 from Knuth; see page 135, and the solution is + on page 551. */ + + double p, q, x, u, v; + p = M_E / (a + M_E); + do + { + u = gsl_rng_uniform (r); + v = gsl_rng_uniform_pos (r); + + if (u < p) + { + x = exp ((1 / a) * log (v)); + q = exp (-x); + } + else + { + x = 1 - log (v); + q = exp ((a - 1) * log (x)); + } + } + while (gsl_rng_uniform (r) >= q); + + return x; +} + +double +gsl_ran_gamma_pdf (const double x, const double a, const double b) +{ + if (x < 0) + { + return 0 ; + } + else if (x == 0) + { + if (a == 1) + return 1/b ; + else + return 0 ; + } + else if (a == 1) + { + return exp(-x/b)/b ; + } + else + { + double p; + double lngamma = gsl_sf_lngamma (a); + p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b; + return p; + } +} + + +/* New version based on Marsaglia and Tsang, "A Simple Method for + * generating gamma variables", ACM Transactions on Mathematical + * Software, Vol 26, No 3 (2000), p363-372. + * + * Implemented by J.D.Lamb@btinternet.com, minor modifications for GSL + * by Brian Gough + */ + +double +gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b) +{ + return gsl_ran_gamma (r, a, b); +} + +double +gsl_ran_gamma (const gsl_rng * r, const double a, const double b) +{ + /* assume a > 0 */ + + if (a < 1) + { + double u = gsl_rng_uniform_pos (r); + return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a); + } + + { + double x, v, u; + double d = a - 1.0 / 3.0; + double c = (1.0 / 3.0) / sqrt (d); + + while (1) + { + do + { + x = gsl_ran_gaussian_ziggurat (r, 1.0); + v = 1.0 + c * x; + } + while (v <= 0); + + v = v * v * v; + u = gsl_rng_uniform_pos (r); + + if (u < 1 - 0.0331 * x * x * x * x) + break; + + if (log (u) < 0.5 * x * x + d * (1 - v + log (v))) + break; + } + + return b * d * v; + } +} |