/[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.8 by bjg, Sat Mar 1 20:21:18 2003 UTC revision 1.1.2.9 by jstover, Sun Mar 9 02:51:01 2003 UTC
# Line 160  static double tdist_pdf (const double x, Line 160  static double tdist_pdf (const double x,
160   * this to be true, the method below can be replaced with that   * this to be true, the method below can be replaced with that
161   * of the incomplete beta function.   * of the incomplete beta function.
162   */   */
163  int gsl_cdf_t_e ( double t, double nu, gsl_cdf_result * result,  /* int gsl_cdf_t_e ( double t, double nu) */
164                    gsl_cdf_tail_t tail )  /* { */
165    /*   int rc; */
166    /*   int stopval; */
167    /*   double b; */
168    /*   double ck; */
169    /*   double q; */
170    /*   double y; */
171    /*   double diff; */
172    /*   double ckp2; */
173    /*   gsl_sf_result lg1; */
174    /*   gsl_sf_result lg2; */
175    /*   double num; */
176    /*   int k; */
177    /*   double u; */
178    /*   double p; */
179    /*   int i; */
180    /*   int inu = (int) nu; */
181    
182    /*   /\* */
183    /*    * First approximation is also found in Abramowitz */
184    /*    * and Stegun. This is used only for small values of t */
185    /*    * and small degrees of freedom. */
186    /*    *\/ */
187    /*   y = t / sqrt(nu); */
188    /*   if ( fabs(t) < 4.0 ) */
189    /*     { */
190    /*       diff = fabs(nu-1.0); */
191    /*       if ( diff < GSL_DBL_EPSILON ) */
192    /*      { */
193    /*        q = M_PI_2 * atan(y); */
194    /*      } */
195    /*       else if ( (nu < 21.0) && (nu > 1.0)) */
196    /*      { */
197    /*        ckp2 = 1.0; */
198    /*        b = 1.0 + t*t/nu; */
199              
200    /*        for ( k = inu-2; k > 1; k-=2) */
201    /*          { */
202    /*            ck = 1 + (ckp2 * (k-1)) / ((double)k*b); */
203    /*            ckp2 = ck; */
204    /*          } */
205    /*        if ( (inu%2) == 0 ) */
206    /*          { */
207    /*            q = ck * y / sqrt(b); */
208    /*          } */
209    /*        else */
210    /*          { */
211    /*            q = M_PI_2 * (atan(y) + ck * y / b); */
212    /*          } */
213    /*        if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */
214    /*          { */
215    /*            result->val = (1.0 - q)/2; */
216    /*          } */
217    /*        else */
218    /*          { */
219    /*            result->val = (1.0+q)/2; */
220    /*          } */
221    /*      } */
222    /*       else if ( nu >= 21.0 ) */
223    /*      { */
224              /*
225               * Use the Cornish-Fisher expansion to find
226               * a point u such that gsl_cdf_gauss(u) = tcdf(t).
227               * Approximate the t cdf with gsl_cdf_gauss.
228               */
229    /*        u = cornish_fisher ( t, nu ); */
230    /*        q = gsl_cdf_gauss_Q ( u ); */
231    /*        if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */
232    /*          { */
233    /*            result->val = q; */
234    /*          } */
235    /*        else */
236    /*          { */
237    /*            result->val = 1.0 - q; */
238    /*          } */
239    /*      } */
240    /*       else */
241    /*      { */
242    /*        GSL_ERROR_VAL ("degrees of freedom are less than 1", */
243    /*                       GSL_EINVAL, nu ); */
244    /*      } */
245    /*     } */
246    /*   return rc; */
247    /* } */
248    
249    /*
250     * First approximation is also found in Abramowitz
251     * and Stegun. This is used only for small values of t
252     * and small degrees of freedom.
253     */
254    static double gsl_cdf_t_smallx ( const double t, const double nu )
255  {  {
256    int rc;    double y;
257    int stopval;    double q;
258      int i;
259    double b;    double b;
260      double diff;
261      double ckp2;
262    double ck;    double ck;
263      int k;
264      int inu = (int) nu;
265    
266      y = t / sqrt(nu);
267      diff = fabs(nu-1.0);
268      if ( diff < GSL_DBL_EPSILON )
269        {
270          q = M_PI_2 * atan(y);
271        }
272      else if ( (nu < 21.0) && (nu > 1.0))
273        {
274          ckp2 = 1.0;
275          b = 1.0 + t*t/nu;
276              
277          for ( k = inu-2; k > 1; k-=2)
278            {
279              ck = 1.0 + (ckp2 * (k-1)) / ((double)k*b);
280              ckp2 = ck;
281            }
282          if ( (inu%2) == 0 )
283            {
284              q = ck * y / sqrt(b);
285            }
286          else
287            {
288              q = M_PI_2 * (atan(y) + ck * y / b);
289            }
290    /*       if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */
291    /*      { */
292    /*        result->val = (1.0 - q)/2; */
293    /*      } */
294    /*       else */
295    /*      { */
296    /*        result->val = (1.0+q)/2; */
297    /*      } */
298        }
299      return q;
300    }
301    /*
302     * Use the Cornish-Fisher expansion to find
303     * a point u such that gsl_cdf_gauss(u) = tcdf(t).
304     * Approximate the t cdf with gsl_cdf_gauss.
305     */
306    static double t_cornish_fisher ( const double x, const double nu )
307    {
308      double y;
309    double q;    double q;
310      double u;
311    
312      y = x / sqrt(nu);
313      u = cornish_fisher ( x, nu );
314      q = gsl_cdf_gauss_Q ( u);
315    /*   if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */
316    /*     { */
317    /*       result->val = q; */
318    /*     } */
319    /*   else */
320    /*     { */
321    /*       result->val = 1.0 - q; */
322    /*     } */
323      return q;
324    }
325    /*
326     * Series approximation for t > 4.0. This needs to be fixed;
327     * it shouldn't subtract the result from 1.0. A better way is
328     * to use two different series expansions. Figuring this out
329     * means rummaging through Fisher's paper in Metron, v5, 1926,
330     * "Expansion of Student's integral in powers of n^{-1}."
331     */
332    static double normal_approx ( const double x, const double nu )
333    {
334    double y;    double y;
335      double num;
336    double diff;    double diff;
337    double ckp2;    double q;
338      int i;
339      int rc;
340    
341    gsl_sf_result lg1;    gsl_sf_result lg1;
342    gsl_sf_result lg2;    gsl_sf_result lg2;
   double num;  
   int k;  
   double u;  
   double p;  
   int i;  
   int inu = (int) nu;  
343    
344    /*    y = 1/sqrt(1+x*x/nu);
345     * First approximation is also found in Abramowitz    num = 1.0;
346     * and Stegun. This is used only for small values of t    q = 0.0;
347     * and small degrees of freedom.    diff = 2*GSL_DBL_EPSILON;
348     */    for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2)
349    y = t / sqrt(nu);      {
350    if ( fabs(t) < 4.0 )        diff = q;
351          num *= y*y*(i-1)/i;
352          q += num / (nu+i);
353          diff = q - diff;
354        }
355      q += 1/nu;
356      rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1);
357      if ( rc == GSL_SUCCESS )
358      {      {
359        diff = fabs(nu-1.0);        rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2);
360        if ( diff < GSL_DBL_EPSILON )        if ( rc != GSL_SUCCESS )
361          {          {
362            q = M_PI_2 * atan(y);            return rc;
363          }          }
364        else if ( (nu < 21.0) && (nu > 1.0))      }
365      diff = (lg2.val) - (lg1.val);
366      q *= pow(y,nu) * exp (diff) / sqrt (M_PI);
367    
368    /*   if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */
369    /*     { */
370    /*       result->val = q; */
371    /*     } */
372    /*   else if ( tail == GSL_CDF_LOWER ) */
373    /*     { */
374    /*       result->val = 1.0 - q; */
375    /*     } */
376        
377      return q;
378    }
379    
380    double gsl_cdf_t_P ( const double x, const double nu )
381    {
382      double val = 0.0;
383      double q;
384      double absx = fabs(x);
385      
386      if ( fabs(x) < 4.0 )
387        {
388          q = gsl_cdf_t_smallx (x, nu );
389          if ( x <= 0.0 )
390          {          {
391            ckp2 = 1.0;            val = (1.0 - q) / 2.0 ;
           b = 1.0 + t*t/nu;  
             
           for ( k = inu-2; k > 1; k-=2)  
             {  
               ck = 1 + (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);  
             }  
           if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0)))  
             {  
               result->val = (1.0 - q)/2;  
             }  
           else  
             {  
               result->val = (1.0+q)/2;  
             }  
392          }          }
393        else if ( nu >= 21.0 )        else
394          {          {
395            /*            val  = (1.0 + q ) / 2.0;
396             * Use the Cornish-Fisher expansion to find          }
397             * a point u such that gsl_cdf_gauss(u) = tcdf(t).      }
398             * Approximate the t cdf with gsl_cdf_gauss.    else if ( nu >= 21.0 )
399             */      {
400            u = cornish_fisher ( t, nu );        q = t_cornish_fisher ( x, nu );
401            q = gsl_cdf_gauss ( u, GSL_CDF_UPPER );        if ( x <= 0.0)
402            if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0)))          {
403              {            val = q;
               result->val = q;  
             }  
           else  
             {  
               result->val = 1.0 - q;  
             }  
404          }          }
405        else        else
406          {          {
407            GSL_ERROR_VAL ("degrees of freedom are less than 1",            val = 1.0 - q;
                          GSL_EINVAL, nu );  
408          }          }
409      }      }
410    else    else
411      {      {
412        /*        val = normal_approx ( x, nu );
413         * Series approximation for t > 4.0. This needs to be fixed;        if ( x > 0.0 )
414         * it shouldn't subtract the result from 1.0. A better way is          {
415         * to use two different series expansions. Figuring this out            val = 1.0 - val;
416         * means rummaging through Fisher's paper in Metron, v5, 1926,          }
417         * "Expansion of Student's integral in powers of n^{-1}."      }
418         */    return val;
419    }
420    double gsl_cdf_t_Q ( const double x, const double nu )
421        y = 1/sqrt(1+t*t/nu);  {
422        num = 1.0;    double val = 0.0;
423        p = 0.0;    double q;
424        diff = 2*GSL_DBL_EPSILON;    double absx = fabs(x);
425        for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2)  
426          {    if ( fabs(x) < 4.0 )
427            diff = p;      {
428            num *= y*y*(i-1)/i;        q = gsl_cdf_t_smallx ( x, nu );
429            p += num / (nu+i);        if ( x > 0.0 )
430            diff = p - diff;          {
431          }            val = (1.0 - q) / 2.0 ;
       p += 1/nu;  
       rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1);  
       if ( rc == GSL_SUCCESS )  
         {  
           rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2);  
           if ( rc != GSL_SUCCESS )  
             {  
               return rc;  
             }  
432          }          }
433        else        else
434          {          {
435            return rc;            val  = (1.0 + q ) / 2.0;
436          }          }
437        diff = (lg2.val) - (lg1.val);      }
438        p *= pow(y,nu) * exp (diff) / sqrt (M_PI);    else if ( nu >= 21.0 )
439        {
440        if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0)))        val = t_cornish_fisher ( x, nu );
441          if ( x <= 0.0)
442          {          {
443            result->val = p;            val = 1.0 - val;
444          }          }
445        else if ( tail == GSL_CDF_LOWER )      }
446      else
447        {
448          val = normal_approx ( x, nu );
449          if ( x <= 0.0 )
450          {          {
451            result->val = 1.0 - p;            val = 1.0 - val;
452          }          }
453      }      }
454    return rc;    return val;
455  }  }

Legend:
Removed from v.1.1.2.8  
changed lines
  Added in v.1.1.2.9

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