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

   1 /* cdf/tdistinv.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 2002 Jason H. Stover.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
  19  */
  20
  21 #include <config.h>
  22 #include <math.h>
  23 #include <gsl/gsl_cdf.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_randist.h>
  26 #include <gsl/gsl_sf_gamma.h>
  27
  28 #include <stdio.h>
  29
  30 static double
  31 inv_cornish_fisher (double z, double nu)
  32 {
  33   double a = 1 / (nu - 0.5);
  34   double b = 48.0 / (a * a);
  35
  36   double cf1 = z * (3 + z * z);
  37   double cf2 = z * (945 + z * z * (360 + z * z * (63 + z * z * 4)));
  38
  39   double y = z - cf1 / b + cf2 / (10 * b * b);
  40
  41   double t = GSL_SIGN (z) * sqrt (nu * expm1 (a * y * y));
  42
  43   return t;
  44 }
  45
  46
  47 double
  48 gsl_cdf_tdist_Pinv (const double P, const double nu)
  49 {
  50   double x, ptail;
  51
  52   if (P == 1.0)
  53     {
  54       return GSL_POSINF;
  55     }
  56   else if (P == 0.0)
  57     {
  58       return GSL_NEGINF;
  59     }
  60
  61   if (nu == 1.0)
  62     {
  63       x = tan (M_PI * (P - 0.5));
  64     }
  65   else if (nu == 2.0)
  66     {
  67       double a = 2 * P - 1;
  68       x = a / sqrt (2 * (1 - a * a));
  69     }
  70
  71   ptail = (P < 0.5) ? P : 1 - P;
  72
  73   if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2))
  74     {
  75       double xg = gsl_cdf_ugaussian_Pinv (P);
  76       x = inv_cornish_fisher (xg, nu);
  77     }
  78   else
  79     {
  80       /* Use an asymptotic expansion of the tail of integral */
  81
  82       double beta = gsl_sf_beta (0.5, nu / 2);
  83
  84       if (P < 0.5)
  85         {
  86           x = -sqrt (nu) * pow (beta * nu * P, -1.0 / nu);
  87         }
  88       else
  89         {
  90           x = sqrt (nu) * pow (beta * nu * (1 - P), -1.0 / nu);
  91         }
  92
  93       /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
  94          for higher order terms. This avoids overestimating x, which
  95          makes the iteration unstable due to the rapidly decreasing
  96          tails of the distribution. */
  97
  98       x /= sqrt (1 + nu / (x * x));
  99     }
 100
 101   {
 102     double dP, phi;
 103     unsigned int n = 0;
 104
 105   start:
 106     dP = P - gsl_cdf_tdist_P (x, nu);
 107     phi = gsl_ran_tdist_pdf (x, nu);
 108
 109     if (dP == 0.0 || n++ > 32)
 110       goto end;
 111
 112     {
 113       double lambda = dP / phi;
 114       double step0 = lambda;
 115       double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);
 116
 117       double step = step0;
 118
 119       if (fabs (step1) < fabs (step0))
 120         {
 121           step += step1;
 122         }
 123
 124       if (P > 0.5 && x + step < 0)
 125         x /= 2;
 126       else if (P < 0.5 && x + step > 0)
 127         x /= 2;
 128       else
 129         x += step;
 130
 131       if (fabs (step) > 1e-10 * fabs (x))
 132         goto start;
 133     }
 134
 135   end:
 136     if (fabs(dP) > GSL_SQRT_DBL_EPSILON * P)
 137       {
 138         GSL_ERROR_VAL("inverse failed to converge", GSL_EFAILED, GSL_NAN);
 139       }
 140
 141     return x;
 142   }
 143 }
 144
 145 double
 146 gsl_cdf_tdist_Qinv (const double Q, const double nu)
 147 {
 148   double x, qtail;
 149
 150   if (Q == 0.0)
 151     {
 152       return GSL_POSINF;
 153     }
 154   else if (Q == 1.0)
 155     {
 156       return GSL_NEGINF;
 157     }
 158
 159   if (nu == 1.0)
 160     {
 161       x = tan (M_PI * (0.5 - Q));
 162     }
 163   else if (nu == 2.0)
 164     {
 165       double a = 2 * (1 - Q) - 1;
 166       x = a / sqrt (2 * (1 - a * a));
 167     }
 168
 169   qtail = (Q < 0.5) ? Q : 1 - Q;
 170
 171   if (sqrt (M_PI * nu / 2) * qtail > pow (0.05, nu / 2))
 172     {
 173       double xg = gsl_cdf_ugaussian_Qinv (Q);
 174       x = inv_cornish_fisher (xg, nu);
 175     }
 176   else
 177     {
 178       /* Use an asymptotic expansion of the tail of integral */
 179
 180       double beta = gsl_sf_beta (0.5, nu / 2);
 181
 182       if (Q < 0.5)
 183         {
 184           x = sqrt (nu) * pow (beta * nu * Q, -1.0 / nu);
 185         }
 186       else
 187         {
 188           x = -sqrt (nu) * pow (beta * nu * (1 - Q), -1.0 / nu);
 189         }
 190
 191       /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
 192          for higher order terms. This avoids overestimating x, which
 193          makes the iteration unstable due to the rapidly decreasing
 194          tails of the distribution. */
 195
 196       x /= sqrt (1 + nu / (x * x));
 197     }
 198
 199   {
 200     double dQ, phi;
 201     unsigned int n = 0;
 202
 203   start:
 204     dQ = Q - gsl_cdf_tdist_Q (x, nu);
 205     phi = gsl_ran_tdist_pdf (x, nu);
 206
 207     if (dQ == 0.0 || n++ > 32)
 208       goto end;
 209
 210     {
 211       double lambda = - dQ / phi;
 212       double step0 = lambda;
 213       double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);
 214
 215       double step = step0;
 216
 217       if (fabs (step1) < fabs (step0))
 218         {
 219           step += step1;
 220         }
 221
 222       if (Q < 0.5 && x + step < 0)
 223         x /= 2;
 224       else if (Q > 0.5 && x + step > 0)
 225         x /= 2;
 226       else
 227         x += step;
 228
 229       if (fabs (step) > 1e-10 * fabs (x))
 230         goto start;
 231     }
 232   }
 233
 234 end:
 235
 236   return x;
 237 }