44 |
|
|
45 |
#define MAXI 40 |
#define MAXI 40 |
46 |
|
|
47 |
|
static double |
48 |
|
poly_eval (const double c[], unsigned int n, double x) |
49 |
|
{ |
50 |
|
unsigned int i; |
51 |
|
double y = c[0] * x; |
52 |
|
|
53 |
|
for (i = 1; i < n; i++) |
54 |
|
{ |
55 |
|
y = x * (y + c[i]); |
56 |
|
} |
57 |
|
|
58 |
|
y += c[n]; |
59 |
|
|
60 |
|
return y; |
61 |
|
} |
62 |
|
|
63 |
static double |
static double |
64 |
cornish_fisher (double t, double n) |
cornish_fisher (double t, double n) |
65 |
{ |
{ |
|
double ret_val = 0.0; |
|
|
double b; |
|
|
double z, zsq; |
|
|
double tmp; |
|
|
double p[6]; |
|
|
int i; |
|
|
|
|
66 |
const double coeffs6[10] = { |
const double coeffs6[10] = { |
67 |
0.265974025974025974026, |
0.265974025974025974026, |
68 |
5.449696969696969696970, |
5.449696969696969696970, |
100 |
85.5 |
85.5 |
101 |
}; |
}; |
102 |
|
|
103 |
tmp = n - 0.5; |
double x = n - 0.5; |
104 |
z = tmp * log (1 + (t * t / n)); |
double b = 48.0 * x * x; |
105 |
z = sqrt (z); |
|
106 |
b = 48.0 * tmp * tmp; |
double z2 = x * log (1 + (x * x / n)); |
107 |
zsq = z * z; |
double z = sqrt (z2); |
108 |
|
|
109 |
p[5] = z * poly_eval(coeffs6, 9, zsq); |
double p5 = z * poly_eval(coeffs6, 9, z2); |
110 |
p[4] = z * poly_eval(coeffs5, 7, zsq); |
double p4 = z * poly_eval(coeffs5, 7, z2); |
111 |
p[3] = z * poly_eval(coeffs4, 5, zsq); |
double p3 = z * poly_eval(coeffs4, 5, z2); |
112 |
p[2] = z * poly_eval(coeffs3, 5, zsq); |
double p2 = z * poly_eval(coeffs3, 3, z2); |
113 |
|
double p1 = z * (z2 + 3.0); |
114 |
p[5] = coeffs6[0] * zsq; |
double p0 = z; |
115 |
for (i = 1; i < 9; i++) |
|
116 |
{ |
double y = p5; |
117 |
p[5] = (p[5] + coeffs6[i]) * zsq; |
y = (y / b) + p4; |
118 |
} |
y = (y / b) + p3; |
119 |
p[5] = z * (p[5] + coeffs6[9]); |
y = (y / b) + p2; |
120 |
|
y = (y / b) + p1; |
121 |
p[4] = coeffs5[0] * zsq; |
y = (y / b) + p0; |
|
for (i = 1; i < 7; i++) |
|
|
{ |
|
|
p[4] = (p[4] + coeffs5[i]) * zsq; |
|
|
} |
|
|
p[4] = z * (p[4] + coeffs5[7]); |
|
|
|
|
|
p[3] = coeffs4[0] * zsq; |
|
|
for (i = 1; i < 5; i++) |
|
|
{ |
|
|
p[3] = (p[3] + coeffs4[i]) * zsq; |
|
|
} |
|
|
p[3] = z * (p[3] + coeffs4[5]); |
|
|
|
|
|
p[2] = coeffs3[0] * zsq; |
|
|
for (i = 0; i < 3; i++) |
|
|
{ |
|
|
p[2] = (p[2] + coeffs3[i]) * zsq; |
|
|
} |
|
|
p[2] = z * (p[2] + coeffs3[1]); /* BJG: should this be a [3]? */ |
|
|
|
|
|
p[1] = z * (zsq + 3.0); |
|
|
p[0] = z; |
|
122 |
|
|
123 |
ret_val = p[5]; |
return y; |
|
for (i = 4; i > -1; i--) |
|
|
{ |
|
|
ret_val = (ret_val / b) + p[i]; |
|
|
} |
|
|
return ret_val; |
|
124 |
} |
} |
125 |
|
|
126 |
static double |
static double |
142 |
* and small degrees of freedom. |
* and small degrees of freedom. |
143 |
*/ |
*/ |
144 |
static double |
static double |
145 |
gsl_cdf_t_smallx (const double t, const double nu) |
t_smallx (const double t, const double nu) |
146 |
{ |
{ |
147 |
double y; |
double x = nu / (nu + t*t); |
148 |
double q; |
double q = 1 - gsl_sf_beta_inc (nu / 2.0, 0.5, x); |
|
int i; |
|
|
double b; |
|
|
double diff; |
|
|
double ckp2; |
|
|
double ck; |
|
|
int k; |
|
|
int inu = (int) nu; |
|
|
|
|
|
y = t / sqrt (nu); |
|
|
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; |
|
149 |
|
|
|
for (k = inu - 2; k > 1; k -= 2) |
|
|
{ |
|
|
ck = 1.0 + (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); |
|
|
} |
|
|
} |
|
150 |
return q; |
return q; |
151 |
} |
} |
152 |
|
|
183 |
double diff; |
double diff; |
184 |
double q; |
double q; |
185 |
int i; |
int i; |
186 |
int rc; |
double lg1, lg2; |
|
|
|
|
gsl_sf_result lg1; |
|
|
gsl_sf_result lg2; |
|
187 |
|
|
188 |
y = 1 / sqrt (1 + x * x / nu); |
y = 1 / sqrt (1 + x * x / nu); |
189 |
num = 1.0; |
num = 1.0; |
197 |
diff = q - diff; |
diff = q - diff; |
198 |
} |
} |
199 |
q += 1 / nu; |
q += 1 / nu; |
200 |
rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1); |
lg1 = gsl_sf_lngamma (nu / 2.0); |
201 |
if (rc == GSL_SUCCESS) |
lg2 = gsl_sf_lngamma ((nu + 1.0) / 2.0); |
202 |
{ |
|
203 |
rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2); |
diff = lg2 - lg1; |
|
if (rc != GSL_SUCCESS) |
|
|
{ |
|
|
return rc; |
|
|
} |
|
|
} |
|
|
diff = (lg2.val) - (lg1.val); |
|
204 |
q *= pow (y, nu) * exp (diff) / sqrt (M_PI); |
q *= pow (y, nu) * exp (diff) / sqrt (M_PI); |
205 |
|
|
206 |
return q; |
return q; |
209 |
double |
double |
210 |
gsl_cdf_t_P (const double x, const double nu) |
gsl_cdf_t_P (const double x, const double nu) |
211 |
{ |
{ |
212 |
double val = 0.0; |
double P; |
|
double q; |
|
|
double absx = fabs (x); |
|
213 |
|
|
214 |
if (fabs (x) < 4.0) |
if (fabs (x) < 4.0) |
215 |
{ |
{ |
216 |
q = gsl_cdf_t_smallx (x, nu); |
double q = t_smallx (fabs(x), nu); |
217 |
if (x <= 0.0) |
|
218 |
|
if (x >= 0.0) |
219 |
{ |
{ |
220 |
val = (1.0 - q) / 2.0; |
P = (1.0 + q) / 2.0; |
221 |
} |
} |
222 |
else |
else |
223 |
{ |
{ |
224 |
val = (1.0 + q) / 2.0; |
P = (1.0 - q) / 2.0; |
225 |
} |
} |
226 |
} |
} |
227 |
else if (nu >= 21.0) |
else if (nu >= 21.0) |
228 |
{ |
{ |
229 |
q = t_cornish_fisher (x, nu); |
double q = t_cornish_fisher (x, nu); |
230 |
|
|
231 |
if (x <= 0.0) |
if (x <= 0.0) |
232 |
{ |
{ |
233 |
val = q; |
P = q; |
234 |
} |
} |
235 |
else |
else |
236 |
{ |
{ |
237 |
val = 1.0 - q; |
P = 1.0 - q; |
238 |
} |
} |
239 |
} |
} |
240 |
else |
else |
241 |
{ |
{ |
242 |
val = normal_approx (x, nu); |
P = normal_approx (x, nu); |
243 |
|
|
244 |
if (x > 0.0) |
if (x > 0.0) |
245 |
{ |
{ |
246 |
val = 1.0 - val; |
P = 1.0 - P; |
247 |
} |
} |
248 |
} |
} |
249 |
return val; |
|
250 |
|
return P; |
251 |
} |
} |
252 |
|
|
253 |
double |
double |
259 |
|
|
260 |
if (fabs (x) < 4.0) |
if (fabs (x) < 4.0) |
261 |
{ |
{ |
262 |
q = gsl_cdf_t_smallx (x, nu); |
q = t_smallx (x, nu); |
263 |
if (x > 0.0) |
if (x > 0.0) |
264 |
{ |
{ |
265 |
val = (1.0 - q) / 2.0; |
val = (1.0 - q) / 2.0; |