32 |
* Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8, |
* Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8, |
33 |
* August 1968, pages 1264-1273. |
* August 1968, pages 1264-1273. |
34 |
*/ |
*/ |
35 |
|
|
36 |
#include <config.h> |
#include <config.h> |
37 |
#include <math.h> |
#include <math.h> |
38 |
#include <gsl/gsl_math.h> |
#include <gsl/gsl_math.h> |
42 |
#include <gsl/gsl_randist.h> |
#include <gsl/gsl_randist.h> |
43 |
#include "gsl-extras.h" |
#include "gsl-extras.h" |
44 |
|
|
45 |
#define BETAINV_INIT_ERR .01 |
#define BETAINV_INIT_ERR .001 |
46 |
#define BETADISTINV_N_TERMS 3 |
#define BETADISTINV_N_TERMS 3 |
47 |
#define BETADISTINV_MAXITER 20 |
#define BETADISTINV_MAXITER 20 |
48 |
|
|
391 |
err = beta_result - p; |
err = beta_result - p; |
392 |
abserr = fabs(err); |
abserr = fabs(err); |
393 |
relerr = abserr / p; |
relerr = abserr / p; |
394 |
while ( relerr > BETAINV_INIT_ERR && n_iter < 100) |
while ( relerr > BETAINV_INIT_ERR) |
395 |
{ |
{ |
396 |
tmp = new_guess_P ( state, lower, upper, |
tmp = new_guess_P ( state, lower, upper, |
397 |
p, a, b); |
p, a, b); |
432 |
result = state; |
result = state; |
433 |
min_err = relerr; |
min_err = relerr; |
434 |
} |
} |
435 |
|
else |
436 |
|
{ |
437 |
|
/* |
438 |
|
* Lagrange polynomial failed to reduce the error. |
439 |
|
* This will happen with a very skewed beta density. |
440 |
|
* Undo previous steps. |
441 |
|
*/ |
442 |
|
state = result; |
443 |
|
beta_result = gsl_cdf_beta_P(state,a,b); |
444 |
|
err = p - beta_result; |
445 |
|
abserr = fabs(err); |
446 |
|
relerr = abserr / p; |
447 |
|
} |
448 |
} |
} |
449 |
|
|
450 |
n_iter = 0; |
n_iter = 0; |
595 |
err = beta_result - q; |
err = beta_result - q; |
596 |
abserr = fabs(err); |
abserr = fabs(err); |
597 |
relerr = abserr / q; |
relerr = abserr / q; |
598 |
while ( relerr > BETAINV_INIT_ERR && n_iter < 100) |
while ( relerr > BETAINV_INIT_ERR) |
599 |
{ |
{ |
600 |
n_iter++; |
n_iter++; |
601 |
tmp = new_guess_Q ( state, lower, upper, |
tmp = new_guess_Q ( state, lower, upper, |
636 |
result = state; |
result = state; |
637 |
min_err = relerr; |
min_err = relerr; |
638 |
} |
} |
639 |
|
else |
640 |
|
{ |
641 |
|
/* |
642 |
|
* Lagrange polynomial failed to reduce the error. |
643 |
|
* This will happen with a very skewed beta density. |
644 |
|
* Undo previous steps. |
645 |
|
*/ |
646 |
|
state = result; |
647 |
|
beta_result = gsl_cdf_beta_P(state,a,b); |
648 |
|
err = q - beta_result; |
649 |
|
abserr = fabs(err); |
650 |
|
relerr = abserr / q; |
651 |
|
} |
652 |
} |
} |
653 |
|
|
654 |
/* |
/* |