167 |
double y; |
double y; |
168 |
double diff; |
double diff; |
169 |
double ckp2; |
double ckp2; |
170 |
|
double lg1; |
171 |
|
double lg2; |
172 |
|
double num; |
173 |
int k; |
int k; |
174 |
double u; |
double u; |
175 |
double p; |
double p; |
247 |
* 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 |
248 |
* to use two different series expansions. |
* to use two different series expansions. |
249 |
*/ |
*/ |
250 |
u = 1+t*t/df; |
|
251 |
q = 2 * tdist_pdf ( t, df) * sqrt( df * u); |
lg1 = gsl_sf_lngamma (df / 2); |
252 |
p = q / df; |
lg2 = gsl_sf_lngamma ((df + 1) / 2); |
253 |
for ( i = 2; i < 20; i++) |
y = 1/sqrt(1+t*t/df); |
254 |
|
num = y; |
255 |
|
p = 0.0; |
256 |
|
for ( i = 2; i < 20; i+=2) |
257 |
{ |
{ |
258 |
q *= (i-1)/(i*u); |
num *= ((i-1)/i); |
259 |
p += q / (df+i); |
p += num / (df+i); |
260 |
} |
} |
261 |
|
p += 1/df; |
262 |
|
p *= pow(y,idf) * exp (lg2 - lg1) / sqrt (M_PI); |
263 |
if ( tail == GSL_CDF_UPPER ) |
if ( tail == GSL_CDF_UPPER ) |
264 |
{ |
{ |
265 |
result->val = p; |
result->val = p; |