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

Diff of /gsl/cdf/t.c

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

revision 1.1.2.12 by bjg, Wed May 28 14:30:35 2003 UTC revision 1.1.2.13 by jstover, Sun Jun 1 13:58:27 2003 UTC
# Line 451  double gsl_cdf_t_Q ( const double x, con Line 451  double gsl_cdf_t_Q ( const double x, con
451   * Invert the T distribution. Uses a method as shown in   * Invert the T distribution. Uses a method as shown in
452   * Statistical Computing, 5.4.2. This method uses an initial   * Statistical Computing, 5.4.2. This method uses an initial
453   * approximation based on Hill's method above, then improves   * approximation based on Hill's method above, then improves
454   * the initial method with a Taylor series.   * the initial method with a Taylor series or a Cornish-Fisher
455     * expansion.
456   */   */
457    static double inv_cornish_fisher ( double z, double n )
458    {
459      double ret_val = 0.0;
460      double b;
461      double zsq;
462      double tmp;
463      double c;
464      double u;
465      double d;
466      double zcube;
467      double zfour;
468      double zfive;
469      double zseven;
470      int i;
471    
472      tmp = n-0.5;
473      b = 48.0 * tmp * tmp;
474      if ( n > 5.0 )
475        {
476          c = 96.36 - 16.0 / tmp - 98.0 / (tmp*tmp) + 20700.0 / (tmp*tmp*tmp*b);
477        }
478      else
479        {
480          c = 0.3 * (n - 4.5) * (z + 0.6);
481        }
482      zsq = z*z;
483      zcube = zsq*z;
484      zfour = zcube*z;
485      zfive = zfour*z;
486      zseven = zfour*zcube;
487      d = n * M_SQRTPI * gsl_sf_gamma ( tmp )  
488        / (2.0*gsl_sf_gamma ( tmp + 0.5 ));
489      u = 10.0 * b * ( b + c - 2.0*z - 7.0 * zsq - 5.0*zcube + 0.05*d*zfour);
490      tmp = 4.0*zseven + 63.0*zfive + 360.0*zcube + 945.0*z;
491      ret_val = z - (zcube+3.0*z)/b + tmp/u;
492      return ret_val;
493    }
494    
495  double gsl_cdf_ut_P_inv (double prob, double nu)  double gsl_cdf_ut_P_inv (double prob, double nu)
496  {  {
497    double initial_result;    double initial_result;
# Line 509  double gsl_cdf_ut_P_inv (double prob, do Line 548  double gsl_cdf_ut_P_inv (double prob, do
548     */     */
549    if ( method_test > 0.05 )    if ( method_test > 0.05 )
550      {      {
551        tmp = prob / 2.0;        tmp = prob;
552        a = nu - 0.5;        a = nu - 0.5;
553        /* gsl_cdf_ugaussian_P_inv(tmp) ? */        /* gsl_cdf_ugaussian_P_inv(tmp) ? */
554        x = gsl_cdf_ugaussian_Q_inv (tmp );        x = gsl_cdf_ugaussian_P_inv (tmp );
555        y = cornish_fisher ( x, nu );        y = inv_cornish_fisher ( x, nu );
556        tsqn = -1.0 + exp ( a*y*y);        tsqn = -1.0 + exp ( a*y*y);
557      }      }
558    else    else
# Line 559  double gsl_cdf_ut_P_inv (double prob, do Line 598  double gsl_cdf_ut_P_inv (double prob, do
598    printf("result = %f\n",result);    printf("result = %f\n",result);
599    return result;    return result;
600  }  }
   

Legend:
Removed from v.1.1.2.12  
changed lines
  Added in v.1.1.2.13

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