* Given a value p in [0..1] of the lower tail area of the Chi^2 distribution * with df degrees of freedom, where ln_gamma_df_2 is ln_gamma(df/2.0), compute * the upper limit on the definite integral from [0..z] that satisfies p, * accurate to 12 decimal places. * * This implementation is based on: * * Best, D.J., D.E. Roberts (1975) Algorithm AS 91: The percentage points of * the Chi^2
| 214 | * points of the Chi^2 distribution. Applied Statistics 40(1):233-235. |
| 215 | */ |
| 216 | static inline double |
| 217 | pt_chi2(double p, double df, double ln_gamma_df_2) { |
| 218 | double e, aa, xx, c, ch, a, q, p1, p2, t, x, b, s1, s2, s3, s4, s5, s6; |
| 219 | unsigned i; |
| 220 | |
| 221 | assert(p >= 0.0 && p < 1.0); |
| 222 | assert(df > 0.0); |
| 223 | |
| 224 | e = 5.0e-7; |
| 225 | aa = 0.6931471805; |
| 226 | |
| 227 | xx = 0.5 * df; |
| 228 | c = xx - 1.0; |
| 229 | |
| 230 | if (df < -1.24 * log(p)) { |
| 231 | /* Starting approximation for small Chi^2. */ |
| 232 | ch = pow(p * xx * exp(ln_gamma_df_2 + xx * aa), 1.0 / xx); |
| 233 | if (ch - e < 0.0) { |
| 234 | return ch; |
| 235 | } |
| 236 | } else { |
| 237 | if (df > 0.32) { |
| 238 | x = pt_norm(p); |
| 239 | /* |
| 240 | * Starting approximation using Wilson and Hilferty |
| 241 | * estimate. |
| 242 | */ |
| 243 | p1 = 0.222222 / df; |
| 244 | ch = df * pow(x * sqrt(p1) + 1.0 - p1, 3.0); |
| 245 | /* Starting approximation for p tending to 1. */ |
| 246 | if (ch > 2.2 * df + 6.0) { |
| 247 | ch = -2.0 * (log(1.0 - p) - c * log(0.5 * ch) + |
| 248 | ln_gamma_df_2); |
| 249 | } |
| 250 | } else { |
| 251 | ch = 0.4; |
| 252 | a = log(1.0 - p); |
| 253 | while (true) { |
| 254 | q = ch; |
| 255 | p1 = 1.0 + ch * (4.67 + ch); |
| 256 | p2 = ch * (6.73 + ch * (6.66 + ch)); |
| 257 | t = -0.5 + (4.67 + 2.0 * ch) / p1 - (6.73 + ch |
| 258 | * (13.32 + 3.0 * ch)) / p2; |
| 259 | ch -= (1.0 - exp(a + ln_gamma_df_2 + 0.5 * ch + |
| 260 | c * aa) * p2 / p1) / t; |
| 261 | if (fabs(q / ch - 1.0) - 0.01 <= 0.0) { |
| 262 | break; |
| 263 | } |
| 264 | } |
| 265 | } |
| 266 | } |
| 267 | |
| 268 | for (i = 0; i < 20; i++) { |
| 269 | /* Calculation of seven-term Taylor series. */ |
| 270 | q = ch; |
| 271 | p1 = 0.5 * ch; |
| 272 | if (p1 < 0.0) { |
| 273 | return -1.0; |