/* * R : A Computer Language for Statistical Data Analysis * Copyright (C) 1998--2022 The R Core Team * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka * based on code (C) 1979 and later Royal Statistical Society * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, a copy is available at * https://www.R-project.org/Licenses/ * * Reference: * Cran, G. W., K. J. Martin and G. E. Thomas (1977). * Remark AS R19 and Algorithm AS 109, * Applied Statistics, 26(1), 111-114. * Remark AS R83 (v.39, 309-310) and the correction (v.40(1) p.236) * have been incorporated in this version. */ #include "nmath.h" #include "dpq.h" /**----------- DEBUGGING ------------- * * make CFLAGS='-DDEBUG_qbeta ...' *MM (w/ Debug, w/o Optimization): if not manually adding -dDEBUG_qbeta , use (cd `R-devel-qbeta-dbg RHOME`/src/nmath ; make qbeta.o; cd ../..; make R) */ #ifdef DEBUG_qbeta # define R_ifDEBUG_printf(...) REprintf(__VA_ARGS__) #else # define R_ifDEBUG_printf(...) #endif #define USE_LOG_X_CUTOFF -5. // --- based on some testing; had = -10 #define n_NEWTON_FREE 4 // --- based on some testing; had = 10 #define MLOGICAL_NA -1 // an "NA_LOGICAL" substitute for Mathlib {only used here, for now} //attribute_hidden static void qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, int swap_01, double log_q_cut, int n_N, double* qb); double qbeta(double alpha, double p, double q, int lower_tail, int log_p) { /* test for admissibility of parameters */ #ifdef IEEE_754 if (ISNAN(p) || ISNAN(q) || ISNAN(alpha)) return p + q + alpha; #endif if(p < 0. || q < 0.) ML_WARN_return_NAN; // allowing p==0 and q==0 <==> treat as one- or two-point mass double qbet[2];// = { qbeta(), 1 - qbeta() } qbeta_raw(alpha, p, q, lower_tail, log_p, // swap_01 , log_q_cut , n_N MLOGICAL_NA, USE_LOG_X_CUTOFF, n_NEWTON_FREE, qbet); return qbet[0]; } static const double #ifdef IEEE_754 // CARE: assumes subnormal numbers, i.e., no underflow at DBL_MIN: DBL_very_MIN = DBL_MIN / 4., DBL_log_v_MIN = M_LN2*(DBL_MIN_EXP - 2), // Too extreme: inaccuracy in pbeta(); e.g for qbeta(0.95, 1e-9, 20): // -> in pbeta() --> bgrat(..... b*z == 0 underflow, hence inaccurate pbeta() /* DBL_very_MIN = 0x0.0000001p-1022, // = 2^-1050 = 2^(-1022 - 28) */ /* DBL_log_v_MIN = -1050. * M_LN2, // = log(DBL_very_MIN) ~= -727.8045 */ // the most extreme -- not ok, as pbeta() then behaves strangely, // e.g., for qbeta(0.95, 1e-8, 20): /* DBL_very_MIN = 0x0.0000000000001p-1022, // = 2^-1074 = 2^(-1022 -52) */ /* DBL_log_v_MIN = -1074. * M_LN2, // = log(DBL_very_MIN) */ DBL_1__eps = 0x1.fffffffffffffp-1; // = 1 - 2^-53 #else // untested : DBL_1__eps = 1 - DBL_EPSILON; // or rather (1 - DBL_EPSILON/2) (??) #endif /* set the exponent of acu to -2r-2 for r digits of accuracy */ /*---- NEW ---- -- still fails for p = 1e11, q=.5*/ #define fpu 3e-308 /* acu_min: Minimal value for accuracy 'acu' which will depend on (a,p); acu_min >= fpu ! */ #define acu_min 1e-300 #define p_lo fpu #define p_hi 1-2.22e-16 #define const1 2.30753 #define const2 0.27061 #define const3 0.99229 #define const4 0.04481 // Returns both qbeta() and its "mirror" 1-qbeta(). Useful notably when qbeta() ~= 1 attribute_hidden void qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, int swap_01, // {TRUE, NA, FALSE}: if NA, algorithm decides swap_tail double log_q_cut, /* if == Inf: return log(qbeta(..)); otherwise, if finite: the bound for switching to log(x)-scale; see use_log_x */ int n_N, // number of "unconstrained" Newton steps before switching to constrained double *qb) // = qb[0:1] = { qbeta(), 1 - qbeta() } { Rboolean swap_choose = (swap_01 == MLOGICAL_NA), swap_tail, log_, give_log_q = (log_q_cut == ML_POSINF), use_log_x = give_log_q, // or u < log_q_cut below warned = FALSE, add_N_step = TRUE; int i_pb, i_inn; double a, la, logbeta, g, h, pp, p_, qq, r, s, t, w, y = -1.; volatile double u, xinbta; // Assuming p >= 0, q >= 0 here ... // Deal with boundary cases here: if(alpha == R_DT_0) { #define return_q_0 \ if(give_log_q) { qb[0] = ML_NEGINF; qb[1] = 0; } \ else { qb[0] = 0; qb[1] = 1; } \ return return_q_0; } if(alpha == R_DT_1) { #define return_q_1 \ if(give_log_q) { qb[0] = 0; qb[1] = ML_NEGINF; } \ else { qb[0] = 1; qb[1] = 0; } \ return return_q_1; } // check alpha {*before* transformation which may lose all accuracy}: if((log_p && alpha > 0) || (!log_p && (alpha < 0 || alpha > 1))) { // alpha is outside R_ifDEBUG_printf("qbeta(alpha=%g, %g, %g, .., log_p=%d): %s%s\n", alpha, p,q, log_p, "alpha not in ", log_p ? "[-Inf, 0]" : "[0,1]"); // ML_WARN_return_NAN : ML_WARNING(ME_DOMAIN, ""); qb[0] = qb[1] = ML_NAN; return; } // p==0, q==0, p = Inf, q = Inf <==> treat as one- or two-point mass if(p == 0 || q == 0 || !R_FINITE(p) || !R_FINITE(q)) { // We know 0 < T(alpha) < 1 : pbeta() is constant and trivial in {0, 1/2, 1} R_ifDEBUG_printf( "qbeta(%g, %g, %g, lower_t=%d, log_p=%d): (p,q)-boundary: trivial\n", alpha, p,q, lower_tail, log_p); if(p == 0 && q == 0) { // point mass 1/2 at each of {0,1} : if(alpha < R_D_half) { return_q_0; } if(alpha > R_D_half) { return_q_1; } // else: alpha == "1/2" #define return_q_half \ if(give_log_q) qb[0] = qb[1] = -M_LN2; \ else qb[0] = qb[1] = 0.5; \ return return_q_half; } else if (p == 0 || p/q == 0) { // point mass 1 at 0 - "flipped around" return_q_0; } else if (q == 0 || q/p == 0) { // point mass 1 at 0 - "flipped around" return_q_1; } // else: p = q = Inf : point mass 1 at 1/2 return_q_half; } /* initialize */ p_ = R_DT_qIv(alpha);/* lower_tail prob (in any case) */ // Conceptually, 0 < p_ < 1 (but can be 0 or 1 because of cancellation!) logbeta = lbeta(p, q); swap_tail = (swap_choose) ? (p_ > 0.5) : swap_01; R_ifDEBUG_printf( "qbeta(%g, %g, %g, lower_t=%d, log_p=%d, swap_01=%d, log_q_cut=%g, n_N=%d):%s\n" " swap_tail=%s :", alpha, p,q, lower_tail, log_p, swap_01, log_q_cut, n_N, (log_p && (p_ == 0. || p_ == 1.)) ? (p_==0.?" p_=0":" p_=1") : "", (swap_tail ? "TRUE": "F")); int n_maybe_swaps = 0; maybe_swap: // change tail; default (swap_01 = NA): afterwards 0 < a <= 1/2 if(swap_tail) { /* change tail, swap copies of {p,q}: p <-> q :*/ a = R_DT_CIv(alpha); // = 1 - p_ , is < 1/2 if(swap_choose) /* la := log(a), but without numerical cancellation: */ la = R_DT_Clog(alpha); pp = q; qq = p; } else { a = p_; la = R_DT_log(alpha); pp = p; qq = q; } n_maybe_swaps++; /* calculate the initial approximation */ /* Desired accuracy for Newton iterations (below) should depend on (a,p) * This is from Remark .. on AS 109, adapted. * However, it's not clear if this is "optimal" for IEEE double prec. * acu = fmax2(acu_min, pow(10., -25. - 5./(pp * pp) - 1./(a * a))); * NEW: 'acu' accuracy NOT for squared adjustment, but simple; * ---- i.e., "new acu" = sqrt(old acu) */ double acu = fmax2(acu_min, pow(10., -13. - 2.5/(pp * pp) - 0.5/(a * a))); // try to catch "extreme left tail" early double tx, u0 = (la + log(pp) + logbeta) / pp, // = log(x_0) rp = pp*(1.-qq)/(pp+1.); static const double log_eps_c = M_LN2 * (1. - DBL_MANT_DIG);// = log(DBL_EPSILON) = -36.04.. t = 0.2; // FIXME: Factor 0.2 is a bit arbitrary; '1' is clearly much too much. Rboolean u0_maybe = (M_LN2 * DBL_MIN_EXP < u0 && u0 < -0.01); /* 1. cannot allow exp(u0) = 0 ==> exp(u1) = exp(u0) = 0 * 2. must: u0 < 0, but too close to 0 <==> x = exp(u0) = 0.99.. */ R_ifDEBUG_printf( " n_maybe_swaps=%d, la=%#8g, u0=%#8g -> u0_maybe=%s (bnd: %g (%g)); ", n_maybe_swaps, la, u0, (u0_maybe ? "TRUE" : "F"), u0_maybe ? (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2. : -0., u0_maybe ? t*log_eps_c - log(fabs(rp)) : -0. ); double u_n = 1.; // to be log(xinbta) <==> xinbta = exp(u_n). 1 is impossible if(u0_maybe && // qq <= 2 && // <--- "arbitrary" // u0 < t*log_eps_c - log(fabs(rp)) && u0 < (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2.) { // TODO: maybe jump here from below, when initial u "fails" ? // L_tail_u: // MM's one-step correction (cheaper than 1 Newton!) rp = rp*exp(u0);// = rp*x0 if(rp > -1.) { u = u0 - log1p(rp)/pp; R_ifDEBUG_printf("u1-u0=%9.3g --> choosing u = u1 = %g\n", u-u0, u); } else { u = u0; R_ifDEBUG_printf("cannot cheaply improve u0\n"); } tx = xinbta = exp(u); use_log_x = TRUE; // or (u < log_q_cut) ?? goto L_Newton; } // y := y_\alpha in AS 64 := Hastings(1955) approximation of qnorm(1 - a) : r = sqrt(-2 * la); y = r - (const1 + const2 * r) / (1. + (const3 + const4 * r) * r); R_ifDEBUG_printf("y_a(AS64)=%g\n ", y); if (pp > 1 && qq > 1) { // use Carter(1947), see AS 109, remark '5.' r = (y * y - 3.) / 6.; s = 1. / (pp + pp - 1.); t = 1. / (qq + qq - 1.); h = 2. / (s + t); w = y * sqrt(h + r) / h - (t - s) * (r + 5. / 6. - 2. / (3. * h)); R_ifDEBUG_printf("p,q > 1 => w=%g <= 300 ? ", w); if(w > 300) { // exp(w+w) is huge or overflows t = w+w + log(qq) - log(pp); // = argument of log1pexp(.) u = // log(xinbta) = - log1p(qq/pp * exp(w+w)) = -log(1 + exp(t)) (t <= 18) ? -log1p(exp(t)) : -t - exp(-t); xinbta = exp(u); } else { xinbta = pp / (pp + qq * exp(w + w)); u = // log(xinbta) - log1p(qq/pp * exp(w+w)); } } else { // use the original AS 64 proposal, Scheffé-Tukey (1944) and Wilson-Hilferty r = qq + qq; /* A slightly more stable version of t := \chi^2_{alpha} of AS 64 * t = 1. / (9. * qq); t = r * R_pow_di(1. - t + y * sqrt(t), 3); */ t = 1. / (3. * sqrt(qq)); // = sqrt(t) of formula above t = r * R_pow_di(1. + t*(-t + y), 3);// = \chi^2_{alpha} of AS 64 s = 4. * pp + r - 2.;// 4p + 2q - 2 = numerator of new t' = s / t = s / chi^2 R_ifDEBUG_printf("min(p,q) <= 1: t=%g, s=%g", t, s); if (t == 0 || (t < 0. && s >= t)) { // cannot use chisq approx // x0 = 1 - { (1-a)*q*B(p,q) } ^{1/q} {AS 65} // xinbta = 1. - exp((log(1-a)+ log(qq) + logbeta) / qq); double l1ma = /* := log(1-a), directly from alpha (as 'la' above); though only seen very small improvements */ swap_tail ? R_DT_log(alpha) : R_DT_Clog(alpha); R_ifDEBUG_printf( " t <= 0 => AS65: log1p(-a)=%7g, better l1ma=%.15g, relD=%g\n", log1p(-a), l1ma, log1p(-a)/l1ma - 1); double xx = (l1ma + log(qq) + logbeta) / qq; R_ifDEBUG_printf(" xx = (l1ma + log(qq) + logbeta) / qq = %.10g; ", xx); if(xx <= 0.) { xinbta = -expm1(xx); u = R_Log1_Exp (xx);// = log(xinbta) = log(1 - exp(...A...)) R_ifDEBUG_printf(" xx <= 0 useful ==> u, xinbta\n"); } else { // xx > 0 ==> 1 - e^xx < 0 .. is nonsense R_ifDEBUG_printf(" xx > 0: xinbta:= 1-e^xx < 0 'nonsense'; "); // Try MM's one-step correction (or else u0) double r_ = rp*exp(u0); if(r_ > -1.) { u = u0 - log1p(r_)/pp; R_ifDEBUG_printf("u:= u0 - log1p(r_)/pp = %g\n"); } else { u = u0; R_ifDEBUG_printf("u:= u0 = %g\n"); } xinbta = exp(u); } } else { t = s / t; R_ifDEBUG_printf(" t > 0 or s < t < 0: new t := s/t = %g ( > 1 ?)\n", t); if (t <= 1.) { // cannot use chisq, either u = u0; xinbta = exp(u); } else { // (1+x0)/(1-x0) = t, solved for x0 : xinbta = 1. - 2. / (t + 1.); u = log1p(-2. / (t + 1.)); } } } // Problem: If initial u is completely wrong, we make a wrong decision here if(swap_choose && // vvvv/ why -exp(*)? u on log-x scale! Swapping can be very good, but needs smart if(..) (( swap_tail && u >= -exp( log_q_cut)) || // ==> "swap back" (!swap_tail && u >= -exp(4*log_q_cut) && pp / qq < 1000.) // ==> "swap now" )) { if(swap_tail) R_ifDEBUG_printf("\"swap back\" as u = %g >= -exp(log_q_cut);", u); else R_ifDEBUG_printf("\"swap now\" as u = %g >= -exp(4*log_q_cut) && pp/qq = %g < 1000;", u, pp/qq); // reverse swap (and typically use_log_x) swap_tail = !swap_tail; if(swap_tail) { // "swap now" (much less easily) a = R_DT_CIv(alpha); // needed ? la = R_DT_Clog(alpha); pp = q; qq = p; } else { // "swap back" : a = p_; la = R_DT_log(alpha); pp = p; qq = q; } R_ifDEBUG_printf(" ==> la = %g\n", la); // we could redo computations above, but this should be stable u = R_Log1_Exp(u); xinbta = exp(u); /* Careful: "swap now" should not fail if 1) the above initial xinbta is "completely wrong" 2) The correction step can go outside (u_n > 0 ==> e^u > 1 is illegal) e.g., for qbeta(0.2066, 0.143891, 0.05) */ } if(!use_log_x) use_log_x = (u < log_q_cut);// <==> xinbta = e^u < exp(log_q_cut) Rboolean bad_u = !R_FINITE(u), bad_init = bad_u || xinbta > p_hi; R_ifDEBUG_printf(" -> u = %g, e^u = xinbta = %.16g, (Newton acu=%g%s%s%s)\n", u, xinbta, acu, (bad_u ? ", ** bad u **" : ""), ((bad_init && !bad_u) ? ", ** bad_init **" : ""), (use_log_x ? ", on u = LOG(x) SCALE" : "")); tx = xinbta; // keeping "original initial x" (for now) if(bad_u || u < log_q_cut) { /* e.g. qbeta(0.21, .001, 0.05) try "left border" quickly, i.e., try at smallest positive number: */ w = pbeta_raw(DBL_very_MIN, pp, qq, TRUE, log_p); if(w > (log_p ? la : a)) { R_ifDEBUG_printf( " quantile is left of %g: boundary \"convergence\"\n", DBL_very_MIN); if(log_p || fabs(w - a) < fabs(0 - a)) { // DBL_very_MIN is better than 0 tx = DBL_very_MIN; u_n = DBL_log_v_MIN;// = log(DBL_very_MIN) } else { tx = 0.; u_n = ML_NEGINF; } use_log_x = log_p; add_N_step = FALSE; goto L_return; } else { R_ifDEBUG_printf(" pbeta(%g, %g, %g, T, log) = %g <= %g (= %s) --> continuing\n", DBL_very_MIN, pp,qq, w, (log_p ? la : a), (log_p ? "la" : "a")); if(u < DBL_log_v_MIN) { u = DBL_log_v_MIN;// = log(DBL_very_MIN) xinbta = DBL_very_MIN; } } } /* Sometimes the approximation is negative (and == 0 is also not "ok") */ if (bad_init && !(use_log_x && tx > 0)) { if(u == ML_NEGINF) { R_ifDEBUG_printf(" u = -Inf;"); u = M_LN2 * DBL_MIN_EXP; xinbta = DBL_MIN; } else { R_ifDEBUG_printf(" bad_init: u=%g, xinbta=%g;", u,xinbta); xinbta = (xinbta > 1.1) // i.e. "way off" ? 0.5 // otherwise, keep the respective boundary: : ((xinbta < p_lo) ? exp(u) : p_hi); if(bad_u) u = log(xinbta); // otherwise: not changing "potentially better" u than the above } R_ifDEBUG_printf(" -> (partly)new u=%g, xinbta=%g\n", u,xinbta); } L_Newton: /* -------------------------------------------------------------------- * Solve for x by a modified Newton-Raphson method, using pbeta_raw() */ r = 1 - pp; t = 1 - qq; double wprev = 0., prev = 1., adj = 1.; // -Wall if(use_log_x) { // find log(xinbta) -- work in u := log(x) scale // if(bad_init && tx > 0) xinbta = tx;// may have been better for (i_pb=0; i_pb < 1000; i_pb++) { // using log_p == TRUE unconditionally here /* FIXME: if exp(u) = xinbta underflows to 0, * want different formula pbeta_log(u, ..) */ y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, TRUE); /* w := Newton step size for L(u) = log F(e^u) =!= 0; u := log(x) * = (L(.) - la) / L'(.); L'(u)= (F'(e^u) * e^u ) / F(e^u) * = (L(.) - la)*F(.) / {F'(e^u) * e^u } = * = (L(.) - la) * e^L(.) * e^{-log F'(e^u) - u} * = ( y - la) * e^{ y - u -log F'(e^u)} and -log F'(x)= -log f(x) = - -logbeta + (1-p) log(x) + (1-q) log(1-x) = logbeta + (1-p) u + (1-q) log(1-e^u) */ w = (y == ML_NEGINF) // y = -Inf well possible: we are on log scale! ? 0. : (y - la) * exp(y - u + logbeta + r * u + t * R_Log1_Exp(u)); R_ifDEBUG_printf("N(i=%2d): u=%#20.16g, lnpb(e^u)=%#15.11g, w=%#12.7g", i_pb, u, y, w); if(!R_FINITE(w)) { // what we should do is ==> go back, get better starting value R_ifDEBUG_printf(" -- not finite --> %s\n", (n_maybe_swaps <= 1) ? "goto maybe_swap" : "give up"); if(n_maybe_swaps <= 1) goto maybe_swap; /* else was 'break;' ... but rather give up returning NaN directly as in "normal scale" Newton */ ML_WARNING(ME_DOMAIN, ""); qb[0] = qb[1] = ML_NAN; return; } if (i_pb >= n_N && w * wprev <= 0.) prev = fmax2(fabs(adj),fpu); R_ifDEBUG_printf(", %s prev=%g,", (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev); g = 1; for (i_inn=0; i_inn < 1000; i_inn++) { adj = g * w; // safe guard (here, from the very beginning) if (fabs(adj) < prev) { u_n = u - adj; // u_{n+1} = u_n - g*w if (u_n <= 0.) { // <==> 0 < xinbta := e^u <= 1 if (prev <= acu || fabs(w) <= acu) { R_ifDEBUG_printf( " it{in}=%d, -adj=%g, %s <= acu ==> convergence\n", i_inn, -adj, (prev <= acu) ? "prev" : "|w|"); goto L_converged; } // if (u_n != ML_NEGINF && u_n != 1) break; } } g /= 3; } // (cancellation in (u_n -u) => may differ from adj: double D = fmin2(fabs(adj), fabs(u_n - u)); /* R_ifDEBUG_printf(" delta(u)=%g\n", u_n - u); */ R_ifDEBUG_printf(" it{in}=%d, d.(u)=%6.3g, D/|.|=%.3g\n", i_inn, u_n - u, D/fabs(u_n + u)); if (D <= 4e-16 * fabs(u_n + u)) goto L_converged; u = u_n; xinbta = exp(u); wprev = w; } // for(i ) } else { // "normal scale" Newton for (i_pb=0; i_pb < 1000; i_pb++) { y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, log_p); // delta{y} : d_y = y - (log_p ? la : a); /* w := Newton step size (F(.) - a) / F'(.) or, * -- log: (lF - la) / (F' / F) = exp(lF) * (lF - la) / F' */ w = log_p ? (y - la) * exp(y + logbeta + r * log(xinbta) + t * log1p(-xinbta)) : (y - a) * exp( logbeta + r * log(xinbta) + t * log1p(-xinbta)); if(!R_FINITE(w)) { // what we should do is ==> go back, get better starting value R_ifDEBUG_printf( "N(i=%2d): x0=%#19.15g, pb(x0)=%#15.11g, w=%#12.7g -- is not finite --> %s\n", i_pb, xinbta, y, w, (n_maybe_swaps <= 2) ? "goto maybe_swap" : "give up"); if(n_maybe_swaps <= 2) { if(!log_p && n_maybe_swaps == 2) use_log_x = TRUE; // try now if(!log_p || n_maybe_swaps <= 1) goto maybe_swap ; } /* else was 'break;' ... but rather give up returning NaN directly as in "normal scale" Newton */ ML_WARNING(ME_DOMAIN, ""); qb[0] = qb[1] = ML_NAN; return; } if (i_pb >= n_N && w * wprev <= 0.) prev = fmax2(fabs(adj),fpu); R_ifDEBUG_printf( "N(i=%2d): x0=%#19.15g, pb(x0)=%#15.11g, w=%#12.7g, %s prev=%g,", i_pb, xinbta, y, w, (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev); g = 1; for (i_inn=0; i_inn < 1000;i_inn++) { adj = g * w; // take full Newton steps at the beginning; only then safe guard: if (i_pb < n_N || fabs(adj) < prev) { tx = xinbta - adj; // x_{n+1} = x_n - g*w if (0. <= tx && tx <= 1.) { if (prev <= acu || fabs(w) <= acu) { R_ifDEBUG_printf(" it{in}=%d, delta(x)=%g, %s <= acu ==> convergence\n", i_inn, -adj, (prev <= acu) ? "prev" : "|w|"); goto L_converged; } if (tx != 0. && tx != 1) break; } } g /= 3; } R_ifDEBUG_printf(" it{in}=%d, delta(x)=%g\n", i_inn, tx - xinbta); if (fabs(tx - xinbta) <= 4e-16 * (tx + xinbta)) // "<=" : (.) == 0 goto L_converged; xinbta = tx; if(tx == 0) // "we have lost" break; wprev = w; } // for( i_pb ..) } // end{else : normal scale Newton} /*-- NOT converged: Iteration count --*/ warned = TRUE; ML_WARNING(ME_PRECISION, "qbeta"); L_converged: log_ = log_p || use_log_x; R_ifDEBUG_printf(" %s: Final delta(y) = %g%s\n", warned ? "_NO_ convergence" : "converged", y - (log_ ? la : a), (log_ ? " (log_)" : "")); if((log_ && y == ML_NEGINF) || (!log_ && y == 0)) { // stuck at left, try if smallest positive number is "better" w = pbeta_raw(DBL_very_MIN, pp, qq, TRUE, log_); if(log_ || fabs(w - a) <= fabs(y - a)) { tx = DBL_very_MIN; u_n = DBL_log_v_MIN;// = log(DBL_very_MIN) } add_N_step = FALSE; // not trying to do better anymore } else if(!warned && (log_ ? fabs(y - la) > 3 : fabs(y - a) > 1e-4)) { if(!(log_ && y == ML_NEGINF && // e.g. qbeta(-1e-10, .2, .03, log=TRUE) cannot get accurate ==> do NOT warn pbeta_raw(DBL_1__eps, // = 1 - eps pp, qq, TRUE, TRUE) > la + 2)) MATHLIB_WARNING2( // low accuracy for more platform independent output: "qbeta(a, *) =: x0 with |pbeta(x0,*%s) - alpha| = %.5g is not accurate", (log_ ? ", log_" : ""), fabs(y - (log_ ? la : a))); } L_return: // use if (use_log_x) u_n else tx {and u_n is on log scale} if(give_log_q) { // {currently not used from R's qbeta()} ==> use_log_x , too if(!use_log_x) // (see if claim above is true) MATHLIB_WARNING( "qbeta() L_return, u_n=%g; give_log_q=TRUE but use_log_x=FALSE -- please report!", u_n); double r = R_Log1_Exp(u_n); if(swap_tail) { qb[0] = r; qb[1] = u_n; } else { qb[0] = u_n; qb[1] = r; } } else { if(use_log_x) { if(add_N_step) { /* add one last Newton step on original x scale, e.g., for qbeta(2^-98, 0.125, 2^-96) */ if(u_n != 1.) // u_n has been computed above xinbta = exp(u_n); y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, log_p); w = log_p ? (y - la) * exp(y + logbeta + r * log(xinbta) + t * log1p(-xinbta)) : (y - a) * exp( logbeta + r * log(xinbta) + t * log1p(-xinbta)); if(R_FINITE(w)) { tx = xinbta - w; R_ifDEBUG_printf(" Final Newton correction(non-log scale):\n" " xinbta=%.16g, y=%g, w=-Delta(x)=%g. \n=> new x=%.16g\n", xinbta, y, w, tx); } else { // Newton step w cannot be used R_ifDEBUG_printf(" Final Newton correction(non-log scale), canNOT use step 'w':\n" " xinbta=%.16g, y=%g, w=-Delta(x)=%g\n=> keep x=xinbta\n", xinbta, y, w); tx = xinbta; } } else { if(swap_tail) { qb[0] = -expm1(u_n); qb[1] = exp (u_n); } else { qb[0] = exp (u_n); qb[1] = -expm1(u_n); } return; } } if(swap_tail) { qb[0] = 1 - tx; qb[1] = tx; } else { qb[0] = tx; qb[1] = 1 - tx; } } return; }