qnorm_R.h

Go to the documentation of this file.

/*
 *  Mathlib : A C Library of Special Functions
 *  Copyright (C) 2000--2023 The R Core Team
 *  Copyright (C) 1998       Ross Ihaka
 *  based on AS 241 (C) 1988 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/
 *
 *  SYNOPSIS
 *
 *  double qnorm5(double p, double mu, double sigma,
 *            int lower_tail, int log_p)
 *            {qnorm (..) is synonymous and preferred inside R}
 *
 *  DESCRIPTION
 *
 *  Compute the quantile function for the normal distribution.
 *
 *  The algorithm AS 241 of Wichura is used,
 *      and has been improved for the very extreme tail (and log_p=TRUE)
 *
 *  REFERENCE
 *
 *      Wichura, M.J. (1988).
 *      Algorithm AS 241: The Percentage Points of the Normal Distribution.
 *      Applied Statistics, 37, 477-484.
 *
 *      Maechler, M. (2022). Asymptotic tail formulas for gaussian quantiles;
 *      https://CRAN.R-project.org/package=DPQ/vignettes/qnorm-asymp.pdf
 */
 
 
// Define here to avoid depending on R
 
#include <cmath>
 
namespace glm {
 
#ifndef ML_NEGINF
#define ML_NEGINF std::numeric_limits<double>::min()
#endif
 
 
#ifndef ML_POSINF   
#define ML_POSINF std::numeric_limits<double>::max()
#endif
 
#ifndef NAN   
#define NAN std::numeric_limits<double>::quiet_NaN()
#endif
 
 
#ifndef R_D_Lval 
#define R_D_Lval(p) (lower_tail ? (p) : (0.5 - (p) + 0.5))  /*  p  */
#endif
 
 
 
#ifndef R_DT_qIv 
#define R_DT_qIv(p) (log_p ? (lower_tail ? exp(p) : - expm1(p)) \
                   : R_D_Lval(p))
#endif
 
 
/* Do the boundaries exactly for q*() functions :
 * Often  _LEFT_ = ML_NEGINF , and very often _RIGHT_ = ML_POSINF;
 *
 * R_Q_P01_boundaries(p, _LEFT_, _RIGHT_)  :<==>
 *
 *     R_Q_P01_check(p);
 *     if (p == R_DT_0) return _LEFT_ ;
 *     if (p == R_DT_1) return _RIGHT_;
 *
 * the following implementation should be more efficient (less tests):
 */
#ifndef R_Q_P01_boundaries 
#define R_Q_P01_boundaries(p, _LEFT_, _RIGHT_)      \
    if (log_p) {                    \
    if(p > 0)                   \
        NAN;              \
    if(p == 0) /* upper bound*/         \
        return lower_tail ? _RIGHT_ : _LEFT_;   \
    if(p == ML_NEGINF)              \
        return lower_tail ? _LEFT_ : _RIGHT_;   \
    }                           \
    else { /* !log_p */                 \
    if(p < 0 || p > 1)              \
        NAN;              \
    if(p == 0)                  \
        return lower_tail ? _LEFT_ : _RIGHT_;   \
    if(p == 1)                  \
        return lower_tail ? _RIGHT_ : _LEFT_;   \
    }
#endif    
 
static double qnorm(double p, double mu, double sigma, int lower_tail, int log_p){
    double p_, q, r, val;
 
// #ifdef IEEE_754
    if (std::isnan(p) || std::isnan(mu) || std::isnan(sigma))
    return p + mu + sigma;
// #endif
    R_Q_P01_boundaries(p, ML_NEGINF, ML_POSINF);
 
    if(sigma  < 0)  NAN;
    if(sigma == 0)  return mu;
 
    p_ = R_DT_qIv(p);/* real lower_tail prob. p */
    q = p_ - 0.5;
 
// #ifdef DEBUG_qnorm
//     REprintf("qnorm(p=%10.7g, m=%g, s=%g, l.t.= %d, log= %d): q = %g\n",
//       p,mu,sigma, lower_tail, log_p, q);
// #endif
 
 
/*-- use AS 241 --- */
/* double ppnd16_(double *p, long *ifault)*/
/*      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
 
        Produces the normal deviate Z corresponding to a given lower
        tail area of P; Z is accurate to about 1 part in 10**16.
 
        (original fortran code used PARAMETER(..) for the coefficients
         and provided hash codes for checking them...)
*/
    if (fabs(q) <= .425) {/* |p~ - 0.5| <= .425  <==> 0.075 <= p~ <= 0.925 */
        r = .180625 - q * q; // = .425^2 - q^2  >= 0
    val =
            q * (((((((r * 2509.0809287301226727 +
                       33430.575583588128105) * r + 67265.770927008700853) * r +
                     45921.953931549871457) * r + 13731.693765509461125) * r +
                   1971.5909503065514427) * r + 133.14166789178437745) * r +
                 3.387132872796366608)
            / (((((((r * 5226.495278852854561 +
                     28729.085735721942674) * r + 39307.89580009271061) * r +
                   21213.794301586595867) * r + 5394.1960214247511077) * r +
                 687.1870074920579083) * r + 42.313330701600911252) * r + 1.);
    }
    else { /* closer than 0.075 from {0,1} boundary :
        *  r := log(p~);  p~ = min(p, 1-p) < 0.075 :  */
    double lp;
    if(log_p && ((lower_tail && q <= 0) || (!lower_tail && q > 0))) {
        lp = p;
    } else {
        lp = log( (q > 0) ? R_DT_CIv(p) /* 1-p */ : p_ /* = R_DT_Iv(p) ^=  p */);
    }
    // r = sqrt( - log(min(p,1-p)) )  <==>  min(p, 1-p) = exp( - r^2 ) :
        r = sqrt(-lp);
#ifdef DEBUG_qnorm
    REprintf("\t close to 0 or 1: r = %7g\n", r);
#endif
        if (r <= 5.) { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */
            r += -1.6;
            val = (((((((r * 7.7454501427834140764e-4 +
                       .0227238449892691845833) * r + .24178072517745061177) *
                     r + 1.27045825245236838258) * r +
                    3.64784832476320460504) * r + 5.7694972214606914055) *
                  r + 4.6303378461565452959) * r +
                 1.42343711074968357734)
                / (((((((r *
                         1.05075007164441684324e-9 + 5.475938084995344946e-4) *
                        r + .0151986665636164571966) * r +
                       .14810397642748007459) * r + .68976733498510000455) *
                     r + 1.6763848301838038494) * r +
                    2.05319162663775882187) * r + 1.);
        }
    else if(r <= 27) { /* p is very close to  0 or 1: r in (5, 27] :
        *  r >   5 <==> min(p,1-p)  < exp(-25) = 1.3888..e-11
        *  r <= 27 <==> min(p,1-p) >= exp(-27^2) = exp(-729) ~= 2.507972e-317
        * i.e., we are just barely in the range where min(p, 1-p) has not yet underflowed to zero.
        */
        // Wichura, p.478: minimax rational approx R_3(t) is for 5 <= t <= 27  (t :== r)
            r += -5.;
            val = (((((((r * 2.01033439929228813265e-7 +
                       2.71155556874348757815e-5) * r +
                      .0012426609473880784386) * r + .026532189526576123093) *
                    r + .29656057182850489123) * r +
                   1.7848265399172913358) * r + 5.4637849111641143699) *
                 r + 6.6579046435011037772)
                / (((((((r *
                         2.04426310338993978564e-15 + 1.4215117583164458887e-7)*
                        r + 1.8463183175100546818e-5) * r +
                       7.868691311456132591e-4) * r + .0148753612908506148525)
                     * r + .13692988092273580531) * r +
                    .59983220655588793769) * r + 1.);
        }
        else { // r > 27: p is *really* close to 0 or 1 .. practically only when log_p =TRUE
        if(r >= 6.4e8) { // p is *very extremely* close to 0 or 1
        // Using the asymptotical formula ("0-th order"): qn = sqrt(2*s)
        val = r * M_SQRT2;
        } else {
        double s2 = -ldexp(lp, 1), // = -2*lp = 2s
            x2 = s2 - log(M_2PI * s2); // = xs_1
        // if(r >= 36000.)  # <==> s >= 36000^2   use x2 = xs_1  above
        if(r < 36000.) {
            x2 = s2 - log(M_2PI * x2) - 2./(2. + x2); // == xs_2
            if(r < 840.) { // 27 < r < 840
            x2 = s2 - log(M_2PI * x2) + 2*log1p(- (1 - 1/(4 + x2))/(2. + x2)); // == xs_3
            if(r < 109.) { // 27 < r < 109
              x2 = s2 - log(M_2PI * x2) +
                  2*log1p(- (1 - (1 - 5/(6 + x2))/(4. + x2))/(2. + x2)); // == xs_4
              if(r < 55.) { // 27 < r < 55
                x2 = s2 - log(M_2PI * x2) +
                  2*log1p(- (1 - (1 - (5 - 9/(8. + x2))/(6. + x2))/(4. + x2))/(2. + x2)); // == xs_5
              }
            }
            }
        }
                val = sqrt(x2);
        }
    }
    if(q < 0.0)
        val = -val;
    }
    return mu + sigma * val;
}
 
};