/[pspp]/pspp/lib/gsl-extras/betadistinv.c
ViewVC logotype

Diff of /pspp/lib/gsl-extras/betadistinv.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by blp, Thu Mar 10 04:07:30 2005 UTC revision 1.3 by blp, Thu Apr 14 17:52:05 2005 UTC
# Line 32  Line 32 
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>
# Line 41  Line 42 
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    
# Line 390  gslextras_cdf_beta_Pinv ( const double p Line 391  gslextras_cdf_beta_Pinv ( const double p
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);
# Line 431  gslextras_cdf_beta_Pinv ( const double p Line 432  gslextras_cdf_beta_Pinv ( const double p
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;
# Line 581  gslextras_cdf_beta_Qinv (double q, doubl Line 595  gslextras_cdf_beta_Qinv (double q, doubl
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,
# Line 622  gslextras_cdf_beta_Qinv (double q, doubl Line 636  gslextras_cdf_beta_Qinv (double q, doubl
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    /*    /*

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26