/[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.5 by bjg, Sat Jul 12 22:00:18 2003 UTC revision 1.1.2.6 by bjg, Mon Jul 14 12:07:33 2003 UTC
# Line 46  gsl_cdf_gamma_Pinv (double p, double a, Line 46  gsl_cdf_gamma_Pinv (double p, double a,
46    
47    if (p < 0.05)    if (p < 0.05)
48      {      {
49        double x0 = exp ((lgamma (b) + log (p)) / b);        double x0 = exp ((lgamma (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 = -log (1 - p) + lgamma (b);        double x0 = -log (1 - p) + lgamma (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 (b)) ? b : sqrt (b) * xg + b;        double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;
61        x = x0;        x = x0;
62      }      }
63    
# Line 71  gsl_cdf_gamma_Pinv (double p, double a, Line 71  gsl_cdf_gamma_Pinv (double p, double a,
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, a, 1.0);
75      phi = gsl_ran_gamma_pdf (x, b, 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;
# Line 81  gsl_cdf_gamma_Pinv (double p, double a, Line 81  gsl_cdf_gamma_Pinv (double p, double a,
81    
82      {      {
83        double step0 = lambda;        double step0 = lambda;
84        double step1 = -((b - 1) / x - 1) * lambda * lambda / 4.0;        double step1 = -((a - 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))
# Line 101  gsl_cdf_gamma_Pinv (double p, double a, Line 101  gsl_cdf_gamma_Pinv (double p, double a,
101    }    }
102    
103  end:  end:
104    return a * x;    return b * x;
105  }  }
106    
107  double  double
# Line 124  gsl_cdf_gamma_Qinv (double p, double a, Line 124  gsl_cdf_gamma_Qinv (double p, double a,
124    
125    if (p < 0.05)    if (p < 0.05)
126      {      {
127        double x0 = -log (p) + lgamma (b);        double x0 = -log (p) + lgamma (a);
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 (a) + log1p (-p)) / 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 (p);
138        double x0 = (xg < -sqrt (b)) ? b : sqrt (b) * xg + b;        double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;
139        x = x0;        x = x0;
140      }      }
141    
# Line 149  gsl_cdf_gamma_Qinv (double p, double a, Line 149  gsl_cdf_gamma_Qinv (double p, double a,
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, a, 1.0));
153      phi = gsl_ran_gamma_pdf (x, 1.0, b);      phi = gsl_ran_gamma_pdf (x, a, 1.0);
154    
155      if (dp == 0.0)      if (dp == 0.0)
156        goto end;        goto end;
# Line 159  gsl_cdf_gamma_Qinv (double p, double a, Line 159  gsl_cdf_gamma_Qinv (double p, double a,
159    
160      {      {
161        double step0 = lambda;        double step0 = lambda;
162        double step1 = -((b - 1) / x - 1) * lambda * lambda / 4.0;        double step1 = -((a - 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))
# Line 179  gsl_cdf_gamma_Qinv (double p, double a, Line 179  gsl_cdf_gamma_Qinv (double p, double a,
179    }    }
180    
181  end:  end:
182    return a * x;    return b * x;
183  }  }

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

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