htswanalysis/MACS/lib/gsl/gsl-1.11/cdf/gammainv.c

   1 /* cdf/gammainv.c
   2  *
   3  * Copyright (C) 2003, 2007 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <math.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_randist.h>
  25 #include <gsl/gsl_sf_gamma.h>
  26
  27 #include <stdio.h>
  28
  29 double
  30 gsl_cdf_gamma_Pinv (const double P, const double a, const double b)
  31 {
  32   double x;
  33
  34   if (P == 1.0)
  35     {
  36       return GSL_POSINF;
  37     }
  38   else if (P == 0.0)
  39     {
  40       return 0.0;
  41     }
  42
  43   /* Consider, small, large and intermediate cases separately.  The
  44      boundaries at 0.05 and 0.95 have not been optimised, but seem ok
  45      for an initial approximation. */
  46
  47   if (P < 0.05)
  48     {
  49       double x0 = exp ((gsl_sf_lngamma (a) + log (P)) / a);
  50       x = x0;
  51     }
  52   else if (P > 0.95)
  53     {
  54       double x0 = -log1p (-P) + gsl_sf_lngamma (a);
  55       x = x0;
  56     }
  57   else
  58     {
  59       double xg = gsl_cdf_ugaussian_Pinv (P);
  60       double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;
  61       x = x0;
  62     }
  63
  64   /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
  65      to an improved value of x (Abramowitz & Stegun, 3.6.6)
  66
  67      where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
  68    */
  69
  70   {
  71     double lambda, dP, phi;
  72     unsigned int n = 0;
  73
  74   start:
  75     dP = P - gsl_cdf_gamma_P (x, a, 1.0);
  76     phi = gsl_ran_gamma_pdf (x, a, 1.0);
  77
  78     if (dP == 0.0 || n++ > 32)
  79       goto end;
  80
  81     lambda = dP / GSL_MAX (2 * fabs (dP / x), phi);
  82
  83     {
  84       double step0 = lambda;
  85       double step1 = -((a - 1) / x - 1) * lambda * lambda / 4.0;
  86
  87       double step = step0;
  88       if (fabs (step1) < fabs (step0))
  89         step += step1;
  90
  91       if (x + step > 0)
  92         x += step;
  93       else
  94         {
  95           x /= 2.0;
  96         }
  97
  98       if (fabs (step0) > 1e-10 * x)
  99         goto start;
 100     }
 101
 102   end:
 103     if (fabs(dP) > GSL_SQRT_DBL_EPSILON * P)
 104       {
 105         GSL_ERROR_VAL("inverse failed to converge", GSL_EFAILED, GSL_NAN);
 106       }
 107
 108     return b * x;
 109   }
 110 }
 111
 112 double
 113 gsl_cdf_gamma_Qinv (const double Q, const double a, const double b)
 114 {
 115   double x;
 116
 117   if (Q == 1.0)
 118     {
 119       return 0.0;
 120     }
 121   else if (Q == 0.0)
 122     {
 123       return GSL_POSINF;
 124     }
 125
 126   /* Consider, small, large and intermediate cases separately.  The
 127      boundaries at 0.05 and 0.95 have not been optimised, but seem ok
 128      for an initial approximation. */
 129
 130   if (Q < 0.05)
 131     {
 132       double x0 = -log (Q) + gsl_sf_lngamma (a);
 133       x = x0;
 134     }
 135   else if (Q > 0.95)
 136     {
 137       double x0 = exp ((gsl_sf_lngamma (a) + log1p (-Q)) / a);
 138       x = x0;
 139     }
 140   else
 141     {
 142       double xg = gsl_cdf_ugaussian_Qinv (Q);
 143       double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;
 144       x = x0;
 145     }
 146
 147   /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
 148      to an improved value of x (Abramowitz & Stegun, 3.6.6)
 149
 150      where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
 151    */
 152
 153   {
 154     double lambda, dQ, phi;
 155     unsigned int n = 0;
 156
 157   start:
 158     dQ = Q - gsl_cdf_gamma_Q (x, a, 1.0);
 159     phi = gsl_ran_gamma_pdf (x, a, 1.0);
 160
 161     if (dQ == 0.0 || n++ > 32)
 162       goto end;
 163
 164     lambda = -dQ / GSL_MAX (2 * fabs (dQ / x), phi);
 165
 166     {
 167       double step0 = lambda;
 168       double step1 = -((a - 1) / x - 1) * lambda * lambda / 4.0;
 169
 170       double step = step0;
 171       if (fabs (step1) < fabs (step0))
 172         step += step1;
 173
 174       if (x + step > 0)
 175         x += step;
 176       else
 177         {
 178           x /= 2.0;
 179         }
 180
 181       if (fabs (step0) > 1e-10 * x)
 182         goto start;
 183     }
 184
 185   }
 186
 187 end:
 188   return b * x;
 189 }