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 |
} |
} |
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 |
} |
} |
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; |
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 |
} |
} |
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 |
} |
} |
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; |