Apply patch from Jason Stover <jstover@sdf.lonestar.org> 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 (file)
index 0000000..282522d
--- /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 (file)
index 0000000..6e462c6
--- /dev/null
+++ b/lib/gsl-extras/ChangeLog
@@ -0,0 +1,11 @@
+Thu Apr 14 10:48:17 2005  Ben Pfaff  <blp@gnu.org>
+
+       * betadistinv.c: Apply patch from Jason Stover
+       <jstover@sdf.lonestar.org> 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 0b721cac759c1158601a872d35ab1ca6cbe09761..76feb0d8d2ca8f82b797685c0ea9269672be210c 100644 (file)
--- 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.
  */
  * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8,
  * August 1968, pages 1264-1273.
  */

+

 #include <config.h>
 #include <math.h>
 #include <gsl/gsl_math.h>
 #include <config.h>
 #include <math.h>
 #include <gsl/gsl_math.h>


@@ -41,7 +42,7 @@
 #include <gsl/gsl_randist.h>
 #include "gsl-extras.h"
 
 #include <gsl/gsl_randist.h>
 #include "gsl-extras.h"
 

-#define BETAINV_INIT_ERR .01
+#define BETAINV_INIT_ERR .001

 #define BETADISTINV_N_TERMS 3
 #define BETADISTINV_MAXITER 20
 
 #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;
   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);
     {
       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;
        }
          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;
     }
 
   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;
   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, 
     {
       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;
        }
          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 e19d5f1ec61bd68ad63b6853bc29ce34642d4a69..ea138295dc13ad851790931594efd91229b18b7b 100644 (file)
--- 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 
      .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      .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.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 
      .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 2e014815dced6e2736b27bc2f437d8f9680bf0a9..ab76a257c36074c1f2848242d714a159733a4681 100644 (file)
--- 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     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 
      .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