/[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.3 by bjg, Tue Jul 1 15:02:03 2003 UTC revision 1.1.2.4 by bjg, Sat Jul 12 20:57:55 2003 UTC
# Line 114  gsl_cdf_gamma_Qinv (double p, double a, Line 114  gsl_cdf_gamma_Qinv (double p, double a,
114        return GSL_POSINF;        return GSL_POSINF;
115      }      }
116    
117      if (p < 0.05)
118        {
119          double x0 = -log(p) + lgamma(b);
120          x = x0;
121        }
122      else if (p > 0.95)
123        {
124          double x0 = exp((lgamma(b) + log1p(-p))/b);
125          x = x0;
126        }
127      else
128        {
129          double xg = gsl_cdf_ugaussian_Qinv (p);
130          double x0 = (xg < -sqrt(b)) ? b : sqrt (b) * xg + b;
131          x = x0;
132        }
133    
134      /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
135         to an improved value of x (Abramowitz & Stegun, 3.6.6)
136    
137         where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
138      */
139    
140      {
141        double lambda, dp, phi;
142    
143      start:
144        dp = -(p - gsl_cdf_gamma_Q (x, 1.0, b)) ;
145        phi = gsl_ran_gamma_pdf (x, 1.0, b);
146    
147        if (dp == 0.0)
148          goto end;
149    
150        lambda = dp / GSL_MAX(2*fabs(dp/x), phi);
151    
152        {
153          double step0 = lambda;
154          double step1 =  -((b-1)/x - 1)*lambda*lambda/4.0;
155        
156          double step = step0;
157          if (fabs(step1) < fabs(step0))
158            step += step1;
159          
160          if (x + step > 0)
161            x += step;
162          else
163            {
164              x /= 2.0;
165            }
166          
167          if (fabs(step0) > 1e-10*x)
168              goto start;
169        }
170        
171      }
172    
173    return x;   end:
174      return a * x;
175  }  }

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

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