# include # include # include # include "asa091.h" /******************************************************************************/ double ppchi2 ( double p, double v, double g, int *ifault ) /******************************************************************************/ /* Purpose: ppchi2() evaluates the percentage points of the Chi-squared PDF. Discussion Incorporates the suggested changes in AS R85 (vol.40(1), pages 233-5, 1991) which should eliminate the need for the limited range for P, though these limits have not been removed from the routine. Licensing: This code is distributed under the MIT license. Modified: 05 June 2013 Author: Original FORTRAN77 version by Donald Best, DE Roberts. C version by John Burkardt. Reference: Donald Best, DE Roberts, Algorithm AS 91: The Percentage Points of the Chi-Squared Distribution, Applied Statistics, Volume 24, Number 3, 1975, pages 385-390. Parameters: Input, double P, value of the chi-squared cumulative probability density function. 0.000002 <= P <= 0.999998. Input, double V, the parameter of the chi-squared probability density function. 0 < V. Input, double G, the value of log ( Gamma ( V / 2 ) ). Output, int *IFAULT, is nonzero if an error occurred. 0, no error. 1, P is outside the legal range. 2, V is not positive. 3, an error occurred in GAMMAD. 4, the result is probably as accurate as the machine will allow. Output, double PPCHI2, the value of the chi-squared random deviate with the property that the probability that a chi-squared random deviate with parameter V is less than or equal to PPCHI2 is P. */ { double a; double aa = 0.6931471806; double b; double c; double c1 = 0.01; double c2 = 0.222222; double c3 = 0.32; double c4 = 0.4; double c5 = 1.24; double c6 = 2.2; double c7 = 4.67; double c8 = 6.66; double c9 = 6.73; double c10 = 13.32; double c11 = 60.0; double c12 = 70.0; double c13 = 84.0; double c14 = 105.0; double c15 = 120.0; double c16 = 127.0; double c17 = 140.0; double c18 = 175.0; double c19 = 210.0; double c20 = 252.0; double c21 = 264.0; double c22 = 294.0; double c23 = 346.0; double c24 = 420.0; double c25 = 462.0; double c26 = 606.0; double c27 = 672.0; double c28 = 707.0; double c29 = 735.0; double c30 = 889.0; double c31 = 932.0; double c32 = 966.0; double c33 = 1141.0; double c34 = 1182.0; double c35 = 1278.0; double c36 = 1740.0; double c37 = 2520.0; double c38 = 5040.0; double ch; double e = 0.5E-06; int i; int if1; int maxit = 20; double pmax = 0.999998; double pmin = 0.000002; double p1; double p2; double q; double s1; double s2; double s3; double s4; double s5; double s6; double t; double value; double x; double xx; /* Test arguments and initialize. */ value = - 1.0; if ( p < pmin || pmax < p ) { *ifault = 1; return value; } if ( v <= 0.0 ) { *ifault = 2; return value; } *ifault = 0; xx = 0.5 * v; c = xx - 1.0; /* Starting approximation for small chi-squared */ if ( v < - c5 * log ( p ) ) { ch = pow ( p * xx * exp ( g + xx * aa ), 1.0 / xx ); if ( ch < e ) { value = ch; return value; } } /* Starting approximation for V less than or equal to 0.32 */ else if ( v <= c3 ) { ch = c4; a = log ( 1.0 - p ); for ( ; ; ) { q = ch; p1 = 1.0 + ch * ( c7 + ch ); p2 = ch * (c9 + ch * ( c8 + ch ) ); t = - 0.5 + (c7 + 2.0 * ch ) / p1 - ( c9 + ch * ( c10 + 3.0 * ch ) ) / p2; ch = ch - ( 1.0 - exp ( a + g + 0.5 * ch + c * aa ) * p2 / p1) / t; if ( fabs ( q / ch - 1.0 ) <= c1 ) { break; } } } else { /* Call to algorithm AS 111 - note that P has been tested above. AS 241 could be used as an alternative. */ x = ppnd ( p, ifault ); /* Starting approximation using Wilson and Hilferty estimate */ p1 = c2 / v; ch = v * pow ( x * sqrt ( p1 ) + 1.0 - p1, 3 ); /* Starting approximation for P tending to 1. */ if ( c6 * v + 6.0 < ch ) { ch = - 2.0 * ( log ( 1.0 - p ) - c * log ( 0.5 * ch ) + g ); } } /* Call to algorithm AS 239 and calculation of seven term Taylor series */ for ( i = 1; i <= maxit; i++ ) { q = ch; p1 = 0.5 * ch; p2 = p - gammad ( p1, xx, &if1 ); if ( if1 != 0 ) { *ifault = 3; return value; } t = p2 * exp ( xx * aa + g + p1 - c * log ( ch ) ); b = t / ch; a = 0.5 * t - b * c; s1 = ( c19 + a * ( c17 + a * ( c14 + a * ( c13 + a * ( c12 + c11 * a ))))) / c24; s2 = ( c24 + a * ( c29 + a * ( c32 + a * ( c33 + c35 * a )))) / c37; s3 = ( c19 + a * ( c25 + a * ( c28 + c31 * a ))) / c37; s4 = ( c20 + a * ( c27 + c34 * a) + c * ( c22 + a * ( c30 + c36 * a ))) / c38; s5 = ( c13 + c21 * a + c * ( c18 + c26 * a )) / c37; s6 = ( c15 + c * ( c23 + c16 * c )) / c38; ch = ch + t * ( 1.0 + 0.5 * t * s1 - b * c * ( s1 - b * ( s2 - b * ( s3 - b * ( s4 - b * ( s5 - b * s6 )))))); if ( e < fabs ( q / ch - 1.0 ) ) { value = ch; return value; } } *ifault = 4; value = ch; return value; }