/[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.16 by bjg, Sun Jun 22 21:20:36 2003 UTC revision 1.1.2.17 by bjg, Mon Jun 23 14:05:03 2003 UTC
# Line 44  Line 44 
44    
45  #include "beta_inc.c"  #include "beta_inc.c"
46    
47    #if 0
48  /*  /*
49   * First approximation is also found in Abramowitz   * First approximation is also found in Abramowitz
50   * and Stegun. This is used only for small values of t   * and Stegun. This is used only for small values of t
# Line 52  Line 53 
53  static double  static double
54  t_small (const double t, const double nu)  t_small (const double t, const double nu)
55  {  {
56    /* Final result is always a small amount to be added to or    /* Compute beta_inc when t^2 < n. Final result is always a small
57       subtracted from 0.5 so we only need absolute accuracy of       amount to be added to or subtracted from 0.5 so we only need
58       O(DBL_EPS) -- this information could be used to optimise inside       absolute accuracy of O(DBL_EPS) -- this information could be used
59       beta_inc (currently it isn't) */       to optimise inside beta_inc (currently it isn't) */
60    
61    double u = (t * t) / nu;    double u = (t * t) / nu;
62    double eps = u / (1 + u);    double eps = u / (1 + u);
# Line 67  t_small (const double t, const double nu Line 68  t_small (const double t, const double nu
68  static double  static double
69  t_large (const double t, const double nu)  t_large (const double t, const double nu)
70  {  {
71    /* Final result may be needed to full relative precision, depending    /* Compute beta_inc when t^2 > n.  Final result may be needed to full
72       on whether we are computing the upper or lower tail.  Again this       relative precision, depending on whether we are computing the
73       information could be used to optimise inside beta_inc.  */       upper or lower tail.  Again this information could be used to
74         optimise inside beta_inc.  */
75    
76    double v = nu / (t * t);    double v = nu / (t * t);
77    double eps = v / (1 + v);    double eps = v / (1 + v);
# Line 78  t_large (const double t, const double nu Line 80  t_large (const double t, const double nu
80    return q;    return q;
81  }  }
82    
83    #endif
84    
85  static double  static double
86  poly_eval (const double c[], unsigned int n, double x)  poly_eval (const double c[], unsigned int n, double x)
87  {  {
# Line 223  gsl_cdf_t_P (const double x, const doubl Line 227  gsl_cdf_t_P (const double x, const doubl
227    
228    double x2 = x * x;    double x2 = x * x;
229    
230    if (nu <= 20.0)    if (1 || nu <= 20.0)
231      {      {
232        if (x2 < nu)        if (x2 < nu)
233          {          {
234            double q = t_small (fabs (x), nu);            double u = x2 / nu;
235              double eps = u / (1 + u);
236    
237            if (x >= 0.0)            if (x >= 0)
238              {              {
239                P = 0.5 * (1.0 + q);                P = beta_inc_AXPY (0.5, 0.5, 0.5, nu / 2.0, eps);
240              }              }
241            else            else
242              {              {
243                P = 0.5 * (1.0 - q);                P = beta_inc_AXPY (-0.5, 0.5, 0.5, nu / 2.0, eps);
244              }              }
245          }          }
246        else        else
247          {          {
248            double q = t_large (fabs (x), nu);            double v = nu / (x * x);
249              double eps = v / (1 + v);
250              
251              if (x >= 0)
252                {
253                  P = beta_inc_AXPY (-0.5, 1.0, nu / 2.0, 0.5, eps);
254                }
255              else
256                {
257                  P = beta_inc_AXPY (0.5, 0.0,  nu / 2.0, 0.5, eps);
258                }
259            }
260        }
261    #if 0
262      else
263        {
264          if (0 && x2 < nu)
265            {
266              double q = t_cornish_fisher (x, nu);
267              
268              if (x <= 0.0)
269                {
270                  P = q;
271                }
272              else
273                {
274                  P = 1.0 - q;
275                }
276            }
277          else if (x2 > nu)
278            {
279              double q = t_large (fabs (x), nu);
280              
281            if (x >= 0.0)            if (x >= 0.0)
282              {              {
283                P = 1 - 0.5 * q;                P = 1 - 0.5 * q;
# Line 252  gsl_cdf_t_P (const double x, const doubl Line 288  gsl_cdf_t_P (const double x, const doubl
288              }              }
289          }          }
290      }      }
291    else  #endif
     {  
       double q = t_cornish_fisher (x, nu);  
   
       if (x <= 0.0)  
         {  
           P = q;  
         }  
       else  
         {  
           P = 1.0 - q;  
         }  
     }  
292  #if 0  #if 0
293    {    {
294      P = normal_approx (x, nu);      P = normal_approx (x, nu);
# Line 279  gsl_cdf_t_P (const double x, const doubl Line 303  gsl_cdf_t_P (const double x, const doubl
303    return P;    return P;
304  }  }
305    
306    #if 0
307  double  double
308  gsl_cdf_t_Q (const double x, const double nu)  gsl_cdf_t_Q (const double x, const double nu)
309  {  {
# Line 495  gsl_cdf_ut_P_inv (double prob, double nu Line 520  gsl_cdf_ut_P_inv (double prob, double nu
520    printf ("result = %f\n", result);    printf ("result = %f\n", result);
521    return result;    return result;
522  }  }
523    #endif

Legend:
Removed from v.1.1.2.16  
changed lines
  Added in v.1.1.2.17

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