/[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.14 by bjg, Fri Jun 20 09:19:45 2003 UTC revision 1.1.2.15 by bjg, Fri Jun 20 14:19:56 2003 UTC
# Line 44  Line 44 
44    
45  #define MAXI 40  #define MAXI 40
46    
47    static double
48    poly_eval (const double c[], unsigned int n, double x)
49    {
50      unsigned int i;
51      double y = c[0] * x;
52    
53      for (i = 1; i < n; i++)
54        {
55          y = x * (y + c[i]);
56        }
57    
58      y += c[n];
59    
60      return y;
61    }
62    
63  static double  static double
64  cornish_fisher (double t, double n)  cornish_fisher (double t, double n)
65  {  {
   double ret_val = 0.0;  
   double b;  
   double z, zsq;  
   double tmp;  
   double p[6];  
   int i;  
   
66    const double coeffs6[10] = {    const double coeffs6[10] = {
67      0.265974025974025974026,      0.265974025974025974026,
68      5.449696969696969696970,      5.449696969696969696970,
# Line 91  cornish_fisher (double t, double n) Line 100  cornish_fisher (double t, double n)
100      85.5      85.5
101    };    };
102    
103    tmp = n - 0.5;    double x = n - 0.5;
104    z = tmp * log (1 + (t * t / n));    double b = 48.0 * x * x;
105    z = sqrt (z);    
106    b = 48.0 * tmp * tmp;    double z2 = x * log (1 + (x * x / n));
107    zsq = z * z;    double z = sqrt (z2);
108      
109    p[5] = z * poly_eval(coeffs6, 9, zsq);    double p5 = z * poly_eval(coeffs6, 9, z2);
110    p[4] = z * poly_eval(coeffs5, 7, zsq);    double p4 = z * poly_eval(coeffs5, 7, z2);
111    p[3] = z * poly_eval(coeffs4, 5, zsq);    double p3 = z * poly_eval(coeffs4, 5, z2);
112    p[2] = z * poly_eval(coeffs3, 5, zsq);    double p2 = z * poly_eval(coeffs3, 3, z2);
113      double p1 = z * (z2 + 3.0);
114    p[5] = coeffs6[0] * zsq;    double p0 = z;
115    for (i = 1; i < 9; i++)    
116      {    double y = p5;
117        p[5] = (p[5] + coeffs6[i]) * zsq;    y = (y / b) + p4;
118      }    y = (y / b) + p3;
119    p[5] = z * (p[5] + coeffs6[9]);    y = (y / b) + p2;
120      y = (y / b) + p1;
121    p[4] = coeffs5[0] * zsq;    y = (y / b) + p0;
   for (i = 1; i < 7; i++)  
     {  
       p[4] = (p[4] + coeffs5[i]) * zsq;  
     }  
   p[4] = z * (p[4] + coeffs5[7]);  
   
   p[3] = coeffs4[0] * zsq;  
   for (i = 1; i < 5; i++)  
     {  
       p[3] = (p[3] + coeffs4[i]) * zsq;  
     }  
   p[3] = z * (p[3] + coeffs4[5]);  
   
   p[2] = coeffs3[0] * zsq;  
   for (i = 0; i < 3; i++)  
     {  
       p[2] = (p[2] + coeffs3[i]) * zsq;  
     }  
   p[2] = z * (p[2] + coeffs3[1]);  /* BJG: should this be a [3]? */  
   
   p[1] = z * (zsq + 3.0);  
   p[0] = z;  
122    
123    ret_val = p[5];    return y;
   for (i = 4; i > -1; i--)  
     {  
       ret_val = (ret_val / b) + p[i];  
     }  
   return ret_val;  
124  }  }
125    
126  static double  static double
# Line 160  tdist_pdf (const double x, const double Line 142  tdist_pdf (const double x, const double
142   * and small degrees of freedom.   * and small degrees of freedom.
143   */   */
144  static double  static double
145  gsl_cdf_t_smallx (const double t, const double nu)  t_smallx (const double t, const double nu)
146  {  {
147    double y;    double x = nu / (nu + t*t);
148    double q;    double q =  1 - gsl_sf_beta_inc (nu / 2.0, 0.5, x);
   int i;  
   double b;  
   double diff;  
   double ckp2;  
   double ck;  
   int k;  
   int inu = (int) nu;  
   
   y = t / sqrt (nu);  
   diff = fabs (nu - 1.0);  
   if (diff < GSL_DBL_EPSILON)  
     {  
       q = M_PI_2 * atan (y);  
     }  
   else if ((nu < 21.0) && (nu > 1.0))  
     {  
       ckp2 = 1.0;  
       b = 1.0 + t * t / nu;  
149    
       for (k = inu - 2; k > 1; k -= 2)  
         {  
           ck = 1.0 + (ckp2 * (k - 1)) / ((double) k * b);  
           ckp2 = ck;  
         }  
       if ((inu % 2) == 0)  
         {  
           q = ck * y / sqrt (b);  
         }  
       else  
         {  
           q = M_PI_2 * (atan (y) + ck * y / b);  
         }  
     }  
150    return q;    return q;
151  }  }
152    
# Line 233  normal_approx (const double x, const dou Line 183  normal_approx (const double x, const dou
183    double diff;    double diff;
184    double q;    double q;
185    int i;    int i;
186    int rc;    double lg1, lg2;
   
   gsl_sf_result lg1;  
   gsl_sf_result lg2;  
187    
188    y = 1 / sqrt (1 + x * x / nu);    y = 1 / sqrt (1 + x * x / nu);
189    num = 1.0;    num = 1.0;
# Line 250  normal_approx (const double x, const dou Line 197  normal_approx (const double x, const dou
197        diff = q - diff;        diff = q - diff;
198      }      }
199    q += 1 / nu;    q += 1 / nu;
200    rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1);    lg1 = gsl_sf_lngamma (nu / 2.0);
201    if (rc == GSL_SUCCESS)    lg2 = gsl_sf_lngamma ((nu + 1.0) / 2.0);
202      {  
203        rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2);    diff = lg2 - lg1;
       if (rc != GSL_SUCCESS)  
         {  
           return rc;  
         }  
     }  
   diff = (lg2.val) - (lg1.val);  
204    q *= pow (y, nu) * exp (diff) / sqrt (M_PI);    q *= pow (y, nu) * exp (diff) / sqrt (M_PI);
205    
206    return q;    return q;
# Line 268  normal_approx (const double x, const dou Line 209  normal_approx (const double x, const dou
209  double  double
210  gsl_cdf_t_P (const double x, const double nu)  gsl_cdf_t_P (const double x, const double nu)
211  {  {
212    double val = 0.0;    double P;
   double q;  
   double absx = fabs (x);  
213    
214    if (fabs (x) < 4.0)    if (fabs (x) < 4.0)
215      {      {
216        q = gsl_cdf_t_smallx (x, nu);        double q = t_smallx (fabs(x), nu);
217        if (x <= 0.0)  
218          if (x >= 0.0)
219          {          {
220            val = (1.0 - q) / 2.0;            P = (1.0 + q) / 2.0;
221          }          }
222        else        else
223          {          {
224            val = (1.0 + q) / 2.0;            P = (1.0 - q) / 2.0;
225          }          }
226      }      }
227    else if (nu >= 21.0)    else if (nu >= 21.0)
228      {      {
229        q = t_cornish_fisher (x, nu);        double q = t_cornish_fisher (x, nu);
230    
231        if (x <= 0.0)        if (x <= 0.0)
232          {          {
233            val = q;            P = q;
234          }          }
235        else        else
236          {          {
237            val = 1.0 - q;            P = 1.0 - q;
238          }          }
239      }      }
240    else    else
241      {      {
242        val = normal_approx (x, nu);        P = normal_approx (x, nu);
243    
244        if (x > 0.0)        if (x > 0.0)
245          {          {
246            val = 1.0 - val;            P = 1.0 - P;
247          }          }
248      }      }
249    return val;  
250      return P;
251  }  }
252    
253  double  double
# Line 316  gsl_cdf_t_Q (const double x, const doubl Line 259  gsl_cdf_t_Q (const double x, const doubl
259    
260    if (fabs (x) < 4.0)    if (fabs (x) < 4.0)
261      {      {
262        q = gsl_cdf_t_smallx (x, nu);        q = t_smallx (x, nu);
263        if (x > 0.0)        if (x > 0.0)
264          {          {
265            val = (1.0 - q) / 2.0;            val = (1.0 - q) / 2.0;

Legend:
Removed from v.1.1.2.14  
changed lines
  Added in v.1.1.2.15

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