/[gsl]/gsl/cdf/tdist.c
ViewVC logotype

Diff of /gsl/cdf/tdist.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.4 by bjg, Sat Jun 28 18:08:40 2003 UTC revision 1.1.2.5 by bjg, Sun Jun 29 13:33:36 2003 UTC
# Line 1  Line 1 
1  /* cdf/t.c  /* cdf/tdist.c
2   *   *
3   * Copyright (C) 2002 Jason H. Stover.   * Copyright (C) 2002 Jason H. Stover.
4   *   *
# Line 263  gsl_cdf_tdist_Q (const double x, const d Line 263  gsl_cdf_tdist_Q (const double x, const d
263    
264    
265    
 #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  

Legend:
Removed from v.1.1.2.4  
changed lines
  Added in v.1.1.2.5

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26