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
