/* cdf/tdistinv.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. */ #if 0 /* * 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) { 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 gsl_cdf_ut_Pinv (double prob, double nu) { 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; double x; double w; double c2; double c3; double c4; if (prob < 0.0) { return GSL_EDOM; } if (prob > 1.0) { return GSL_EDOM; } if (nu < 0.0) { return GSL_EDOM; } if (fabs (prob) < GSL_DBL_EPSILON) { return GSL_POSINF; } if (fabs (1.0 - prob) < GSL_DBL_EPSILON) { return GSL_NEGINF; } printf ("prob is %f\t", prob); tmp = nu / 2.0; d = tmp * M_SQRTPI * gsl_sf_gamma (tmp) / gsl_sf_gamma (tmp + 0.5); method_test = gsl_max (d * prob, 2.0 / nu); /* * There are two possible initial approximations. * Which is used depends on prob and nu. */ if (method_test > 0.05) { tmp = prob; a = nu - 0.5; /* gsl_cdf_ugaussian_Pinv(tmp) ? */ x = gsl_cdf_ugaussian_Pinv (tmp); y = inv_cornish_fisher (x, nu); tsqn = -1.0 + exp (a * y * y); } else { z = pow (prob * d, 1 / tmp); zz = z * z; zzz = z * z * z; 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; } #endif