/[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.11 by jstover, Wed Apr 30 01:23:28 2003 UTC revision 1.1.2.12 by bjg, Wed May 28 14:30:35 2003 UTC
# Line 447  double gsl_cdf_t_Q ( const double x, con Line 447  double gsl_cdf_t_Q ( const double x, con
447      }      }
448    return val;    return val;
449  }  }
450    /*
451     * Invert the T distribution. Uses a method as shown in
452     * Statistical Computing, 5.4.2. This method uses an initial
453     * approximation based on Hill's method above, then improves
454     * the initial method with a Taylor series.
455     */
456    double gsl_cdf_ut_P_inv (double prob, double nu)
457    {
458      double initial_result;
459      double result;
460      double d;
461      double tmp;
462      double tmp2;
463      double z;
464      double zz;
465      double zzz;
466      double y;
467      double tsqn;
468      double a;
469      double method_test;
470      double psi;
471      double psi_prime;
472      double d_psi_dt;
473      double d_psiprime_dt;
474      double tcdf;
475      double x;
476      double w;
477      double c2;
478      double c3;
479      double c4;
480    
481      if ( prob < 0.0 )
482        {
483          return GSL_EDOM;
484        }
485      if ( prob > 1.0 )
486        {
487          return GSL_EDOM;
488        }
489      if ( nu < 0.0 )
490        {
491          return GSL_EDOM;
492        }
493      if ( fabs(prob) < GSL_DBL_EPSILON )
494        {
495          return GSL_POSINF;
496        }
497      if ( fabs ( 1.0 - prob ) < GSL_DBL_EPSILON )
498        {
499          return GSL_NEGINF;
500        }
501      printf("prob is %f\t",prob);
502      tmp = nu / 2.0;
503      d = tmp * M_SQRTPI * gsl_sf_gamma ( tmp )  
504        / gsl_sf_gamma ( tmp + 0.5 );
505      method_test = gsl_max ( d*prob, 2.0 / nu);
506      /*
507       * There are two possible initial approximations.
508       * Which is used depends on prob and nu.
509       */
510      if ( method_test > 0.05 )
511        {
512          tmp = prob / 2.0;
513          a = nu - 0.5;
514          /* gsl_cdf_ugaussian_P_inv(tmp) ? */
515          x = gsl_cdf_ugaussian_Q_inv (tmp );
516          y = cornish_fisher ( x, nu );
517          tsqn = -1.0 + exp ( a*y*y);
518        }
519      else
520        {
521          z = pow ( prob * d, 1 / tmp );
522          zz = z * z;
523          zzz = z*z*z;
524          tsqn = 1/z + (nu + 1.0)*(-1.0 + z/(2.0 * (nu + 4.0)) +
525                                 nu * zz / (3.0 * (nu+2.0)*(nu+6.0))
526                                 + nu*(nu+3.0)*(2.0*nu*nu+9.0*nu-2.0)*zzz
527                                 /(8.0 * (nu + 2.0) * (nu + 2.0)
528                                   * (nu + 4.0) * (nu + 4.0) *
529                                   (nu + 8.0)))
530            /(nu + 2);
531    
532        }
533      initial_result = sqrt(tsqn * nu);
534      printf("initial_result = %f \t",initial_result);
535      tcdf = gsl_cdf_t_Q( initial_result, nu );
536      tmp = 0.5 * ( tcdf - prob );
537      tmp2 = tdist_pdf ( initial_result, nu );
538      w = tmp / tmp2;
539      psi = initial_result * (nu+1.0) /
540        (nu + initial_result * initial_result);
541      psi_prime = (nu+1.0) * (nu - initial_result * initial_result) /
542        ((nu + initial_result*initial_result) *
543         (nu + initial_result * initial_result));
544      d_psi_dt = 2 * initial_result * nu * (nu + 1.0) /
545        ((nu + initial_result * initial_result) *
546         (nu + initial_result * initial_result));
547      d_psiprime_dt = -2.0 * initial_result * (nu + 1.0) *
548        (3.0 * nu - initial_result * initial_result) /
549        ((nu + initial_result * initial_result) *
550         (nu + initial_result * initial_result) *
551         (nu + initial_result * initial_result));
552      c2 = psi/2.0;
553      c3 = (2.0 * psi * psi + psi_prime) / 6.0;
554      c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) +
555        4.0 * psi * d_psi_dt + d_psiprime_dt;
556      c4 /= 24.0;
557      result = initial_result + w + c2 * w * w +
558        c3 * w * w * w + c4 * w * w * w * w;
559      printf("result = %f\n",result);
560      return result;
561    }
562    

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

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