*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/*
* 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 <gsl/gsl_errno.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
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);
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;
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,
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;
+ }
}
/*