26 |
|
|
27 |
#include <stdio.h> |
#include <stdio.h> |
28 |
|
|
29 |
|
static double |
30 |
|
inv_cornish_fisher (double z, double nu) |
31 |
|
{ |
32 |
|
double a = 1 / (nu - 0.5); |
33 |
|
double b = 48.0 / (a * a); |
34 |
|
|
35 |
|
double cf1 = z * (3 + z * z); |
36 |
|
double cf2 = z * (945 + z * z * (360 + z * z * (63 + z * z * 4))); |
37 |
|
|
38 |
|
double y = z - cf1 / b + cf2 / (10 * b * b); |
39 |
|
|
40 |
|
double t = GSL_SIGN (z) * sqrt (nu * expm1 (a * y * y)); |
41 |
|
|
42 |
|
return t; |
43 |
|
} |
44 |
|
|
45 |
|
|
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; |
50 |
|
|
51 |
if (p == 1.0) |
if (p == 1.0) |
52 |
{ |
{ |
53 |
return GSL_POSINF; |
return GSL_POSINF; |
57 |
return GSL_NEGINF; |
return GSL_NEGINF; |
58 |
} |
} |
59 |
|
|
60 |
if (p < 0.05) |
|
61 |
|
if (p < 0.05) |
62 |
{ |
{ |
63 |
|
double beta = gsl_sf_beta (0.5, nu / 2); |
64 |
|
x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
65 |
} |
} |
66 |
else if (p > 0.95) |
else if (p > 0.95) |
67 |
{ |
{ |
68 |
|
double beta = gsl_sf_beta (0.5, nu / 2); |
69 |
|
x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
70 |
} |
} |
71 |
else |
else |
72 |
{ |
{ |
73 |
double xg = gsl_cdf_ugaussian_Pinv (p); |
double xg = gsl_cdf_ugaussian_Pinv (p); |
74 |
double x0 = xg + xg*(xg*xg+1.0)/(4.0*nu) ; |
x = inv_cornish_fisher (xg, nu); |
|
x = x0; |
|
75 |
} |
} |
76 |
|
|
77 |
|
|
78 |
{ |
{ |
79 |
double dp, phi; |
double dp, phi; |
80 |
|
|
86 |
goto end; |
goto end; |
87 |
|
|
88 |
{ |
{ |
89 |
double step = dp / phi; |
double lambda = dp / phi; |
90 |
|
double step0 = lambda; |
91 |
|
double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0); |
92 |
|
|
93 |
|
double step = step0; |
94 |
|
|
95 |
|
if (fabs (step1) < fabs (step0)) |
96 |
|
step += step1; |
97 |
|
|
98 |
|
printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu)); |
99 |
|
|
|
printf("x=%.18e step=%.18e\n", x, step); |
|
100 |
x += step; |
x += step; |
101 |
|
|
102 |
if (fabs (step) > 1e-10 * fabs(x)) |
if (fabs (step) > 1e-10 * fabs (x)) |
103 |
goto start; |
goto start; |
104 |
} |
} |
105 |
} |
} |
106 |
|
|
109 |
return x; |
return x; |
110 |
} |
} |
111 |
|
|
112 |
#if 0 |
double |
113 |
/* |
gsl_cdf_tdist_Qinv (const double p, const double nu) |
|
* Invert the T distribution. Uses a method as shown in |
|
|
* Statistical Computing, 5.4.2. This method uses an initial |
|
|
* approximation based on Hill's method above, then improves |
|
|
* the initial method with a Taylor series or a Cornish-Fisher |
|
|
* expansion. |
|
|
*/ |
|
|
static double |
|
|
inv_cornish_fisher (double z, double n) |
|
|
{ |
|
|
double ret_val = 0.0; |
|
|
double b; |
|
|
double zsq; |
|
|
double tmp; |
|
|
double c; |
|
|
double u; |
|
|
double d; |
|
|
double zcube; |
|
|
double zfour; |
|
|
double zfive; |
|
|
double zseven; |
|
|
int i; |
|
|
|
|
|
tmp = n - 0.5; |
|
|
b = 48.0 * tmp * tmp; |
|
|
if (n > 5.0) |
|
|
{ |
|
|
c = |
|
|
96.36 - 16.0 / tmp - 98.0 / (tmp * tmp) + |
|
|
20700.0 / (tmp * tmp * tmp * b); |
|
|
} |
|
|
else |
|
|
{ |
|
|
c = 0.3 * (n - 4.5) * (z + 0.6); |
|
|
} |
|
|
zsq = z * z; |
|
|
zcube = zsq * z; |
|
|
zfour = zcube * z; |
|
|
zfive = zfour * z; |
|
|
zseven = zfour * zcube; |
|
|
d = n * M_SQRTPI * gsl_sf_gamma (tmp) / (2.0 * gsl_sf_gamma (tmp + 0.5)); |
|
|
u = |
|
|
10.0 * b * (b + c - 2.0 * z - 7.0 * zsq - 5.0 * zcube + 0.05 * d * zfour); |
|
|
tmp = 4.0 * zseven + 63.0 * zfive + 360.0 * zcube + 945.0 * z; |
|
|
ret_val = z - (zcube + 3.0 * z) / b + tmp / u; |
|
|
return ret_val; |
|
|
} |
|
|
|
|
|
static double |
|
|
tdist_pdf (const double x, const double nu) |
|
114 |
{ |
{ |
|
double p; |
|
|
|
|
|
double lg1 = gsl_sf_lngamma (nu / 2); |
|
|
double lg2 = gsl_sf_lngamma ((nu + 1) / 2); |
|
|
|
|
|
p = ((exp (lg2 - lg1) / sqrt (M_PI * nu)) |
|
|
* pow ((1 + x * x / nu), -(nu + 1) / 2)); |
|
|
return p; |
|
|
} |
|
|
|
|
|
double initial_result; |
|
|
double result; |
|
|
double d; |
|
|
double tmp; |
|
|
double tmp2; |
|
|
double z; |
|
|
double zz; |
|
|
double zzz; |
|
|
double y; |
|
|
double tsqn; |
|
|
double a; |
|
|
double method_test; |
|
|
double psi; |
|
|
double psi_prime; |
|
|
double d_psi_dt; |
|
|
double d_psiprime_dt; |
|
|
double tcdf; |
|
115 |
double x; |
double x; |
|
double w; |
|
|
double c2; |
|
|
double c3; |
|
|
double c4; |
|
116 |
|
|
117 |
if (prob < 0.0) |
if (p == 0.0) |
118 |
{ |
{ |
119 |
return GSL_EDOM; |
return GSL_POSINF; |
120 |
} |
} |
121 |
if (prob > 1.0) |
else if (p == 1.0) |
122 |
{ |
{ |
123 |
return GSL_EDOM; |
return GSL_NEGINF; |
124 |
} |
} |
125 |
if (nu < 0.0) |
|
126 |
|
{ |
127 |
|
double xg = gsl_cdf_ugaussian_Qinv (p); |
128 |
|
x = inv_cornish_fisher (xg, nu); |
129 |
|
} |
130 |
|
|
131 |
|
if (x > nu) |
132 |
{ |
{ |
133 |
return GSL_EDOM; |
double beta = gsl_sf_beta (0.5, nu / 2); |
134 |
|
x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
135 |
} |
} |
136 |
if (fabs (prob) < GSL_DBL_EPSILON) |
else if (x < -nu) |
137 |
{ |
{ |
138 |
return GSL_POSINF; |
double beta = gsl_sf_beta (0.5, nu / 2); |
139 |
|
x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
140 |
} |
} |
141 |
if (fabs (1.0 - prob) < GSL_DBL_EPSILON) |
|
142 |
|
{ |
143 |
|
double dp, phi; |
144 |
|
|
145 |
|
start: |
146 |
|
dp = -(p - gsl_cdf_tdist_Q (x, nu)); |
147 |
|
phi = gsl_ran_tdist_pdf (x, nu); |
148 |
|
|
149 |
|
if (dp == 0.0) |
150 |
|
goto end; |
151 |
|
|
152 |
{ |
{ |
153 |
return GSL_NEGINF; |
double lambda = dp / phi; |
154 |
} |
double step0 = lambda; |
155 |
printf ("prob is %f\t", prob); |
double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0); |
156 |
tmp = nu / 2.0; |
|
157 |
d = tmp * M_SQRTPI * gsl_sf_gamma (tmp) / gsl_sf_gamma (tmp + 0.5); |
double step = step0; |
158 |
method_test = gsl_max (d * prob, 2.0 / nu); |
|
159 |
/* |
if (fabs (step1) < fabs (step0)) |
160 |
* There are two possible initial approximations. |
step += step1; |
161 |
* Which is used depends on prob and nu. |
|
162 |
*/ |
printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu)); |
163 |
if (method_test > 0.05) |
|
164 |
{ |
x += step; |
165 |
tmp = prob; |
|
166 |
a = nu - 0.5; |
if (fabs (step) > 1e-10 * fabs (x)) |
167 |
/* gsl_cdf_ugaussian_Pinv(tmp) ? */ |
goto start; |
|
x = gsl_cdf_ugaussian_Pinv (tmp); |
|
|
y = inv_cornish_fisher (x, nu); |
|
|
tsqn = -1.0 + exp (a * y * y); |
|
168 |
} |
} |
169 |
else |
} |
170 |
{ |
|
171 |
z = pow (prob * d, 1 / tmp); |
end: |
172 |
zz = z * z; |
|
173 |
zzz = z * z * z; |
return x; |
|
tsqn = 1 / z + (nu + 1.0) * (-1.0 + z / (2.0 * (nu + 4.0)) + |
|
|
nu * zz / (3.0 * (nu + 2.0) * (nu + 6.0)) |
|
|
+ nu * (nu + 3.0) * (2.0 * nu * nu + |
|
|
9.0 * nu - |
|
|
2.0) * zzz / (8.0 * |
|
|
(nu + |
|
|
2.0) * |
|
|
(nu + |
|
|
2.0) * |
|
|
(nu + |
|
|
4.0) * |
|
|
(nu + |
|
|
4.0) * |
|
|
(nu + |
|
|
8.0))) |
|
|
/ (nu + 2); |
|
|
|
|
|
} |
|
|
initial_result = sqrt (tsqn * nu); |
|
|
printf ("initial_result = %f \t", initial_result); |
|
|
tcdf = gsl_cdf_tdist_Q (initial_result, nu); |
|
|
tmp = 0.5 * (tcdf - prob); |
|
|
tmp2 = tdist_pdf (initial_result, nu); |
|
|
w = tmp / tmp2; |
|
|
psi = initial_result * (nu + 1.0) / (nu + initial_result * initial_result); |
|
|
psi_prime = (nu + 1.0) * (nu - initial_result * initial_result) / |
|
|
((nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result)); |
|
|
d_psi_dt = 2 * initial_result * nu * (nu + 1.0) / |
|
|
((nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result)); |
|
|
d_psiprime_dt = -2.0 * initial_result * (nu + 1.0) * |
|
|
(3.0 * nu - initial_result * initial_result) / |
|
|
((nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result)); |
|
|
c2 = psi / 2.0; |
|
|
c3 = (2.0 * psi * psi + psi_prime) / 6.0; |
|
|
c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) + |
|
|
4.0 * psi * d_psi_dt + d_psiprime_dt; |
|
|
c4 /= 24.0; |
|
|
result = initial_result + w + c2 * w * w + |
|
|
c3 * w * w * w + c4 * w * w * w * w; |
|
|
printf ("result = %f\n", result); |
|
|
return result; |
|
174 |
} |
} |
|
#endif |
|