/* cdf/t.c * * Copyright (C) 2002 Jason H. Stover. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. */ /* * Computes the Student's t cumulative distribution function using * the method detailed in * * W.J. Kennedy and J.E. Gentle, "Statistical Computing." 1980. * Marcel Dekker. ISBN 0-8247-6898-1. * * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions * of Cornish-Fisher type." Annals of Mathematical Statistics, * vol. 39, 1264-1273. 1968. * * G.W. Hill. "Algorithm 395: Student's t-distribution," Communications * of the ACM, volume 13, number 10, page 617. October 1970. * * G.W. Hill, "Remark on algorithm 395: Student's t-distribution," * Transactions on Mathematical Software, volume 7, number 2, page * 247. June 1982. */ /* * Currently this approximation is giving the wrong answers. * Not sure why; Hill's method may have to be dropped in favor * of use of the incomplete beta function. (Or perhaps it just * needs my own errors to be fixed.) */ #include #include #include static double cornish_fisher ( double t, double n ) { double ret_val = 0.0; double b; double z,zsq; double tmp; double p[6]; int i; const double coeffs6[10] = { 0.265974025974025974026, 5.449696969696969696970, 122.20294372294372294372, 2354.7298701298701298701, 37625.00902597402597403, 486996.1392857142857143, 4960870.65, 37978595.55, 201505390.875, 622437908.625}; const double coeffs5[8] = { 0.2742857142857142857142, 4.499047619047619047619, 78.45142857142857142857, 1118.710714285714285714, 12387.6, 101024.55, 559494.0, 1764959.625}; const double coeffs4[6] ={ 0.3047619047619047619048, 3.752380952380952380952, 46.67142857142857142857, 427.5, 2587.5, 8518.5}; const double coeffs3[4] = { .4, 3.3, 24.0, 85.5}; tmp = n-0.5; z = tmp*log(1+(t*t/n)); z = sqrt(z); b = 48.0 * sqrt(tmp); zsq = z*z; p[5] = coeffs6[0] * zsq; for ( i = 1; i < 9; i++) { p[5] = (p[5] + coeffs6[i])*zsq; } p[5] = z * (p[5] + coeffs6[9]); p[4] = coeffs5[0] * zsq; 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]); p[1] = z * (zsq +3.0); p[0] = z; ret_val = p[5]; for( i = 4; i > -1; i-- ) { ret_val = (ret_val/b)+p[i]; } return ret_val; } static double tdist_pdf (const double x, const double nu) { 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; } /* * This approximation uses different methods depending on the * degrees of freedom (df) 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 df, gsl_cdf_result * result, gsl_cdf_tail_t tail ) { int rc; int stopval; double b; double ck; double q; double y; double diff; double ckp2; int k; double u; double p; int i; /* * 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(df); if ( fabs(t) < 4.0 ) { diff = fabs(df-1.0); if ( diff < GSL_DBL_EPSILON ) { q = M_PI_2 * atan(y); } else if ( (df < 21.0) && (df > 1.0)) { ckp2 = 1.0; b = 1.0 + t*t/df; for ( k = (int)(df-2); k > 1; k-=2) { ck = 1 + (ckp2 * (k-1)) / ((double)k*b); ckp2 = ck; } if ( ((int)df%2) == 0 ) { q = ck * y / sqrt(b); } else { q = M_PI_2 * atan(y) + ck * y / b; } if ( tail == GSL_CDF_UPPER ) { result->val = 1.0 - q; } else { result->val = q; } } else if ( df >= 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, df ); q = gsl_cdf_gauss ( u, GSL_CDF_UPPER ); if ( tail == GSL_CDF_UPPER ) { result->val = q; } else { result->val = 1.0 - q; } } else { GSL_ERROR_VAL ("degrees of freedom are less than 1", GSL_EINVAL, df ); } } else { /* * Series approximation for t > 4.0. This needs to be fixed; * it shouldn't subtract the result from 1.0. A better way is * to use two different series expansions. */ u = 1+t*t/df; q = 2 * tdist_pdf ( t, df) * sqrt( df * u); p = q / df; for ( i = 2; i < 20; i++) { q *= (i-1)/(i*u); p += q / (df+i); } if ( tail == GSL_CDF_UPPER ) { result->val = p; } else { result->val = 1.0 - p; } } return rc; }