From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 8731 invoked by alias); 20 Feb 2008 18:04:20 -0000 Received: (qmail 8723 invoked by uid 22791); 20 Feb 2008 18:04:19 -0000 X-Spam-Check-By: sourceware.org Received: from mail.network-theory.co.uk (HELO mail.network-theory.co.uk) (66.199.228.187) by sourceware.org (qpsmtpd/0.31) with ESMTP; Wed, 20 Feb 2008 18:03:58 +0000 Date: Wed, 20 Feb 2008 18:04:00 -0000 Message-ID: <871w77jyie.wl%bjg@network-theory.co.uk> From: Brian Gough To: gsl-discuss@sourceware.org Subject: Re: Examining bug #21838 In-Reply-To: <20080215003043.GA7666@SDF.LONESTAR.ORG> References: <47B4C827.1040004@googlemail.com> <20080215003043.GA7666@SDF.LONESTAR.ORG> User-Agent: Wanderlust/2.14.0 (Africa) Emacs/22.1 Mule/5.0 (SAKAKI) MIME-Version: 1.0 (generated by SEMI 1.14.6 - "Maruoka") Content-Type: text/plain; charset=US-ASCII X-Message-Mac: c8a3b6da33a5eaf0fddff19e23f3a5c4 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2008-q1/txt/msg00035.txt.bz2 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 + 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