/[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.7 by bjg, Tue Jul 15 17:56:47 2003 UTC revision 1.1.2.8 by bjg, Sun Jul 20 15:02:50 2003 UTC
# Line 27  Line 27 
27  #include <stdio.h>  #include <stdio.h>
28    
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;    double x;
33    
34    if (p == 1.0)    if (P == 1.0)
35      {      {
36        return GSL_POSINF;        return GSL_POSINF;
37      }      }
38    else if (p == 0.0)    else if (P == 0.0)
39      {      {
40        return 0.0;        return 0.0;
41      }      }
# Line 44  gsl_cdf_gamma_Pinv (double p, double a, Line 44  gsl_cdf_gamma_Pinv (double p, double a,
44       boundaries at 0.05 and 0.95 have not been optimised, but seem ok       boundaries at 0.05 and 0.95 have not been optimised, but seem ok
45       for an initial approximation. */       for an initial approximation. */
46    
47    if (p < 0.05)    if (P < 0.05)
48      {      {
49        double x0 = exp ((lgamma (a) + log (p)) / a);        double x0 = exp ((gsl_sf_lngamma (a) + log (P)) / a);
50        x = x0;        x = x0;
51      }      }
52    else if (p > 0.95)    else if (P > 0.95)
53      {      {
54        double x0 = -log1p (-p) + lgamma (a);        double x0 = -log1p (-P) + gsl_sf_lngamma (a);
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 (a)) ? a : sqrt (a) * xg + a;        double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;
61        x = x0;        x = x0;
62      }      }
# Line 68  gsl_cdf_gamma_Pinv (double p, double a, Line 68  gsl_cdf_gamma_Pinv (double p, double a,
68     */     */
69    
70    {    {
71      double lambda, dp, phi;      double lambda, dP, phi;
72    
73    start:    start:
74      dp = p - gsl_cdf_gamma_P (x, a, 1.0);      dP = P - gsl_cdf_gamma_P (x, a, 1.0);
75      phi = gsl_ran_gamma_pdf (x, a, 1.0);      phi = gsl_ran_gamma_pdf (x, a, 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;
# Line 105  end: Line 105  end:
105  }  }
106    
107  double  double
108  gsl_cdf_gamma_Qinv (double p, double a, double b)  gsl_cdf_gamma_Qinv (double Q, double a, double b)
109  {  {
110    double x;    double x;
111    
112    if (p == 1.0)    if (Q == 1.0)
113      {      {
114        return 0.0;        return 0.0;
115      }      }
116    else if (p == 0.0)    else if (Q == 0.0)
117      {      {
118        return GSL_POSINF;        return GSL_POSINF;
119      }      }
# Line 122  gsl_cdf_gamma_Qinv (double p, double a, Line 122  gsl_cdf_gamma_Qinv (double p, double a,
122       boundaries at 0.05 and 0.95 have not been optimised, but seem ok       boundaries at 0.05 and 0.95 have not been optimised, but seem ok
123       for an initial approximation. */       for an initial approximation. */
124    
125    if (p < 0.05)    if (Q < 0.05)
126      {      {
127        double x0 = -log (p) + lgamma (a);        double x0 = -log (Q) + gsl_sf_lngamma (a);
128        x = x0;        x = x0;
129      }      }
130    else if (p > 0.95)    else if (Q > 0.95)
131      {      {
132        double x0 = exp ((lgamma (a) + log1p (-p)) / a);        double x0 = exp ((gsl_sf_lngamma (a) + log1p (-Q)) / a);
133        x = x0;        x = x0;
134      }      }
135    else    else
136      {      {
137        double xg = gsl_cdf_ugaussian_Qinv (p);        double xg = gsl_cdf_ugaussian_Qinv (Q);
138        double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;        double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;
139        x = x0;        x = x0;
140      }      }
# Line 146  gsl_cdf_gamma_Qinv (double p, double a, Line 146  gsl_cdf_gamma_Qinv (double p, double a,
146     */     */
147    
148    {    {
149      double lambda, dp, phi;      double lambda, dQ, phi;
150    
151    start:    start:
152      dp = -(p - gsl_cdf_gamma_Q (x, a, 1.0));      dQ = Q - gsl_cdf_gamma_Q (x, a, 1.0);
153      phi = gsl_ran_gamma_pdf (x, a, 1.0);      phi = gsl_ran_gamma_pdf (x, a, 1.0);
154    
155      if (dp == 0.0)      if (dQ == 0.0)
156        goto end;        goto end;
157    
158      lambda = dp / GSL_MAX (2 * fabs (dp / x), phi);      lambda = -dQ / GSL_MAX (2 * fabs (dQ / x), phi);
159    
160      {      {
161        double step0 = lambda;        double step0 = lambda;

Legend:
Removed from v.1.1.2.7  
changed lines
  Added in v.1.1.2.8

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