46 |
double |
double |
47 |
gsl_cdf_tdist_Pinv (const double p, const double nu) |
gsl_cdf_tdist_Pinv (const double p, const double nu) |
48 |
{ |
{ |
49 |
double x; |
double x, ptail; |
50 |
|
|
51 |
if (p == 1.0) |
if (p == 1.0) |
52 |
{ |
{ |
57 |
return GSL_NEGINF; |
return GSL_NEGINF; |
58 |
} |
} |
59 |
|
|
60 |
|
if (nu == 1.0) |
|
if (p < 0.05) |
|
61 |
{ |
{ |
62 |
double beta = gsl_sf_beta (0.5, nu / 2); |
x = tan (M_PI * (p - 0.5)); |
|
x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
|
63 |
} |
} |
64 |
else if (p > 0.95) |
else if (nu == 2.0) |
65 |
{ |
{ |
66 |
double beta = gsl_sf_beta (0.5, nu / 2); |
double a = 2 * p - 1; |
67 |
x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
x = a / sqrt (2 * (1 - a * a)); |
68 |
} |
} |
69 |
else |
|
70 |
|
ptail = (p < 0.5) ? p : 1 - p; |
71 |
|
|
72 |
|
if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2)) |
73 |
{ |
{ |
74 |
double xg = gsl_cdf_ugaussian_Pinv (p); |
double xg = gsl_cdf_ugaussian_Pinv (p); |
75 |
x = inv_cornish_fisher (xg, nu); |
x = inv_cornish_fisher (xg, nu); |
76 |
} |
} |
77 |
|
else |
78 |
|
{ |
79 |
|
/* Use an asymptotic expansion of the tail of integral */ |
80 |
|
|
81 |
|
double beta = gsl_sf_beta (0.5, nu / 2); |
82 |
|
|
83 |
|
if (p < 0.5) |
84 |
|
{ |
85 |
|
x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
86 |
|
} |
87 |
|
else |
88 |
|
{ |
89 |
|
x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
90 |
|
} |
91 |
|
|
92 |
|
/* Correct nu -> nu/(1+nu/x^2) in the leading term to account |
93 |
|
for higher order terms. This avoids overestimating x, which |
94 |
|
makes the iteration unstable due to the rapidly decreasing |
95 |
|
tails of the distribution. */ |
96 |
|
|
97 |
|
x /= sqrt (1 + nu / (x * x)); |
98 |
|
} |
99 |
|
|
100 |
{ |
{ |
101 |
double dp, phi; |
double dp, phi; |
115 |
double step = step0; |
double step = step0; |
116 |
|
|
117 |
if (fabs (step1) < fabs (step0)) |
if (fabs (step1) < fabs (step0)) |
118 |
step += step1; |
{ |
119 |
|
step += step1; |
120 |
printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu)); |
} |
121 |
|
|
122 |
x += step; |
if (p > 0.5 && x + step < 0) |
123 |
|
x /= 2; |
124 |
|
else if (p < 0.5 && x + step > 0) |
125 |
|
x /= 2; |
126 |
|
else |
127 |
|
x += step; |
128 |
|
|
129 |
if (fabs (step) > 1e-10 * fabs (x)) |
if (fabs (step) > 1e-10 * fabs (x)) |
130 |
goto start; |
goto start; |
139 |
double |
double |
140 |
gsl_cdf_tdist_Qinv (const double p, const double nu) |
gsl_cdf_tdist_Qinv (const double p, const double nu) |
141 |
{ |
{ |
142 |
double x; |
double x, ptail; |
143 |
|
|
144 |
if (p == 0.0) |
if (p == 0.0) |
145 |
{ |
{ |
150 |
return GSL_NEGINF; |
return GSL_NEGINF; |
151 |
} |
} |
152 |
|
|
153 |
{ |
if (nu == 1.0) |
154 |
double xg = gsl_cdf_ugaussian_Qinv (p); |
{ |
155 |
x = inv_cornish_fisher (xg, nu); |
x = tan (M_PI * (0.5 - p)); |
156 |
} |
} |
157 |
|
else if (nu == 2.0) |
158 |
|
{ |
159 |
|
double a = 2 * (1 - p) - 1; |
160 |
|
x = a / sqrt (2 * (1 - a * a)); |
161 |
|
} |
162 |
|
|
163 |
|
ptail = (p < 0.5) ? p : 1 - p; |
164 |
|
|
165 |
if (x > nu) |
if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2)) |
166 |
{ |
{ |
167 |
double beta = gsl_sf_beta (0.5, nu / 2); |
double xg = gsl_cdf_ugaussian_Qinv (p); |
168 |
x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
x = inv_cornish_fisher (xg, nu); |
169 |
} |
} |
170 |
else if (x < -nu) |
else |
171 |
{ |
{ |
172 |
|
/* Use an asymptotic expansion of the tail of integral */ |
173 |
|
|
174 |
double beta = gsl_sf_beta (0.5, nu / 2); |
double beta = gsl_sf_beta (0.5, nu / 2); |
175 |
x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
|
176 |
|
if (p < 0.5) |
177 |
|
{ |
178 |
|
x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
179 |
|
} |
180 |
|
else |
181 |
|
{ |
182 |
|
x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
183 |
|
} |
184 |
|
|
185 |
|
/* Correct nu -> nu/(1+nu/x^2) in the leading term to account |
186 |
|
for higher order terms. This avoids overestimating x, which |
187 |
|
makes the iteration unstable due to the rapidly decreasing |
188 |
|
tails of the distribution. */ |
189 |
|
|
190 |
|
x /= sqrt (1 + nu / (x * x)); |
191 |
} |
} |
192 |
|
|
193 |
{ |
{ |
208 |
double step = step0; |
double step = step0; |
209 |
|
|
210 |
if (fabs (step1) < fabs (step0)) |
if (fabs (step1) < fabs (step0)) |
211 |
step += step1; |
{ |
212 |
|
step += step1; |
213 |
printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu)); |
} |
214 |
|
|
215 |
x += step; |
if (p < 0.5 && x + step < 0) |
216 |
|
x /= 2; |
217 |
|
else if (p > 0.5 && x + step > 0) |
218 |
|
x /= 2; |
219 |
|
else |
220 |
|
x += step; |
221 |
|
|
222 |
if (fabs (step) > 1e-10 * fabs (x)) |
if (fabs (step) > 1e-10 * fabs (x)) |
223 |
goto start; |
goto start; |