44 |
|
|
45 |
#define MAXI 40 |
#define MAXI 40 |
46 |
|
|
47 |
static double cornish_fisher ( double t, double n ) |
static double |
48 |
|
cornish_fisher (double t, double n) |
49 |
{ |
{ |
50 |
double ret_val = 0.0; |
double ret_val = 0.0; |
51 |
double b; |
double b; |
52 |
double z,zsq; |
double z, zsq; |
53 |
double tmp; |
double tmp; |
54 |
double p[6]; |
double p[6]; |
55 |
int i; |
int i; |
64 |
4960870.65, |
4960870.65, |
65 |
37978595.55, |
37978595.55, |
66 |
201505390.875, |
201505390.875, |
67 |
622437908.625}; |
622437908.625 |
68 |
|
}; |
69 |
const double coeffs5[8] = { |
const double coeffs5[8] = { |
70 |
0.2742857142857142857142, |
0.2742857142857142857142, |
71 |
4.499047619047619047619, |
4.499047619047619047619, |
74 |
12387.6, |
12387.6, |
75 |
101024.55, |
101024.55, |
76 |
559494.0, |
559494.0, |
77 |
1764959.625}; |
1764959.625 |
78 |
const double coeffs4[6] ={ |
}; |
79 |
|
const double coeffs4[6] = { |
80 |
0.3047619047619047619048, |
0.3047619047619047619048, |
81 |
3.752380952380952380952, |
3.752380952380952380952, |
82 |
46.67142857142857142857, |
46.67142857142857142857, |
83 |
427.5, |
427.5, |
84 |
2587.5, |
2587.5, |
85 |
8518.5}; |
8518.5 |
86 |
|
}; |
87 |
const double coeffs3[4] = { |
const double coeffs3[4] = { |
88 |
.4, |
.4, |
89 |
3.3, |
3.3, |
90 |
24.0, |
24.0, |
91 |
85.5}; |
85.5 |
92 |
|
}; |
93 |
|
|
94 |
tmp = n-0.5; |
tmp = n - 0.5; |
95 |
z = tmp*log(1+(t*t/n)); |
z = tmp * log (1 + (t * t / n)); |
96 |
z = sqrt(z); |
z = sqrt (z); |
97 |
b = 48.0 * tmp * tmp; |
b = 48.0 * tmp * tmp; |
98 |
zsq = z*z; |
zsq = z * z; |
99 |
|
|
100 |
|
p[5] = z * poly_eval(coeffs6, 9, zsq); |
101 |
|
p[4] = z * poly_eval(coeffs5, 7, zsq); |
102 |
|
p[3] = z * poly_eval(coeffs4, 5, zsq); |
103 |
|
p[2] = z * poly_eval(coeffs3, 5, zsq); |
104 |
|
|
105 |
p[5] = coeffs6[0] * zsq; |
p[5] = coeffs6[0] * zsq; |
106 |
for ( i = 1; i < 9; i++) |
for (i = 1; i < 9; i++) |
107 |
{ |
{ |
108 |
p[5] = (p[5] + coeffs6[i])*zsq; |
p[5] = (p[5] + coeffs6[i]) * zsq; |
109 |
} |
} |
110 |
p[5] = z * (p[5] + coeffs6[9]); |
p[5] = z * (p[5] + coeffs6[9]); |
111 |
|
|
112 |
p[4] = coeffs5[0] * zsq; |
p[4] = coeffs5[0] * zsq; |
113 |
for(i = 1; i < 7; i++ ) |
for (i = 1; i < 7; i++) |
114 |
{ |
{ |
115 |
p[4] = (p[4] + coeffs5[i]) * zsq; |
p[4] = (p[4] + coeffs5[i]) * zsq; |
116 |
} |
} |
117 |
p[4] = z * ( p[4] + coeffs5[7] ); |
p[4] = z * (p[4] + coeffs5[7]); |
118 |
|
|
119 |
p[3] = coeffs4[0] * zsq; |
p[3] = coeffs4[0] * zsq; |
120 |
for( i = 1; i < 5; i++) |
for (i = 1; i < 5; i++) |
121 |
{ |
{ |
122 |
p[3] = (p[3]+coeffs4[i]) * zsq; |
p[3] = (p[3] + coeffs4[i]) * zsq; |
123 |
} |
} |
124 |
p[3] = z * (p[3] + coeffs4[5]); |
p[3] = z * (p[3] + coeffs4[5]); |
125 |
|
|
126 |
p[2] = coeffs3[0] * zsq; |
p[2] = coeffs3[0] * zsq; |
127 |
for(i = 0; i < 3; i++) |
for (i = 0; i < 3; i++) |
128 |
{ |
{ |
129 |
p[2] = (p[2] + coeffs3[i]) * zsq; |
p[2] = (p[2] + coeffs3[i]) * zsq; |
130 |
} |
} |
131 |
p[2] = z * (p[2] + coeffs3[1]); |
p[2] = z * (p[2] + coeffs3[1]); /* BJG: should this be a [3]? */ |
132 |
|
|
133 |
p[1] = z * (zsq +3.0); |
p[1] = z * (zsq + 3.0); |
134 |
p[0] = z; |
p[0] = z; |
135 |
|
|
136 |
ret_val = p[5]; |
ret_val = p[5]; |
137 |
for( i = 4; i > -1; i-- ) |
for (i = 4; i > -1; i--) |
138 |
{ |
{ |
139 |
ret_val = (ret_val/b)+p[i]; |
ret_val = (ret_val / b) + p[i]; |
140 |
} |
} |
141 |
return ret_val; |
return ret_val; |
142 |
} |
} |
143 |
|
|
144 |
static double tdist_pdf (const double x, const double nu) |
static double |
145 |
|
tdist_pdf (const double x, const double nu) |
146 |
{ |
{ |
147 |
double p; |
double p; |
148 |
|
|
155 |
} |
} |
156 |
|
|
157 |
/* |
/* |
|
* This approximation uses different methods depending on the |
|
|
* degrees of freedom (nu) and the argument t. An alternate method |
|
|
* could use the incomplete beta function. I didn't choose that |
|
|
* method because Kennedy and Gentle state that it is usually inferior to |
|
|
* the method shown here. But "Statistical Computing" was written |
|
|
* a long time ago, and a recent method for computing the |
|
|
* incomplete beta function may be superior. If experience shows |
|
|
* this to be true, the method below can be replaced with that |
|
|
* of the incomplete beta function. |
|
|
*/ |
|
|
/* int gsl_cdf_t_e ( double t, double nu) */ |
|
|
/* { */ |
|
|
/* int rc; */ |
|
|
/* int stopval; */ |
|
|
/* double b; */ |
|
|
/* double ck; */ |
|
|
/* double q; */ |
|
|
/* double y; */ |
|
|
/* double diff; */ |
|
|
/* double ckp2; */ |
|
|
/* gsl_sf_result lg1; */ |
|
|
/* gsl_sf_result lg2; */ |
|
|
/* double num; */ |
|
|
/* int k; */ |
|
|
/* double u; */ |
|
|
/* double p; */ |
|
|
/* int i; */ |
|
|
/* int inu = (int) nu; */ |
|
|
|
|
|
/* /\* */ |
|
|
/* * First approximation is also found in Abramowitz */ |
|
|
/* * and Stegun. This is used only for small values of t */ |
|
|
/* * and small degrees of freedom. */ |
|
|
/* *\/ */ |
|
|
/* y = t / sqrt(nu); */ |
|
|
/* if ( fabs(t) < 4.0 ) */ |
|
|
/* { */ |
|
|
/* diff = fabs(nu-1.0); */ |
|
|
/* if ( diff < GSL_DBL_EPSILON ) */ |
|
|
/* { */ |
|
|
/* q = M_PI_2 * atan(y); */ |
|
|
/* } */ |
|
|
/* else if ( (nu < 21.0) && (nu > 1.0)) */ |
|
|
/* { */ |
|
|
/* ckp2 = 1.0; */ |
|
|
/* b = 1.0 + t*t/nu; */ |
|
|
|
|
|
/* for ( k = inu-2; k > 1; k-=2) */ |
|
|
/* { */ |
|
|
/* ck = 1 + (ckp2 * (k-1)) / ((double)k*b); */ |
|
|
/* ckp2 = ck; */ |
|
|
/* } */ |
|
|
/* if ( (inu%2) == 0 ) */ |
|
|
/* { */ |
|
|
/* q = ck * y / sqrt(b); */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* q = M_PI_2 * (atan(y) + ck * y / b); */ |
|
|
/* } */ |
|
|
/* if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
|
|
/* { */ |
|
|
/* result->val = (1.0 - q)/2; */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* result->val = (1.0+q)/2; */ |
|
|
/* } */ |
|
|
/* } */ |
|
|
/* else if ( nu >= 21.0 ) */ |
|
|
/* { */ |
|
|
/* |
|
|
* Use the Cornish-Fisher expansion to find |
|
|
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
|
|
* Approximate the t cdf with gsl_cdf_gauss. |
|
|
*/ |
|
|
/* u = cornish_fisher ( t, nu ); */ |
|
|
/* q = gsl_cdf_gauss_Q ( u ); */ |
|
|
/* if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
|
|
/* { */ |
|
|
/* result->val = q; */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* result->val = 1.0 - q; */ |
|
|
/* } */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* GSL_ERROR_VAL ("degrees of freedom are less than 1", */ |
|
|
/* GSL_EINVAL, nu ); */ |
|
|
/* } */ |
|
|
/* } */ |
|
|
/* return rc; */ |
|
|
/* } */ |
|
|
|
|
|
/* |
|
158 |
* First approximation is also found in Abramowitz |
* First approximation is also found in Abramowitz |
159 |
* and Stegun. This is used only for small values of t |
* and Stegun. This is used only for small values of t |
160 |
* and small degrees of freedom. |
* and small degrees of freedom. |
161 |
*/ |
*/ |
162 |
static double gsl_cdf_t_smallx ( const double t, const double nu ) |
static double |
163 |
|
gsl_cdf_t_smallx (const double t, const double nu) |
164 |
{ |
{ |
165 |
double y; |
double y; |
166 |
double q; |
double q; |
172 |
int k; |
int k; |
173 |
int inu = (int) nu; |
int inu = (int) nu; |
174 |
|
|
175 |
y = t / sqrt(nu); |
y = t / sqrt (nu); |
176 |
diff = fabs(nu-1.0); |
diff = fabs (nu - 1.0); |
177 |
if ( diff < GSL_DBL_EPSILON ) |
if (diff < GSL_DBL_EPSILON) |
178 |
{ |
{ |
179 |
q = M_PI_2 * atan(y); |
q = M_PI_2 * atan (y); |
180 |
} |
} |
181 |
else if ( (nu < 21.0) && (nu > 1.0)) |
else if ((nu < 21.0) && (nu > 1.0)) |
182 |
{ |
{ |
183 |
ckp2 = 1.0; |
ckp2 = 1.0; |
184 |
b = 1.0 + t*t/nu; |
b = 1.0 + t * t / nu; |
185 |
|
|
186 |
for ( k = inu-2; k > 1; k-=2) |
for (k = inu - 2; k > 1; k -= 2) |
187 |
{ |
{ |
188 |
ck = 1.0 + (ckp2 * (k-1)) / ((double)k*b); |
ck = 1.0 + (ckp2 * (k - 1)) / ((double) k * b); |
189 |
ckp2 = ck; |
ckp2 = ck; |
190 |
} |
} |
191 |
if ( (inu%2) == 0 ) |
if ((inu % 2) == 0) |
192 |
{ |
{ |
193 |
q = ck * y / sqrt(b); |
q = ck * y / sqrt (b); |
194 |
} |
} |
195 |
else |
else |
196 |
{ |
{ |
197 |
q = M_PI_2 * (atan(y) + ck * y / b); |
q = M_PI_2 * (atan (y) + ck * y / b); |
198 |
} |
} |
|
/* if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
|
|
/* { */ |
|
|
/* result->val = (1.0 - q)/2; */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* result->val = (1.0+q)/2; */ |
|
|
/* } */ |
|
199 |
} |
} |
200 |
return q; |
return q; |
201 |
} |
} |
202 |
|
|
203 |
/* |
/* |
204 |
* Use the Cornish-Fisher expansion to find |
* Use the Cornish-Fisher expansion to find |
205 |
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
206 |
* Approximate the t cdf with gsl_cdf_gauss. |
* Approximate the t cdf with gsl_cdf_gauss. |
207 |
*/ |
*/ |
208 |
static double t_cornish_fisher ( const double x, const double nu ) |
static double |
209 |
|
t_cornish_fisher (const double x, const double nu) |
210 |
{ |
{ |
211 |
double y; |
double y; |
212 |
double q; |
double q; |
213 |
double u; |
double u; |
214 |
|
|
215 |
y = x / sqrt(nu); |
y = x / sqrt (nu); |
216 |
u = cornish_fisher ( x, nu ); |
u = cornish_fisher (x, nu); |
217 |
q = gsl_cdf_ugaussian_Q ( u); |
q = gsl_cdf_ugaussian_Q (u); |
|
/* if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
|
|
/* { */ |
|
|
/* result->val = q; */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* result->val = 1.0 - q; */ |
|
|
/* } */ |
|
218 |
return q; |
return q; |
219 |
} |
} |
220 |
|
|
221 |
/* |
/* |
222 |
* Series approximation for t > 4.0. This needs to be fixed; |
* Series approximation for t > 4.0. This needs to be fixed; |
223 |
* it shouldn't subtract the result from 1.0. A better way is |
* it shouldn't subtract the result from 1.0. A better way is |
225 |
* means rummaging through Fisher's paper in Metron, v5, 1926, |
* means rummaging through Fisher's paper in Metron, v5, 1926, |
226 |
* "Expansion of Student's integral in powers of n^{-1}." |
* "Expansion of Student's integral in powers of n^{-1}." |
227 |
*/ |
*/ |
228 |
static double normal_approx ( const double x, const double nu ) |
static double |
229 |
|
normal_approx (const double x, const double nu) |
230 |
{ |
{ |
231 |
double y; |
double y; |
232 |
double num; |
double num; |
238 |
gsl_sf_result lg1; |
gsl_sf_result lg1; |
239 |
gsl_sf_result lg2; |
gsl_sf_result lg2; |
240 |
|
|
241 |
y = 1/sqrt(1+x*x/nu); |
y = 1 / sqrt (1 + x * x / nu); |
242 |
num = 1.0; |
num = 1.0; |
243 |
q = 0.0; |
q = 0.0; |
244 |
diff = 2*GSL_DBL_EPSILON; |
diff = 2 * GSL_DBL_EPSILON; |
245 |
for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2) |
for (i = 2; (i < MAXI) && (diff > GSL_DBL_EPSILON); i += 2) |
246 |
{ |
{ |
247 |
diff = q; |
diff = q; |
248 |
num *= y*y*(i-1)/i; |
num *= y * y * (i - 1) / i; |
249 |
q += num / (nu+i); |
q += num / (nu + i); |
250 |
diff = q - diff; |
diff = q - diff; |
251 |
} |
} |
252 |
q += 1/nu; |
q += 1 / nu; |
253 |
rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1); |
rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1); |
254 |
if ( rc == GSL_SUCCESS ) |
if (rc == GSL_SUCCESS) |
255 |
{ |
{ |
256 |
rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2); |
rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2); |
257 |
if ( rc != GSL_SUCCESS ) |
if (rc != GSL_SUCCESS) |
258 |
{ |
{ |
259 |
return rc; |
return rc; |
260 |
} |
} |
261 |
} |
} |
262 |
diff = (lg2.val) - (lg1.val); |
diff = (lg2.val) - (lg1.val); |
263 |
q *= pow(y,nu) * exp (diff) / sqrt (M_PI); |
q *= pow (y, nu) * exp (diff) / sqrt (M_PI); |
264 |
|
|
|
/* if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
|
|
/* { */ |
|
|
/* result->val = q; */ |
|
|
/* } */ |
|
|
/* else if ( tail == GSL_CDF_LOWER ) */ |
|
|
/* { */ |
|
|
/* result->val = 1.0 - q; */ |
|
|
/* } */ |
|
|
|
|
265 |
return q; |
return q; |
266 |
} |
} |
267 |
|
|
268 |
double gsl_cdf_t_P ( const double x, const double nu ) |
double |
269 |
|
gsl_cdf_t_P (const double x, const double nu) |
270 |
{ |
{ |
271 |
double val = 0.0; |
double val = 0.0; |
272 |
double q; |
double q; |
273 |
double absx = fabs(x); |
double absx = fabs (x); |
274 |
|
|
275 |
if ( fabs(x) < 4.0 ) |
if (fabs (x) < 4.0) |
276 |
{ |
{ |
277 |
q = gsl_cdf_t_smallx (x, nu ); |
q = gsl_cdf_t_smallx (x, nu); |
278 |
if ( x <= 0.0 ) |
if (x <= 0.0) |
279 |
{ |
{ |
280 |
val = (1.0 - q) / 2.0 ; |
val = (1.0 - q) / 2.0; |
281 |
} |
} |
282 |
else |
else |
283 |
{ |
{ |
284 |
val = (1.0 + q ) / 2.0; |
val = (1.0 + q) / 2.0; |
285 |
} |
} |
286 |
} |
} |
287 |
else if ( nu >= 21.0 ) |
else if (nu >= 21.0) |
288 |
{ |
{ |
289 |
q = t_cornish_fisher ( x, nu ); |
q = t_cornish_fisher (x, nu); |
290 |
if ( x <= 0.0) |
if (x <= 0.0) |
291 |
{ |
{ |
292 |
val = q; |
val = q; |
293 |
} |
} |
298 |
} |
} |
299 |
else |
else |
300 |
{ |
{ |
301 |
val = normal_approx ( x, nu ); |
val = normal_approx (x, nu); |
302 |
if ( x > 0.0 ) |
if (x > 0.0) |
303 |
{ |
{ |
304 |
val = 1.0 - val; |
val = 1.0 - val; |
305 |
} |
} |
306 |
} |
} |
307 |
return val; |
return val; |
308 |
} |
} |
309 |
double gsl_cdf_t_Q ( const double x, const double nu ) |
|
310 |
|
double |
311 |
|
gsl_cdf_t_Q (const double x, const double nu) |
312 |
{ |
{ |
313 |
double val = 0.0; |
double val = 0.0; |
314 |
double q; |
double q; |
315 |
double absx = fabs(x); |
double absx = fabs (x); |
316 |
|
|
317 |
if ( fabs(x) < 4.0 ) |
if (fabs (x) < 4.0) |
318 |
{ |
{ |
319 |
q = gsl_cdf_t_smallx ( x, nu ); |
q = gsl_cdf_t_smallx (x, nu); |
320 |
if ( x > 0.0 ) |
if (x > 0.0) |
321 |
{ |
{ |
322 |
val = (1.0 - q) / 2.0 ; |
val = (1.0 - q) / 2.0; |
323 |
} |
} |
324 |
else |
else |
325 |
{ |
{ |
326 |
val = (1.0 + q ) / 2.0; |
val = (1.0 + q) / 2.0; |
327 |
} |
} |
328 |
} |
} |
329 |
else if ( nu >= 21.0 ) |
else if (nu >= 21.0) |
330 |
{ |
{ |
331 |
val = t_cornish_fisher ( x, nu ); |
val = t_cornish_fisher (x, nu); |
332 |
if ( x <= 0.0) |
if (x <= 0.0) |
333 |
{ |
{ |
334 |
val = 1.0 - val; |
val = 1.0 - val; |
335 |
} |
} |
336 |
} |
} |
337 |
else |
else |
338 |
{ |
{ |
339 |
val = normal_approx ( x, nu ); |
val = normal_approx (x, nu); |
340 |
if ( x <= 0.0 ) |
if (x <= 0.0) |
341 |
{ |
{ |
342 |
val = 1.0 - val; |
val = 1.0 - val; |
343 |
} |
} |
344 |
} |
} |
345 |
return val; |
return val; |
346 |
} |
} |
347 |
|
|
348 |
/* |
/* |
349 |
* Invert the T distribution. Uses a method as shown in |
* Invert the T distribution. Uses a method as shown in |
350 |
* Statistical Computing, 5.4.2. This method uses an initial |
* Statistical Computing, 5.4.2. This method uses an initial |
352 |
* the initial method with a Taylor series or a Cornish-Fisher |
* the initial method with a Taylor series or a Cornish-Fisher |
353 |
* expansion. |
* expansion. |
354 |
*/ |
*/ |
355 |
static double inv_cornish_fisher ( double z, double n ) |
static double |
356 |
|
inv_cornish_fisher (double z, double n) |
357 |
{ |
{ |
358 |
double ret_val = 0.0; |
double ret_val = 0.0; |
359 |
double b; |
double b; |
368 |
double zseven; |
double zseven; |
369 |
int i; |
int i; |
370 |
|
|
371 |
tmp = n-0.5; |
tmp = n - 0.5; |
372 |
b = 48.0 * tmp * tmp; |
b = 48.0 * tmp * tmp; |
373 |
if ( n > 5.0 ) |
if (n > 5.0) |
374 |
{ |
{ |
375 |
c = 96.36 - 16.0 / tmp - 98.0 / (tmp*tmp) + 20700.0 / (tmp*tmp*tmp*b); |
c = |
376 |
|
96.36 - 16.0 / tmp - 98.0 / (tmp * tmp) + |
377 |
|
20700.0 / (tmp * tmp * tmp * b); |
378 |
} |
} |
379 |
else |
else |
380 |
{ |
{ |
381 |
c = 0.3 * (n - 4.5) * (z + 0.6); |
c = 0.3 * (n - 4.5) * (z + 0.6); |
382 |
} |
} |
383 |
zsq = z*z; |
zsq = z * z; |
384 |
zcube = zsq*z; |
zcube = zsq * z; |
385 |
zfour = zcube*z; |
zfour = zcube * z; |
386 |
zfive = zfour*z; |
zfive = zfour * z; |
387 |
zseven = zfour*zcube; |
zseven = zfour * zcube; |
388 |
d = n * M_SQRTPI * gsl_sf_gamma ( tmp ) |
d = n * M_SQRTPI * gsl_sf_gamma (tmp) / (2.0 * gsl_sf_gamma (tmp + 0.5)); |
389 |
/ (2.0*gsl_sf_gamma ( tmp + 0.5 )); |
u = |
390 |
u = 10.0 * b * ( b + c - 2.0*z - 7.0 * zsq - 5.0*zcube + 0.05*d*zfour); |
10.0 * b * (b + c - 2.0 * z - 7.0 * zsq - 5.0 * zcube + 0.05 * d * zfour); |
391 |
tmp = 4.0*zseven + 63.0*zfive + 360.0*zcube + 945.0*z; |
tmp = 4.0 * zseven + 63.0 * zfive + 360.0 * zcube + 945.0 * z; |
392 |
ret_val = z - (zcube+3.0*z)/b + tmp/u; |
ret_val = z - (zcube + 3.0 * z) / b + tmp / u; |
393 |
return ret_val; |
return ret_val; |
394 |
} |
} |
395 |
|
|
396 |
double gsl_cdf_ut_P_inv (double prob, double nu) |
double |
397 |
|
gsl_cdf_ut_P_inv (double prob, double nu) |
398 |
{ |
{ |
399 |
double initial_result; |
double initial_result; |
400 |
double result; |
double result; |
419 |
double c3; |
double c3; |
420 |
double c4; |
double c4; |
421 |
|
|
422 |
if ( prob < 0.0 ) |
if (prob < 0.0) |
423 |
{ |
{ |
424 |
return GSL_EDOM; |
return GSL_EDOM; |
425 |
} |
} |
426 |
if ( prob > 1.0 ) |
if (prob > 1.0) |
427 |
{ |
{ |
428 |
return GSL_EDOM; |
return GSL_EDOM; |
429 |
} |
} |
430 |
if ( nu < 0.0 ) |
if (nu < 0.0) |
431 |
{ |
{ |
432 |
return GSL_EDOM; |
return GSL_EDOM; |
433 |
} |
} |
434 |
if ( fabs(prob) < GSL_DBL_EPSILON ) |
if (fabs (prob) < GSL_DBL_EPSILON) |
435 |
{ |
{ |
436 |
return GSL_POSINF; |
return GSL_POSINF; |
437 |
} |
} |
438 |
if ( fabs ( 1.0 - prob ) < GSL_DBL_EPSILON ) |
if (fabs (1.0 - prob) < GSL_DBL_EPSILON) |
439 |
{ |
{ |
440 |
return GSL_NEGINF; |
return GSL_NEGINF; |
441 |
} |
} |
442 |
printf("prob is %f\t",prob); |
printf ("prob is %f\t", prob); |
443 |
tmp = nu / 2.0; |
tmp = nu / 2.0; |
444 |
d = tmp * M_SQRTPI * gsl_sf_gamma ( tmp ) |
d = tmp * M_SQRTPI * gsl_sf_gamma (tmp) / gsl_sf_gamma (tmp + 0.5); |
445 |
/ gsl_sf_gamma ( tmp + 0.5 ); |
method_test = gsl_max (d * prob, 2.0 / nu); |
|
method_test = gsl_max ( d*prob, 2.0 / nu); |
|
446 |
/* |
/* |
447 |
* There are two possible initial approximations. |
* There are two possible initial approximations. |
448 |
* Which is used depends on prob and nu. |
* Which is used depends on prob and nu. |
449 |
*/ |
*/ |
450 |
if ( method_test > 0.05 ) |
if (method_test > 0.05) |
451 |
{ |
{ |
452 |
tmp = prob; |
tmp = prob; |
453 |
a = nu - 0.5; |
a = nu - 0.5; |
454 |
/* gsl_cdf_ugaussian_P_inv(tmp) ? */ |
/* gsl_cdf_ugaussian_P_inv(tmp) ? */ |
455 |
x = gsl_cdf_ugaussian_P_inv (tmp ); |
x = gsl_cdf_ugaussian_P_inv (tmp); |
456 |
y = inv_cornish_fisher ( x, nu ); |
y = inv_cornish_fisher (x, nu); |
457 |
tsqn = -1.0 + exp ( a*y*y); |
tsqn = -1.0 + exp (a * y * y); |
458 |
} |
} |
459 |
else |
else |
460 |
{ |
{ |
461 |
z = pow ( prob * d, 1 / tmp ); |
z = pow (prob * d, 1 / tmp); |
462 |
zz = z * z; |
zz = z * z; |
463 |
zzz = z*z*z; |
zzz = z * z * z; |
464 |
tsqn = 1/z + (nu + 1.0)*(-1.0 + z/(2.0 * (nu + 4.0)) + |
tsqn = 1 / z + (nu + 1.0) * (-1.0 + z / (2.0 * (nu + 4.0)) + |
465 |
nu * zz / (3.0 * (nu+2.0)*(nu+6.0)) |
nu * zz / (3.0 * (nu + 2.0) * (nu + 6.0)) |
466 |
+ nu*(nu+3.0)*(2.0*nu*nu+9.0*nu-2.0)*zzz |
+ nu * (nu + 3.0) * (2.0 * nu * nu + |
467 |
/(8.0 * (nu + 2.0) * (nu + 2.0) |
9.0 * nu - |
468 |
* (nu + 4.0) * (nu + 4.0) * |
2.0) * zzz / (8.0 * |
469 |
(nu + 8.0))) |
(nu + |
470 |
/(nu + 2); |
2.0) * |
471 |
|
(nu + |
472 |
} |
2.0) * |
473 |
initial_result = sqrt(tsqn * nu); |
(nu + |
474 |
printf("initial_result = %f \t",initial_result); |
4.0) * |
475 |
tcdf = gsl_cdf_t_Q( initial_result, nu ); |
(nu + |
476 |
tmp = 0.5 * ( tcdf - prob ); |
4.0) * |
477 |
tmp2 = tdist_pdf ( initial_result, nu ); |
(nu + |
478 |
|
8.0))) |
479 |
|
/ (nu + 2); |
480 |
|
|
481 |
|
} |
482 |
|
initial_result = sqrt (tsqn * nu); |
483 |
|
printf ("initial_result = %f \t", initial_result); |
484 |
|
tcdf = gsl_cdf_t_Q (initial_result, nu); |
485 |
|
tmp = 0.5 * (tcdf - prob); |
486 |
|
tmp2 = tdist_pdf (initial_result, nu); |
487 |
w = tmp / tmp2; |
w = tmp / tmp2; |
488 |
psi = initial_result * (nu+1.0) / |
psi = initial_result * (nu + 1.0) / (nu + initial_result * initial_result); |
489 |
(nu + initial_result * initial_result); |
psi_prime = (nu + 1.0) * (nu - initial_result * initial_result) / |
490 |
psi_prime = (nu+1.0) * (nu - initial_result * initial_result) / |
((nu + initial_result * initial_result) * |
|
((nu + initial_result*initial_result) * |
|
491 |
(nu + initial_result * initial_result)); |
(nu + initial_result * initial_result)); |
492 |
d_psi_dt = 2 * initial_result * nu * (nu + 1.0) / |
d_psi_dt = 2 * initial_result * nu * (nu + 1.0) / |
493 |
((nu + initial_result * initial_result) * |
((nu + initial_result * initial_result) * |
497 |
((nu + initial_result * initial_result) * |
((nu + initial_result * initial_result) * |
498 |
(nu + initial_result * initial_result) * |
(nu + initial_result * initial_result) * |
499 |
(nu + initial_result * initial_result)); |
(nu + initial_result * initial_result)); |
500 |
c2 = psi/2.0; |
c2 = psi / 2.0; |
501 |
c3 = (2.0 * psi * psi + psi_prime) / 6.0; |
c3 = (2.0 * psi * psi + psi_prime) / 6.0; |
502 |
c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) + |
c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) + |
503 |
4.0 * psi * d_psi_dt + d_psiprime_dt; |
4.0 * psi * d_psi_dt + d_psiprime_dt; |
504 |
c4 /= 24.0; |
c4 /= 24.0; |
505 |
result = initial_result + w + c2 * w * w + |
result = initial_result + w + c2 * w * w + |
506 |
c3 * w * w * w + c4 * w * w * w * w; |
c3 * w * w * w + c4 * w * w * w * w; |
507 |
printf("result = %f\n",result); |
printf ("result = %f\n", result); |
508 |
return result; |
return result; |
509 |
} |
} |