Re: Examining bug #21838

[Brian Gough <bjg@network-theory.co.uk>]

At Fri, 15 Feb 2008 00:30:43 +0000,
Jason Stover wrote:
> The test cases have huge numbers of denominator degrees of freedom
> (\nu_2 in the manual).  As \nu_2 --> infinity, the F density tends to
> a gamma density. So it might be worthwhile to fix this by calling
> gsl_cdf_gamma_[PQ] instead of beta_inc_AXPY() in such cases.

Thanks for the suggestion, I've added a patch for that.

diff --git a/cdf/beta_inc.c b/cdf/beta_inc.c
index 8c086d0..ac1d340 100644
--- a/cdf/beta_inc.c
+++ b/cdf/beta_inc.c
@@ -20,6 +20,8 @@
 /* Author:  G. Jungman */
 /* Modified for cdfs by Brian Gough, June 2003 */
 
+#include <gsl/gsl_sf_gamma.h>
+
 static double
 beta_cont_frac (const double a, const double b, const double x,
                 const double epsabs)
@@ -100,6 +102,8 @@ beta_cont_frac (const double a, const double b, const double x,
    absolute error when A*beta_inc is added to Y. (e.g. if Y >>
    A*beta_inc then the accuracy of beta_inc can be reduced) */
 
+
+
 static double
 beta_inc_AXPY (const double A, const double Y,
                const double a, const double b, const double x)
@@ -112,6 +116,18 @@ beta_inc_AXPY (const double A, const double Y,
     {
       return A * 1 + Y;
     }
+  else if (a > 1e5 && b < 10 && x > a / (a + b))
+    {
+      /* Handle asymptotic regime, large a, small b, x > peak [AS 26.5.17] */
+      double N = a + (b - 1.0) / 2.0;
+      return A * gsl_sf_gamma_inc_Q (b, -N * log (x)) + Y;
+    }
+  else if (b > 1e5 && a < 10 && x < b / (a + b))
+    {
+      /* Handle asymptotic regime, small a, large b, x < peak [AS 26.5.17] */
+      double N = b + (a - 1.0) / 2.0;
+      return A * gsl_sf_gamma_inc_P (a, -N * log1p (-x)) + Y;
+    }
   else
     {
       double ln_beta = gsl_sf_lnbeta (a, b);
-- 
1.5.3.4