From: Ben Pfaff Date: Thu, 14 Apr 2005 17:52:05 +0000 (+0000) Subject: Apply patch from Jason Stover to obtain X-Git-Tag: v0.4.0~121 X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;ds=sidebyside;h=084280e4bfd6b94c26320e72553a9a38201610a3;p=pspp-builds.git Apply patch from Jason Stover to obtain better initial approximation for the inverse beta distribution and tidy up a bit. Update test results for affected values in the beta and F distributions. --- diff --git a/lib/gsl-extras/.cvsignore b/lib/gsl-extras/.cvsignore new file mode 100644 index 00000000..282522db --- /dev/null +++ b/lib/gsl-extras/.cvsignore @@ -0,0 +1,2 @@ +Makefile +Makefile.in diff --git a/lib/gsl-extras/ChangeLog b/lib/gsl-extras/ChangeLog new file mode 100644 index 00000000..6e462c61 --- /dev/null +++ b/lib/gsl-extras/ChangeLog @@ -0,0 +1,11 @@ +Thu Apr 14 10:48:17 2005 Ben Pfaff + + * betadistinv.c: Apply patch from Jason Stover + to obtain better initial approximation + and tidy up a bit. + +---------------------------------------------------------------------- +Local Variables: +mode: change-log +version-control: never +End: diff --git a/lib/gsl-extras/betadistinv.c b/lib/gsl-extras/betadistinv.c index 0b721cac..76feb0d8 100644 --- a/lib/gsl-extras/betadistinv.c +++ b/lib/gsl-extras/betadistinv.c @@ -32,6 +32,7 @@ * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8, * August 1968, pages 1264-1273. */ + #include #include #include @@ -41,7 +42,7 @@ #include #include "gsl-extras.h" -#define BETAINV_INIT_ERR .01 +#define BETAINV_INIT_ERR .001 #define BETADISTINV_N_TERMS 3 #define BETADISTINV_MAXITER 20 @@ -390,7 +391,7 @@ gslextras_cdf_beta_Pinv ( const double p, const double a, const double b) err = beta_result - p; abserr = fabs(err); relerr = abserr / p; - while ( relerr > BETAINV_INIT_ERR && n_iter < 100) + while ( relerr > BETAINV_INIT_ERR) { tmp = new_guess_P ( state, lower, upper, p, a, b); @@ -431,6 +432,19 @@ gslextras_cdf_beta_Pinv ( const double p, const double a, const double b) result = state; min_err = relerr; } + else + { + /* + * Lagrange polynomial failed to reduce the error. + * This will happen with a very skewed beta density. + * Undo previous steps. + */ + state = result; + beta_result = gsl_cdf_beta_P(state,a,b); + err = p - beta_result; + abserr = fabs(err); + relerr = abserr / p; + } } n_iter = 0; @@ -581,7 +595,7 @@ gslextras_cdf_beta_Qinv (double q, double a, double b) err = beta_result - q; abserr = fabs(err); relerr = abserr / q; - while ( relerr > BETAINV_INIT_ERR && n_iter < 100) + while ( relerr > BETAINV_INIT_ERR) { n_iter++; tmp = new_guess_Q ( state, lower, upper, @@ -622,6 +636,19 @@ gslextras_cdf_beta_Qinv (double q, double a, double b) result = state; min_err = relerr; } + else + { + /* + * Lagrange polynomial failed to reduce the error. + * This will happen with a very skewed beta density. + * Undo previous steps. + */ + state = result; + beta_result = gsl_cdf_beta_P(state,a,b); + err = q - beta_result; + abserr = fabs(err); + relerr = abserr / q; + } } /* diff --git a/tests/expressions/randist/beta.out b/tests/expressions/randist/beta.out index e19d5f1e..ea138295 100644 --- a/tests/expressions/randist/beta.out +++ b/tests/expressions/randist/beta.out @@ -208,21 +208,21 @@ .90 1.50 2.00 .76 .90 .79 .90 1.50 2.50 .68 .90 .75 .90 1.50 3.00 .62 .90 .74 - .99 .75 .25 1.00 .98 20860.76 + .99 .75 .25 1.00 .99 164255.7 .99 .75 .50 1.00 .99 34.83 .99 .75 1.00 .99 .99 .75 .99 .75 1.50 .94 .99 .26 .99 .75 2.00 .88 .99 .17 .99 .75 2.50 .81 .99 .13 .99 .75 3.00 .75 .99 .12 - .99 1.00 .25 1.00 .98 46067.00 + .99 1.00 .25 1.00 .99 250000.0 .99 1.00 .50 1.00 .99 50.00 .99 1.00 1.00 .99 .99 1.00 .99 1.00 1.50 .95 .99 .32 .99 1.00 2.00 .90 .99 .20 .99 1.00 2.50 .84 .99 .16 .99 1.00 3.00 .78 .99 .14 - .99 1.50 .25 1.00 .98 67836.74 + .99 1.50 .25 1.00 .99 428406.6 .99 1.50 .50 1.00 .99 81.05 .99 1.50 1.00 .99 .99 1.49 .99 1.50 1.50 .97 .99 .45 diff --git a/tests/expressions/randist/f.out b/tests/expressions/randist/f.out index 2e014815..ab76a257 100644 --- a/tests/expressions/randist/f.out +++ b/tests/expressions/randist/f.out @@ -169,7 +169,7 @@ .99 5.00 1.00 5763.65 .99 .00 .99 5.00 2.00 99.30 .99 .00 .99 5.00 5.00 10.97 .99 .00 - .99 5.00 10.00 10.00 1.00 .00 + .99 5.00 10.00 5.64 .99 .01 .99 10.00 1.00 6055.85 .99 .00 .99 10.00 2.00 99.40 .99 .00 .99 10.00 5.00 10.05 .99 .00