htswanalysis/MACS/lib/gsl/gsl-1.11/randist/gamma.c

   1 /* randist/gamma.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <math.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_sf_gamma.h>
  24 #include <gsl/gsl_rng.h>
  25 #include <gsl/gsl_randist.h>
  26
  27 static double gamma_large (const gsl_rng * r, const double a);
  28 static double gamma_frac (const gsl_rng * r, const double a);
  29
  30 /* The Gamma distribution of order a>0 is defined by:
  31
  32    p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx
  33
  34    for x>0.  If X and Y are independent gamma-distributed random
  35    variables of order a1 and a2 with the same scale parameter b, then
  36    X+Y has gamma distribution of order a1+a2.
  37
  38    The algorithms below are from Knuth, vol 2, 2nd ed, p. 129. */
  39
  40 double
  41 gsl_ran_gamma_knuth (const gsl_rng * r, const double a, const double b)
  42 {
  43   /* assume a > 0 */
  44   unsigned int na = floor (a);
  45
  46   if (a == na)
  47     {
  48       return b * gsl_ran_gamma_int (r, na);
  49     }
  50   else if (na == 0)
  51     {
  52       return b * gamma_frac (r, a);
  53     }
  54   else
  55     {
  56       return b * (gsl_ran_gamma_int (r, na) + gamma_frac (r, a - na)) ;
  57     }
  58 }
  59
  60 double
  61 gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a)
  62 {
  63   if (a < 12)
  64     {
  65       unsigned int i;
  66       double prod = 1;
  67
  68       for (i = 0; i < a; i++)
  69         {
  70           prod *= gsl_rng_uniform_pos (r);
  71         }
  72
  73       /* Note: for 12 iterations we are safe against underflow, since
  74          the smallest positive random number is O(2^-32). This means
  75          the smallest possible product is 2^(-12*32) = 10^-116 which
  76          is within the range of double precision. */
  77
  78       return -log (prod);
  79     }
  80   else
  81     {
  82       return gamma_large (r, (double) a);
  83     }
  84 }
  85
  86 static double
  87 gamma_large (const gsl_rng * r, const double a)
  88 {
  89   /* Works only if a > 1, and is most efficient if a is large
  90
  91      This algorithm, reported in Knuth, is attributed to Ahrens.  A
  92      faster one, we are told, can be found in: J. H. Ahrens and
  93      U. Dieter, Computing 12 (1974) 223-246.  */
  94
  95   double sqa, x, y, v;
  96   sqa = sqrt (2 * a - 1);
  97   do
  98     {
  99       do
 100         {
 101           y = tan (M_PI * gsl_rng_uniform (r));
 102           x = sqa * y + a - 1;
 103         }
 104       while (x <= 0);
 105       v = gsl_rng_uniform (r);
 106     }
 107   while (v > (1 + y * y) * exp ((a - 1) * log (x / (a - 1)) - sqa * y));
 108
 109   return x;
 110 }
 111
 112 static double
 113 gamma_frac (const gsl_rng * r, const double a)
 114 {
 115   /* This is exercise 16 from Knuth; see page 135, and the solution is
 116      on page 551.  */
 117
 118   double p, q, x, u, v;
 119   p = M_E / (a + M_E);
 120   do
 121     {
 122       u = gsl_rng_uniform (r);
 123       v = gsl_rng_uniform_pos (r);
 124
 125       if (u < p)
 126         {
 127           x = exp ((1 / a) * log (v));
 128           q = exp (-x);
 129         }
 130       else
 131         {
 132           x = 1 - log (v);
 133           q = exp ((a - 1) * log (x));
 134         }
 135     }
 136   while (gsl_rng_uniform (r) >= q);
 137
 138   return x;
 139 }
 140
 141 double
 142 gsl_ran_gamma_pdf (const double x, const double a, const double b)
 143 {
 144   if (x < 0)
 145     {
 146       return 0 ;
 147     }
 148   else if (x == 0)
 149     {
 150       if (a == 1)
 151         return 1/b ;
 152       else
 153         return 0 ;
 154     }
 155   else if (a == 1)
 156     {
 157       return exp(-x/b)/b ;
 158     }
 159   else
 160     {
 161       double p;
 162       double lngamma = gsl_sf_lngamma (a);
 163       p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b;
 164       return p;
 165     }
 166 }
 167
 168
 169 /* New version based on Marsaglia and Tsang, "A Simple Method for
 170  * generating gamma variables", ACM Transactions on Mathematical
 171  * Software, Vol 26, No 3 (2000), p363-372.
 172  *
 173  * Implemented by J.D.Lamb@btinternet.com, minor modifications for GSL
 174  * by Brian Gough
 175  */
 176
 177 double
 178 gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b)
 179 {
 180   return gsl_ran_gamma (r, a, b);
 181 }
 182
 183 double
 184 gsl_ran_gamma (const gsl_rng * r, const double a, const double b)
 185 {
 186   /* assume a > 0 */
 187
 188   if (a < 1)
 189     {
 190       double u = gsl_rng_uniform_pos (r);
 191       return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
 192     }
 193
 194   {
 195     double x, v, u;
 196     double d = a - 1.0 / 3.0;
 197     double c = (1.0 / 3.0) / sqrt (d);
 198
 199     while (1)
 200       {
 201         do
 202           {
 203             x = gsl_ran_gaussian_ziggurat (r, 1.0);
 204             v = 1.0 + c * x;
 205           }
 206         while (v <= 0);
 207
 208         v = v * v * v;
 209         u = gsl_rng_uniform_pos (r);
 210
 211         if (u < 1 - 0.0331 * x * x * x * x)
 212           break;
 213
 214         if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
 215           break;
 216       }
 217
 218     return b * d * v;
 219   }
 220 }