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

   1 /* cdf/betainv.c
   2  *
   3  * Copyright (C) 2004 Free Software Foundation, Inc.
   4  * Copyright (C) 2006, 2007 Brian Gough
   5  * Written by Jason H. Stover.
   6  * Modified for GSL by Brian Gough
   7  *
   8  * This program is free software; you can redistribute it and/or modify
   9  * it under the terms of the GNU General Public License as published by
  10  * the Free Software Foundation; either version 3 of the License, or (at
  11  * your option) any later version.
  12  *
  13  * This program is distributed in the hope that it will be useful, but
  14  * WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License
  19  * along with this program; if not, write to the Free Software
  20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
  21  */
  22
  23 /*
  24  * Invert the Beta distribution.
  25  *
  26  * References:
  27  *
  28  * Roger W. Abernathy and Robert P. Smith. "Applying Series Expansion
  29  * to the Inverse Beta Distribution to Find Percentiles of the
  30  * F-Distribution," ACM Transactions on Mathematical Software, volume
  31  * 19, number 4, December 1993, pages 474-480.
  32  *
  33  * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions of a
  34  * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39,
  35  * number 8, August 1968, pages 1264-1273.
  36  */
  37
  38 #include <config.h>
  39 #include <math.h>
  40 #include <gsl/gsl_math.h>
  41 #include <gsl/gsl_errno.h>
  42 #include <gsl/gsl_sf_gamma.h>
  43 #include <gsl/gsl_cdf.h>
  44 #include <gsl/gsl_randist.h>
  45
  46 #include "error.h"
  47
  48 static double
  49 bisect (double x, double P, double a, double b, double xtol, double Ptol)
  50 {
  51   double x0 = 0, x1 = 1, Px;
  52
  53   while (fabs(x1 - x0) > xtol) {
  54     Px = gsl_cdf_beta_P (x, a, b);
  55     if (fabs(Px - P) < Ptol) {
  56       /* return as soon as approximation is good enough, including on
  57          the first iteration */
  58       return x;
  59     } else if (Px < P) {
  60       x0 = x;
  61     } else if (Px > P) {
  62       x1 = x;
  63     }
  64     x = 0.5 * (x0 + x1);
  65   }
  66   return x;
  67 }
  68
  69
  70 double
  71 gsl_cdf_beta_Pinv (const double P, const double a, const double b)
  72 {
  73   double x, mean;
  74
  75   if (P < 0.0 || P > 1.0)
  76     {
  77       CDF_ERROR ("P must be in range 0 < P < 1", GSL_EDOM);
  78     }
  79
  80   if (a < 0.0)
  81     {
  82       CDF_ERROR ("a < 0", GSL_EDOM);
  83     }
  84
  85   if (b < 0.0)
  86     {
  87       CDF_ERROR ("b < 0", GSL_EDOM);
  88     }
  89
  90   if (P == 0.0)
  91     {
  92       return 0.0;
  93     }
  94
  95   if (P == 1.0)
  96     {
  97       return 1.0;
  98     }
  99
 100   if (P > 0.5)
 101     {
 102       return gsl_cdf_beta_Qinv (1 - P, a, b);
 103     }
 104
 105   mean = a / (a + b);
 106
 107   if (P < 0.1)
 108     {
 109       /* small x */
 110
 111       double lg_ab = gsl_sf_lngamma (a + b);
 112       double lg_a = gsl_sf_lngamma (a);
 113       double lg_b = gsl_sf_lngamma (b);
 114
 115       double lx = (log (a) + lg_a + lg_b - lg_ab + log (P)) / a;
 116       if (lx <= 0) {
 117         x = exp (lx);             /* first approximation */
 118         x *= pow (1 - x, -(b - 1) / a);   /* second approximation */
 119       } else {
 120         x = mean;
 121       }
 122
 123       if (x > mean)
 124         x = mean;
 125     }
 126   else
 127     {
 128       /* Use expected value as first guess */
 129       x = mean;
 130     }
 131
 132   /* Do bisection to get closer */
 133   x = bisect (x, P, a, b, 0.01, 0.01);
 134
 135   {
 136     double lambda, dP, phi;
 137     unsigned int n = 0;
 138
 139   start:
 140     dP = P - gsl_cdf_beta_P (x, a, b);
 141     phi = gsl_ran_beta_pdf (x, a, b);
 142
 143     if (dP == 0.0 || n++ > 64)
 144       goto end;
 145
 146     lambda = dP / GSL_MAX (2 * fabs (dP / x), phi);
 147
 148     {
 149       double step0 = lambda;
 150       double step1 = -((a - 1) / x - (b - 1) / (1 - x)) * lambda * lambda / 2;
 151
 152       double step = step0;
 153
 154       if (fabs (step1) < fabs (step0))
 155         {
 156           step += step1;
 157         }
 158       else
 159         {
 160           /* scale back step to a reasonable size when too large */
 161           step *= 2 * fabs (step0 / step1);
 162         };
 163
 164       if (x + step > 0 && x + step < 1)
 165         {
 166           x += step;
 167         }
 168       else
 169         {
 170           x = sqrt (x) * sqrt (mean);   /* try a new starting point */
 171         }
 172
 173       if (fabs (step0) > 1e-10 * x)
 174         goto start;
 175     }
 176
 177   end:
 178
 179     if (fabs(dP) > GSL_SQRT_DBL_EPSILON * P)
 180       {
 181         GSL_ERROR_VAL("inverse failed to converge", GSL_EFAILED, GSL_NAN);
 182       }
 183
 184     return x;
 185   }
 186 }
 187
 188 double
 189 gsl_cdf_beta_Qinv (const double Q, const double a, const double b)
 190 {
 191
 192   if (Q < 0.0 || Q > 1.0)
 193     {
 194       CDF_ERROR ("Q must be inside range 0 < Q < 1", GSL_EDOM);
 195     }
 196
 197   if (a < 0.0)
 198     {
 199       CDF_ERROR ("a < 0", GSL_EDOM);
 200     }
 201
 202   if (b < 0.0)
 203     {
 204       CDF_ERROR ("b < 0", GSL_EDOM);
 205     }
 206
 207   if (Q == 0.0)
 208     {
 209       return 1.0;
 210     }
 211
 212   if (Q == 1.0)
 213     {
 214       return 0.0;
 215     }
 216
 217   if (Q > 0.5)
 218     {
 219       return gsl_cdf_beta_Pinv (1 - Q, a, b);
 220     }
 221   else
 222     {
 223       return 1 - gsl_cdf_beta_Pinv (Q, b, a);
 224     };
 225 }