43 |
* needs my own errors to be fixed.) |
* needs my own errors to be fixed.) |
44 |
*/ |
*/ |
45 |
#include <gsl/gsl_cdf.h> |
#include <gsl/gsl_cdf.h> |
46 |
|
#include <gsl/gsl_sf.h> |
47 |
#include <gsl/gsl_math.h> |
#include <gsl/gsl_math.h> |
48 |
#include <gsl/gsl_errno.h> |
#include <gsl/gsl_errno.h> |
49 |
|
|
50 |
|
#define MAXI 40 |
51 |
|
|
52 |
static double cornish_fisher ( double t, double n ) |
static double cornish_fisher ( double t, double n ) |
53 |
{ |
{ |
54 |
double ret_val = 0.0; |
double ret_val = 0.0; |
150 |
|
|
151 |
/* |
/* |
152 |
* This approximation uses different methods depending on the |
* This approximation uses different methods depending on the |
153 |
* degrees of freedom (df) and the argument t. An alternate method |
* degrees of freedom (nu) and the argument t. An alternate method |
154 |
* could use the incomplete beta function. I didn't choose that |
* could use the incomplete beta function. I didn't choose that |
155 |
* method because Kennedy and Gentle state that it is usually inferior to |
* method because Kennedy and Gentle state that it is usually inferior to |
156 |
* the method shown here. But "Statistical Computing" was written |
* the method shown here. But "Statistical Computing" was written |
159 |
* this to be true, the method below can be replaced with that |
* this to be true, the method below can be replaced with that |
160 |
* of the incomplete beta function. |
* of the incomplete beta function. |
161 |
*/ |
*/ |
162 |
int gsl_cdf_t_e ( double t, double df, gsl_cdf_result * result, |
int gsl_cdf_t_e ( double t, double nu, gsl_cdf_result * result, |
163 |
gsl_cdf_tail_t tail ) |
gsl_cdf_tail_t tail ) |
164 |
{ |
{ |
165 |
int rc; |
int rc; |
170 |
double y; |
double y; |
171 |
double diff; |
double diff; |
172 |
double ckp2; |
double ckp2; |
173 |
double lg1; |
gsl_sf_result lg1; |
174 |
double lg2; |
gsl_sf_result lg2; |
175 |
double num; |
double num; |
176 |
int k; |
int k; |
177 |
double u; |
double u; |
178 |
double p; |
double p; |
179 |
int i; |
int i; |
180 |
int idf = (int) df; |
int inu = (int) nu; |
181 |
|
|
182 |
/* |
/* |
183 |
* First approximation is also found in Abramowitz |
* First approximation is also found in Abramowitz |
184 |
* and Stegun. This is used only for small values of t |
* and Stegun. This is used only for small values of t |
185 |
* and small degrees of freedom. |
* and small degrees of freedom. |
186 |
*/ |
*/ |
187 |
y = t / sqrt(df); |
y = t / sqrt(nu); |
188 |
if ( fabs(t) < 4.0 ) |
if ( fabs(t) < 4.0 ) |
189 |
{ |
{ |
190 |
diff = fabs(df-1.0); |
diff = fabs(nu-1.0); |
191 |
if ( diff < GSL_DBL_EPSILON ) |
if ( diff < GSL_DBL_EPSILON ) |
192 |
{ |
{ |
193 |
q = M_PI_2 * atan(y); |
q = M_PI_2 * atan(y); |
194 |
} |
} |
195 |
else if ( (df < 21.0) && (df > 1.0)) |
else if ( (nu < 21.0) && (nu > 1.0)) |
196 |
{ |
{ |
197 |
ckp2 = 1.0; |
ckp2 = 1.0; |
198 |
b = 1.0 + t*t/df; |
b = 1.0 + t*t/nu; |
199 |
|
|
200 |
for ( k = idf-2; k > 1; k-=2) |
for ( k = inu-2; k > 1; k-=2) |
201 |
{ |
{ |
202 |
ck = 1 + (ckp2 * (k-1)) / ((double)k*b); |
ck = 1 + (ckp2 * (k-1)) / ((double)k*b); |
203 |
ckp2 = ck; |
ckp2 = ck; |
204 |
} |
} |
205 |
if ( (idf%2) == 0 ) |
if ( (inu%2) == 0 ) |
206 |
{ |
{ |
207 |
q = ck * y / sqrt(b); |
q = ck * y / sqrt(b); |
208 |
} |
} |
219 |
result->val = (1.0+q)/2; |
result->val = (1.0+q)/2; |
220 |
} |
} |
221 |
} |
} |
222 |
else if ( df >= 21.0 ) |
else if ( nu >= 21.0 ) |
223 |
{ |
{ |
224 |
/* |
/* |
225 |
* Use the Cornish-Fisher expansion to find |
* Use the Cornish-Fisher expansion to find |
226 |
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
227 |
* Approximate the t cdf with gsl_cdf_gauss. |
* Approximate the t cdf with gsl_cdf_gauss. |
228 |
*/ |
*/ |
229 |
u = cornish_fisher ( t, df ); |
u = cornish_fisher ( t, nu ); |
230 |
q = gsl_cdf_gauss ( u, GSL_CDF_UPPER ); |
q = gsl_cdf_gauss ( u, GSL_CDF_UPPER ); |
231 |
if ( tail == GSL_CDF_UPPER ) |
if ( tail == GSL_CDF_UPPER ) |
232 |
{ |
{ |
240 |
else |
else |
241 |
{ |
{ |
242 |
GSL_ERROR_VAL ("degrees of freedom are less than 1", |
GSL_ERROR_VAL ("degrees of freedom are less than 1", |
243 |
GSL_EINVAL, df ); |
GSL_EINVAL, nu ); |
244 |
} |
} |
245 |
} |
} |
246 |
else |
else |
251 |
* to use two different series expansions. |
* to use two different series expansions. |
252 |
*/ |
*/ |
253 |
|
|
254 |
lg1 = gsl_sf_lngamma (df / 2); |
|
255 |
lg2 = gsl_sf_lngamma ((df + 1) / 2); |
y = 1/sqrt(1+t*t/nu); |
|
y = 1/sqrt(1+t*t/df); |
|
256 |
num = y; |
num = y; |
257 |
p = 0.0; |
p = 0.0; |
258 |
for ( i = 2; i < 20; i+=2) |
|
259 |
|
for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2) |
260 |
{ |
{ |
261 |
|
diff = p; |
262 |
num *= ((i-1)/i); |
num *= ((i-1)/i); |
263 |
p += num / (df+i); |
p += num / (nu+i); |
264 |
|
diff = p - diff; |
265 |
|
} |
266 |
|
p = p + 1/nu; |
267 |
|
rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1); |
268 |
|
if ( rc == GSL_SUCCESS ) |
269 |
|
{ |
270 |
|
rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2); |
271 |
|
if ( rc != GSL_SUCCESS ) |
272 |
|
{ |
273 |
|
return rc; |
274 |
|
} |
275 |
|
} |
276 |
|
else |
277 |
|
{ |
278 |
|
return rc; |
279 |
} |
} |
280 |
p += 1/df; |
diff = (lg2.val) - (lg1.val); |
281 |
p *= pow(y,idf) * exp (lg2 - lg1) / sqrt (M_PI); |
p = p * pow(y,nu) * exp (diff) / sqrt (M_PI); |
282 |
|
|
283 |
if ( tail == GSL_CDF_UPPER ) |
if ( tail == GSL_CDF_UPPER ) |
284 |
{ |
{ |
285 |
result->val = p; |
result->val = p; |
286 |
} |
} |
287 |
else |
else if ( tail == GSL_CDF_LOWER ) |
288 |
{ |
{ |
289 |
result->val = 1.0 - p; |
result->val = 1.0 - p; |
290 |
} |
} |