/[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.4 by bjg, Sat Jul 12 20:57:55 2003 UTC revision 1.1.2.5 by bjg, Sat Jul 12 22:00:18 2003 UTC
# Line 40  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      /* 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)    if (p < 0.05)
48      {      {
49        double x0 = exp((lgamma(b) + log(p))/b);        double x0 = exp ((lgamma (b) + log (p)) / b);
50        x = x0;        x = x0;
51      }      }
52    else if (p > 0.95)    else if (p > 0.95)
53      {      {
54        double x0 = -log(1-p) + lgamma(b);        double x0 = -log (1 - p) + lgamma (b);
55        x = x0;        x = x0;
56      }      }
57    else    else
58      {      {
59        double xg = gsl_cdf_ugaussian_Pinv (p);        double xg = gsl_cdf_ugaussian_Pinv (p);
60        double x0 = (xg < -sqrt(b)) ? b : sqrt (b) * xg + b;        double x0 = (xg < -sqrt (b)) ? b : sqrt (b) * xg + b;
61        x = x0;        x = x0;
62      }      }
63    
# Line 61  gsl_cdf_gamma_Pinv (double p, double a, Line 65  gsl_cdf_gamma_Pinv (double p, double a,
65       to an improved value of x (Abramowitz & Stegun, 3.6.6)       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.       where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
68    */     */
69    
70    {    {
71      double lambda, dp, phi;      double lambda, dp, phi;
72    
73    start:    start:
74      dp = p - gsl_cdf_gamma_P (x, 1.0, b) ;      dp = p - gsl_cdf_gamma_P (x, 1.0, b);
75      phi = gsl_ran_gamma_pdf (x, 1.0, b);      phi = gsl_ran_gamma_pdf (x, b, 1.0);
76    
77      if (dp == 0.0)      if (dp == 0.0)
78        goto end;        goto end;
79    
80      lambda = dp / GSL_MAX(2*fabs(dp/x), phi);      lambda = dp / GSL_MAX (2 * fabs (dp / x), phi);
81    
82      {      {
83        double step0 = lambda;        double step0 = lambda;
84        double step1 =  -((b-1)/x - 1)*lambda*lambda/4.0;        double step1 = -((b - 1) / x - 1) * lambda * lambda / 4.0;
85        
86        double step = step0;        double step = step0;
87        if (fabs(step1) < fabs(step0))        if (fabs (step1) < fabs (step0))
88          step += step1;          step += step1;
89          
90        if (x + step > 0)        if (x + step > 0)
91          x += step;          x += step;
92        else        else
93          {          {
94            x /= 2.0;            x /= 2.0;
95          }          }
96          
97        if (fabs(step0) > 1e-10*x)        if (fabs (step0) > 1e-10 * x)
98            goto start;          goto start;
99      }      }
100        
101    }    }
102    
103   end:  end:
104    return a * x;    return a * x;
105  }  }
106    
# Line 114  gsl_cdf_gamma_Qinv (double p, double a, Line 118  gsl_cdf_gamma_Qinv (double p, double a,
118        return GSL_POSINF;        return GSL_POSINF;
119      }      }
120    
121      /* Consider, small, large and intermediate cases separately.  The
122         boundaries at 0.05 and 0.95 have not been optimised, but seem ok
123         for an initial approximation. */
124    
125    if (p < 0.05)    if (p < 0.05)
126      {      {
127        double x0 = -log(p) + lgamma(b);        double x0 = -log (p) + lgamma (b);
128        x = x0;        x = x0;
129      }      }
130    else if (p > 0.95)    else if (p > 0.95)
131      {      {
132        double x0 = exp((lgamma(b) + log1p(-p))/b);        double x0 = exp ((lgamma (b) + log1p (-p)) / b);
133        x = x0;        x = x0;
134      }      }
135    else    else
136      {      {
137        double xg = gsl_cdf_ugaussian_Qinv (p);        double xg = gsl_cdf_ugaussian_Qinv (p);
138        double x0 = (xg < -sqrt(b)) ? b : sqrt (b) * xg + b;        double x0 = (xg < -sqrt (b)) ? b : sqrt (b) * xg + b;
139        x = x0;        x = x0;
140      }      }
141    
# Line 135  gsl_cdf_gamma_Qinv (double p, double a, Line 143  gsl_cdf_gamma_Qinv (double p, double a,
143       to an improved value of x (Abramowitz & Stegun, 3.6.6)       to an improved value of x (Abramowitz & Stegun, 3.6.6)
144    
145       where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.       where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
146    */     */
147    
148    {    {
149      double lambda, dp, phi;      double lambda, dp, phi;
150    
151    start:    start:
152      dp = -(p - gsl_cdf_gamma_Q (x, 1.0, b)) ;      dp = -(p - gsl_cdf_gamma_Q (x, 1.0, b));
153      phi = gsl_ran_gamma_pdf (x, 1.0, b);      phi = gsl_ran_gamma_pdf (x, 1.0, b);
154    
155      if (dp == 0.0)      if (dp == 0.0)
156        goto end;        goto end;
157    
158      lambda = dp / GSL_MAX(2*fabs(dp/x), phi);      lambda = dp / GSL_MAX (2 * fabs (dp / x), phi);
159    
160      {      {
161        double step0 = lambda;        double step0 = lambda;
162        double step1 =  -((b-1)/x - 1)*lambda*lambda/4.0;        double step1 = -((b - 1) / x - 1) * lambda * lambda / 4.0;
163        
164        double step = step0;        double step = step0;
165        if (fabs(step1) < fabs(step0))        if (fabs (step1) < fabs (step0))
166          step += step1;          step += step1;
167          
168        if (x + step > 0)        if (x + step > 0)
169          x += step;          x += step;
170        else        else
171          {          {
172            x /= 2.0;            x /= 2.0;
173          }          }
174          
175        if (fabs(step0) > 1e-10*x)        if (fabs (step0) > 1e-10 * x)
176            goto start;          goto start;
177      }      }
178        
179    }    }
180    
181   end:  end:
182    return a * x;    return a * x;
183  }  }

Legend:
Removed from v.1.1.2.4  
changed lines
  Added in v.1.1.2.5

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