| /* |
| * Copyright (C) 2000-2015 The R Core Team |
| * |
| * Algorithm AS 226 Appl. Statist. (1987) Vol. 36, No. 2 |
| * by Russell V. Lenth |
| * Incorporates modification AS R84 from AS Vol. 39, pp311-2, 1990 |
| * and AS R95 from AS Vol. 44, pp551-2, 1995 |
| * by H. Frick and Min Long Lam. |
| * original (C) Royal Statistical Society 1987, 1990, 1995 |
| * |
| * Returns the cumulative probability of x for the non-central |
| * beta distribution with parameters a, b and non-centrality ncp. |
| * |
| * Auxiliary routines required: |
| * lgamma - log-gamma function |
| * pbeta - incomplete-beta function {nowadays: pbeta_raw() -> bratio()} |
| */ |
| |
| #include "nmath.h" |
| #include "dpq.h" |
| |
| LDOUBLE attribute_hidden |
| pnbeta_raw(double x, double o_x, double a, double b, double ncp) |
| { |
| /* o_x == 1 - x but maybe more accurate */ |
| |
| /* change errmax and itrmax if desired; |
| * original (AS 226, R84) had (errmax; itrmax) = (1e-6; 100) */ |
| const static double errmax = 1.0e-9; |
| const int itrmax = 10000; /* 100 is not enough for pf(ncp=200) |
| see PR#11277 */ |
| |
| double a0, lbeta, c, errbd, x0, temp, tmp_c; |
| int ierr; |
| |
| LDOUBLE ans, ax, gx, q, sumq; |
| |
| if (ncp < 0. || a <= 0. || b <= 0.) ML_ERR_return_NAN; |
| |
| if(x < 0. || o_x > 1. || (x == 0. && o_x == 1.)) return 0.; |
| if(x > 1. || o_x < 0. || (x == 1. && o_x == 0.)) return 1.; |
| |
| c = ncp / 2.; |
| |
| /* initialize the series */ |
| |
| x0 = floor(fmax2(c - 7. * sqrt(c), 0.)); |
| a0 = a + x0; |
| lbeta = lgammafn(a0) + lgammafn(b) - lgammafn(a0 + b); |
| /* temp = pbeta_raw(x, a0, b, TRUE, FALSE), but using (x, o_x): */ |
| bratio(a0, b, x, o_x, &temp, &tmp_c, &ierr, FALSE); |
| |
| gx = exp(a0 * log(x) + b * (x < .5 ? log1p(-x) : log(o_x)) |
| - lbeta - log(a0)); |
| if (a0 > a) |
| q = exp(-c + x0 * log(c) - lgammafn(x0 + 1.)); |
| else |
| q = exp(-c); |
| |
| sumq = 1. - q; |
| ans = ax = q * temp; |
| |
| /* recurse over subsequent terms until convergence is achieved */ |
| double j = floor(x0); // x0 could be billions, and is in package EnvStats |
| do { |
| j++; |
| temp -= (double) gx; |
| gx *= x * (a + b + j - 1.) / (a + j); |
| q *= c / j; |
| sumq -= q; |
| ax = temp * q; |
| ans += ax; |
| errbd = (double)((temp - gx) * sumq); |
| } |
| while (errbd > errmax && j < itrmax + x0); |
| |
| if (errbd > errmax) |
| ML_ERROR(ME_PRECISION, "pnbeta"); |
| if (j >= itrmax + x0) |
| ML_ERROR(ME_NOCONV, "pnbeta"); |
| |
| return ans; |
| } |
| |
| double attribute_hidden |
| pnbeta2(double x, double o_x, double a, double b, double ncp, |
| /* o_x == 1 - x but maybe more accurate */ |
| int lower_tail, int log_p) |
| { |
| LDOUBLE ans = pnbeta_raw(x, o_x, a,b, ncp); |
| |
| /* return R_DT_val(ans), but we want to warn about cancellation here */ |
| if (lower_tail) |
| #ifdef HAVE_LONG_DOUBLE |
| return (double) (log_p ? logl(ans) : ans); |
| #else |
| return log_p ? log(ans) : ans; |
| #endif |
| else { |
| if (ans > 1. - 1e-10) ML_ERROR(ME_PRECISION, "pnbeta"); |
| if (ans > 1.0) ans = 1.0; /* Precaution */ |
| #if defined(HAVE_LONG_DOUBLE) && defined(HAVE_LOG1PL) |
| return (double) (log_p ? log1pl(-ans) : (1. - ans)); |
| #else |
| /* include standalone case */ |
| return (double) (log_p ? log1p((double)-ans) : (1. - ans)); |
| #endif |
| } |
| } |
| |
| double pnbeta(double x, double a, double b, double ncp, |
| int lower_tail, int log_p) |
| { |
| #ifdef IEEE_754 |
| if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(ncp)) |
| return x + a + b + ncp; |
| #endif |
| |
| R_P_bounds_01(x, 0., 1.); |
| return pnbeta2(x, 1-x, a, b, ncp, lower_tail, log_p); |
| } |