/[gsl]/gsl/cdf/gammainv.c
ViewVC logotype

Diff of /gsl/cdf/gammainv.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.2 by bjg, Sun Jun 29 17:31:38 2003 UTC revision 1.1.2.3 by bjg, Tue Jul 1 15:02:03 2003 UTC
# Line 21  Line 21 
21  #include <math.h>  #include <math.h>
22  #include <gsl/gsl_cdf.h>  #include <gsl/gsl_cdf.h>
23  #include <gsl/gsl_math.h>  #include <gsl/gsl_math.h>
24    #include <gsl/gsl_randist.h>
25  #include <gsl/gsl_sf_gamma.h>  #include <gsl/gsl_sf_gamma.h>
26    
27  #include <stdio.h>  #include <stdio.h>
# Line 28  Line 29 
29  double  double
30  gsl_cdf_gamma_Pinv (double p, double a, double b)  gsl_cdf_gamma_Pinv (double p, double a, double b)
31  {  {
32    double x, xc;    double x;
33    
34    if (p == 1.0)    if (p == 1.0)
35      {      {
# Line 39  gsl_cdf_gamma_Pinv (double p, double a, Line 40  gsl_cdf_gamma_Pinv (double p, double a,
40        return 0.0;        return 0.0;
41      }      }
42    
43      if (p < 0.05)
44        {
45          double x0 = exp((lgamma(b) + log(p))/b);
46          x = x0;
47        }
48      else if (p > 0.95)
49        {
50          double x0 = -log(1-p) + lgamma(b);
51          x = x0;
52        }
53      else
54        {
55          double xg = gsl_cdf_ugaussian_Pinv (p);
56          double x0 = (xg < -sqrt(b)) ? b : sqrt (b) * xg + b;
57          x = x0;
58        }
59    
60      /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
61         to an improved value of x (Abramowitz & Stegun, 3.6.6)
62    
63         where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
64      */
65    
66      {
67        double lambda, dp, phi;
68    
69      start:
70        dp = p - gsl_cdf_gamma_P (x, 1.0, b) ;
71        phi = gsl_ran_gamma_pdf (x, 1.0, b);
72    
73        if (dp == 0.0)
74          goto end;
75    
76        lambda = dp / GSL_MAX(2*fabs(dp/x), phi);
77    
78        {
79          double step0 = lambda;
80          double step1 =  -((b-1)/x - 1)*lambda*lambda/4.0;
81        
82          double step = step0;
83          if (fabs(step1) < fabs(step0))
84            step += step1;
85          
86          if (x + step > 0)
87            x += step;
88          else
89            {
90              x /= 2.0;
91            }
92          
93          if (fabs(step0) > 1e-10*x)
94              goto start;
95        }
96        
97      }
98    
99     end:
100    return a * x;    return a * x;
101  }  }
102    

Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.1.2.3

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26