/[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.3 by bjg, Tue Jun 24 17:29:11 2003 UTC revision 1.1.2.4 by bjg, Sat Jun 28 18:08:40 2003 UTC
# Line 60  poly_eval (const double c[], unsigned in Line 60  poly_eval (const double c[], unsigned in
60    return y;    return y;
61  }  }
62    
63    /*
64     * Use the Cornish-Fisher asymptotic expansion to find a point u such
65     * that gsl_cdf_gauss(y) = tcdf(t).
66     *
67     */
68    
69  static double  static double
70  cornish_fisher (double t, double n)  cornish_fisher (double t, double n)
71  {  {
# Line 94  cornish_fisher (double t, double n) Line 100  cornish_fisher (double t, double n)
100      8518.5      8518.5
101    };    };
102    const double coeffs3[4] = {    const double coeffs3[4] = {
103      .4,      0.4,
104      3.3,      3.3,
105      24.0,      24.0,
106      85.5      85.5
107    };    };
108    
109    double m = n - 0.5;    double a = n - 0.5;
110    double b = 48.0 * m * m;    double b = 48.0 * a * a;
111    
112    double z2 = m * log (1 + (t * t / n));    double z2 = a * log (1 + (t * t / n));
113    double z = sqrt (z2);    double z = sqrt (z2);
114    
115    double p5 = z * poly_eval (coeffs6, 9, z2);    double p5 = z * poly_eval (coeffs6, 9, z2);
116    double p4 = z * poly_eval (coeffs5, 7, z2);    double p4 = -z * poly_eval (coeffs5, 7, z2);
117    double p3 = z * poly_eval (coeffs4, 5, z2);    double p3 = z * poly_eval (coeffs4, 5, z2);
118    double p2 = z * poly_eval (coeffs3, 3, z2);    double p2 = -z * poly_eval (coeffs3, 3, z2);
119    double p1 = z * (z2 + 3.0);    double p1 = z * (z2 + 3.0);
120    double p0 = z;    double p0 = z;
121    
# Line 120  cornish_fisher (double t, double n) Line 126  cornish_fisher (double t, double n)
126    y = (y / b) + p1;    y = (y / b) + p1;
127    y = (y / b) + p0;    y = (y / b) + p0;
128    
129    return y;    if (t < 0) y *= -1;
 }  
   
130    
131  /*    return y;
  * 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.  
  *  
  * Valid for: ????  
  *  
  */  
 static double  
 t_cornish_fisher (const double x, const double nu)  
 {  
   double u = cornish_fisher (x, nu);  
   double q = gsl_cdf_ugaussian_Q (u);  
   
   return q;  
132  }  }
133    
134    #if 0
135  /*  /*
136   * Series approximation for t > 4.0. This needs to be fixed;   * Series approximation for t > 4.0. This needs to be fixed;
137   * it shouldn't subtract the result from 1.0. A better way is   * it shouldn't subtract the result from 1.0. A better way is
# Line 181  normal_approx (const double x, const dou Line 172  normal_approx (const double x, const dou
172    
173    return q;    return q;
174  }  }
175    #endif
176    
177  double  double
178  gsl_cdf_tdist_P (const double x, const double nu)  gsl_cdf_tdist_P (const double x, const double nu)
# Line 189  gsl_cdf_tdist_P (const double x, const d Line 181  gsl_cdf_tdist_P (const double x, const d
181    
182    double x2 = x * x;    double x2 = x * x;
183    
184    if (x2 < nu)    if (nu > 30 && x2 < 10 * nu)
185        {
186          double u = cornish_fisher (x, nu);
187          P = gsl_cdf_ugaussian_P (u);
188    
189          return P;
190        }
191    
192     if (x2 < nu)
193      {      {
194        double u = x2 / nu;        double u = x2 / nu;
195        double eps = u / (1 + u);        double eps = u / (1 + u);
# Line 264  gsl_cdf_tdist_Q (const double x, const d Line 264  gsl_cdf_tdist_Q (const double x, const d
264    
265    
266  #if 0  #if 0
   else  
     {  
       if (0 && x2 < nu)  
         {  
           double q = t_cornish_fisher (x, nu);  
             
           if (x <= 0.0)  
             {  
               P = q;  
             }  
           else  
             {  
               P = 1.0 - q;  
             }  
         }  
       else if (x2 > nu)  
         {  
           double q = t_large (fabs (x), nu);  
             
           if (x >= 0.0)  
             {  
               P = 1 - 0.5 * q;  
             }  
           else  
             {  
               P = 0.5 * q;  
             }  
         }  
     }  
 #endif  
   
 #if 0  
   {  
     P = normal_approx (x, nu);  
   
     if (x > 0.0)  
       {  
         P = 1.0 - P;  
       }  
   }  
   
   return P;  
 }  
 #endif  
   
   
   
 #if 0  
267  /*  /*
268   * Invert the T distribution. Uses a method as shown in   * Invert the T distribution. Uses a method as shown in
269   * Statistical Computing, 5.4.2. This method uses an initial   * Statistical Computing, 5.4.2. This method uses an initial
# Line 374  tdist_pdf (const double x, const double Line 326  tdist_pdf (const double x, const double
326  }  }
327    
328  double  double
329  gsl_cdf_ut_P_inv (double prob, double nu)  gsl_cdf_ut_Pinv (double prob, double nu)
330  {  {
331    double initial_result;    double initial_result;
332    double result;    double result;
# Line 431  gsl_cdf_ut_P_inv (double prob, double nu Line 383  gsl_cdf_ut_P_inv (double prob, double nu
383      {      {
384        tmp = prob;        tmp = prob;
385        a = nu - 0.5;        a = nu - 0.5;
386        /* gsl_cdf_ugaussian_P_inv(tmp) ? */        /* gsl_cdf_ugaussian_Pinv(tmp) ? */
387        x = gsl_cdf_ugaussian_P_inv (tmp);        x = gsl_cdf_ugaussian_Pinv (tmp);
388        y = inv_cornish_fisher (x, nu);        y = inv_cornish_fisher (x, nu);
389        tsqn = -1.0 + exp (a * y * y);        tsqn = -1.0 + exp (a * y * y);
390      }      }

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

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