Package 'DPQ'

Description

Computations for approximations and alternatives for the 'DPQ' (Density (pdf), Probability (cdf) and Quantile) functions for probability distributions in R. Primary focus is on (central and non-central) beta, gamma and related distributions such as the chi-squared, F, and t. – For several distribution functions, provide functions implementing formulas from Johnson, Kotz, and Kemp (1992) <doi:10.1002/bimj.4710360207> and Johnson, Kotz, and Balakrishnan (1995) for discrete or continuous distributions respectively. This is for the use of researchers in these numerical approximation implementations, notably for my own use in order to improve standard R pbeta(), qgamma(), ..., etc: {'"dpq"'-functions}.

Details

The DESCRIPTION file:

Index of help topics:

.D_0                    Distribution Utilities "dpq"
Bern                    Bernoulli Numbers
DPQ-package             Density, Probability, Quantile ('DPQ')
                        Computations
Ixpq                    Normalized Incomplete Beta Function "Like"
                        'pbeta()'
M_LN2                   Numerical Utilities - Functions, Constants
algdiv                  Compute log(gamma(b)/gamma(a+b)) when b >= 8
b_chi                   Compute E[chi_nu]/sqrt(nu) useful for t- and
                        chi-Distributions
bd0                     Binomial Deviance - Auxiliary Functions for
                        'dgamma()' Etc
bpser                   'pbeta()' 'bpser' series computation
chebyshevPoly           Chebyshev Polynomial Evaluation
dbinom_raw              R's C Mathlib (Rmath) dbinom_raw() Binomial
                        Probability pure R Function
dgamma.R                Gamma Density Function Alternatives
dhyperBinMolenaar       HyperGeometric (Point) Probabilities via
                        Molenaar's Binomial Approximation
dnbinomR                Pure R Versions of R's C (Mathlib) dnbinom()
                        Negative Binomial Probabilities
dnchisqR                Approximations of the (Noncentral) Chi-Squared
                        Density
dntJKBf1                Non-central t-Distribution Density - Algorithms
                        and Approximations
dpsifn                  Psi Gamma Functions Workhorse from R's API
dtWV                    Asymptotic Noncentral t Distribution Density by
                        Viechtbauer
expm1x                  Accurate exp(x) - 1 - x (for smallish |x|)
format01prec            Format Numbers in [0,1] with "Precise" Result
frexp                   Base-2 Representation and Multiplication of
                        Numbers
gam1d                   Compute 1/Gamma(x+1) - 1 Accurately
gamln1                  Compute log( Gamma(x+1) ) Accurately in [-0.2,
                        1.25]
gammaVer                Gamma Function Versions
hyper2binomP            Transform Hypergeometric Distribution
                        Parameters to Binomial Probability
lbetaM                  (Log) Beta and Ratio of Gammas Approximations
lfastchoose             R versions of Simple Formulas for Logarithmic
                        Binomial Coefficients
lgamma1p                Accurate 'log(gamma(a+1))'
lgammaAsymp             Asymptotic Log Gamma Function
log1mexp                Compute log(1 - exp(-a)) and log(1 + exp(x))
                        Numerically Optimally
log1pmx                 Accurate 'log(1+x) - x' Computation
logcf                   Continued Fraction Approximation of Log-Related
                        Power Series
logspace.add            Logspace Arithmetix - Addition and Subtraction
lssum                   Compute Logarithm of a Sum with Signed Large
                        Summands
lsum                    Properly Compute the Logarithm of a Sum (of
                        Exponentials)
newton                  Simple R level Newton Algorithm, Mostly for
                        Didactical Reasons
p1l1                    Numerically Stable p1l1(t) = (t+1)*log(1+t) - t
pbetaRv1                Pure R Implementation of Old pbeta()
pchisqV                 Wienergerm Approximations to (Non-Central)
                        Chi-squared Probabilities
phyper1molenaar         Molenaar's Normal Approximations to the
                        Hypergeometric Distribution
phyperAllBin            Compute Hypergeometric Probabilities via
                        Binomial Approximations
phyperApprAS152         Normal Approximation to cumulative Hyperbolic
                        Distribution - AS 152
phyperBin.1             HyperGeometric Distribution via Approximate
                        Binomial Distribution
phyperBinMolenaar       HyperGeometric Distribution via Molenaar's
                        Binomial Approximation
phyperIbeta             Pearson's incomplete Beta Approximation to the
                        Hyperbolic Distribution
phyperPeizer            Peizer's Normal Approximation to the Cumulative
                        Hyperbolic
phyperR                 R-only version of R's original phyper()
                        algorithm
phyperR2                Pure R version of R's C level phyper()
phypers                 The Four (4) Symmetric 'phyper()' Calls
pl2curves               Plot 2 Noncentral Distribution Curves for
                        Visual Comparison
pnbetaAppr2             Noncentral Beta Probabilities
pnchi1sq                (Probabilities of Non-Central Chi-squared
                        Distribution for Special Cases
pnchisq                 (Approximate) Probabilities of Non-Central
                        Chi-squared Distribution
pnormAsymp              Asymptotic Approxmation of (Extreme Tail)
                        'pnorm()'
pnormL_LD10             Bounds for 1-Phi(.) - Mill's Ratio related
                        Bounds for pnorm()
pntR                    Non-central t Probability Distribution -
                        Algorithms and Approximations
pow                     X to Power of Y - R C API 'R_pow()'
pow1p                   Accurate (1+x)^y, notably for small |x|
ppoisErr                Direct Computation of 'ppois()' Poisson
                        Distribution Probabilities
pt_Witkovsky_Tab1       Viktor Witosky's Table_1 pt() Examples
qbetaAppr               Compute (Approximate) Quantiles of the Beta
                        Distribution
qbinomR                 Pure R Implementation of R's qbinom() with
                        Tuning Parameters
qchisqAppr              Compute Approximate Quantiles of the
                        Chi-Squared Distribution
qgammaAppr              Compute (Approximate) Quantiles of the Gamma
                        Distribution
qnbinomR                Pure R Implementation of R's qnbinom() with
                        Tuning Parameters
qnchisqAppr             Compute Approximate Quantiles of Noncentral
                        Chi-Squared Distribution
qnormAppr               Approximations to 'qnorm()', i.e., z_alpha
qnormAsymp              Asymptotic Approximation to Outer Tail of
                        qnorm()
qnormR                  Pure R version of R's 'qnorm()' with
                        Diagnostics and Tuning Parameters
qntR                    Pure R Implementation of R's qt() / qnt()
qpoisR                  Pure R Implementation of R's qpois() with
                        Tuning Parameters
qtAppr                  Compute Approximate Quantiles of the
                        (Non-Central) t-Distribution
qtR                     Pure R Implementation of R's C-level
                        t-Distribution Quantiles 'qt()'
qtU                     'uniroot()'-based Computing of t-Distribution
                        Quantiles
r_pois                  Compute Relative Size of i-th term of Poisson
                        Distribution Series
rexpm1                  TOMS 708 Approximation REXP(x) of expm1(x) =
                        exp(x) - 1
stirlerr                Stirling's Error Function - Auxiliary for
                        Gamma, Beta, etc

An important goal is to investigate diverse algorithms and approximations of R's own density (d*()), probability (p*()), and quantile (q*()) functions, notably in “border” cases where the traditional published algorithms have shown to be suboptimal, not quite accurate, or even useless.

Examples are border cases of the beta distribution, or non-central distributions such as the non-central chi-squared and t-distributions.

Author(s)

Principal author and maintainer: Martin Maechler <[email protected]>

See Also

The package DPQmpfr, which builds on this package and on Rmpfr.

Examples

## Show problem in R's non-central t-distrib. density (and approximations):
example(dntJKBf)

Description

Computes

$\code{algdiv(a,b)} := \log \frac{\Gamma(b)}{\Gamma(a+b)} = \log \Gamma(b) - \log\Gamma(a+b) = \code{lgamma(b) - lgamma(a+b)}$

in a numerically stable way.

This is an auxiliary function in R's (TOMS 708) implementation of pbeta(), aka the incomplete beta function ratio.

Usage

algdiv(a, b)

Arguments

Details

Note that this is also useful to compute the Beta function

$B(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}.$

Clearly,

$\log B(a,b) = \log\Gamma(a) + \mathrm{algdiv(a,b)} = \log\Gamma(a) - \mathrm{logQab}(a,b)$

In our ‘../tests/qbeta-dist.R’ we look into computing $\log(p B(p,q))$ accurately for $p \ll q$ .

We are proposing a nice solution there.
How is this related to algdiv() ?

Additionally, we have defined

$Qab = Q_{a,b} := \frac{\Gamma(a+b),\Gamma(b)},$

such that $\code{logQab(a,b)} := \log Qab(a,b)$ fulfills simply

$\code{logQab(a,b)} = - \code{algdiv(a,b)}$

see logQab_asy.

Value

a numeric vector of length max(length(a), length(b)) (if neither is of length 0, in which case the result has length 0 as well).

Author(s)

Didonato, A. and Morris, A., Jr, (1992); algdiv()'s C version from the R sources, authored by the R core team; C and R interface: Martin Maechler

References

Didonato, A. and Morris, A., Jr, (1992) Algorithm 708: Significant digit computation of the incomplete beta function ratios, ACM Transactions on Mathematical Software 18, 360–373.

See Also

gamma, beta; my own logQab_asy().

Examples

Qab <- algdiv(2:3, 8:14)
cbind(a = 2:3, b = 8:14, Qab) # recycling with a warning

## algdiv()  and my  logQab_asy()  give *very* similar results for largish b:
all.equal( -   algdiv(3, 100),
           logQab_asy(3, 100), tolerance=0) # 1.283e-16 !!
(lQab <- logQab_asy(3, 1e10))
## relative error
1 + lQab/ algdiv(3, 1e10) # 0 (64b F 30 Linux; 2019-08-15)

## in-and outside of "certified" argument range {b >= 8}:
a. <- c(1:3, 4*(1:8))/32
b. <- seq(1/4, 20, by=1/4)
ad <- t(outer(a., b., algdiv))
## direct computation:
f.algdiv <- function(a,b) lgamma(b) - lgamma(a+b)
ad.d <- t(outer(a., b., f.algdiv))

matplot (b., ad.d, type = "o", cex=3/4,
         main = quote(log(Gamma(b)/Gamma(a+b)) ~"  vs.  algdiv(a,b)"))
mtext(paste0("a[1:",length(a.),"] = ",
        paste0(paste(head(paste0(formatC(a.*32), "/32")), collapse=", "), ", .., 1")))
matlines(b., ad,   type = "l", lwd=4, lty=1, col=adjustcolor(1:6, 1/2))
abline(v=1, lty=3, col="midnightblue")
# The larger 'b', the more accurate the direct formula wrt algdiv()
all.equal(ad[b. >= 1,], ad.d[b. >= 1,]       )# 1.5e-5
all.equal(ad[b. >= 2,], ad.d[b. >= 2,], tol=0)# 3.9e-9
all.equal(ad[b. >= 4,], ad.d[b. >= 4,], tol=0)# 4.6e-13
all.equal(ad[b. >= 6,], ad.d[b. >= 6,], tol=0)# 3.0e-15
all.equal(ad[b. >= 8,], ad.d[b. >= 8,], tol=0)# 2.5e-15 (not much better)

Description

$b_\chi(\nu) := E[\chi(\nu)] / \sqrt{\nu} = \frac{\sqrt{2/\nu}\Gamma((\nu+1)/2)}{\Gamma(\nu/2)},$

where $\chi(\nu)$ denotes a chi-distributed random variable, i.e., the square of a chi-squared variable, and $\Gamma(z)$ is the Gamma function, gamma() in R.

This is a relatively important auxiliary function when computing with non-central t distribution functions and approximations, specifically see Johnson et al.(1994), p.520, after (31.26a), e.g., our pntJW39().

Its logarithm,

$lb_\chi(\nu) := log\bigl(\frac{\sqrt{2/\nu}\Gamma((\nu+1)/2)}{\Gamma(\nu/2)}\bigr),$

is even easier to compute via lgamma and log, and I have used Maple to derive an asymptotic expansion in $\frac{1}{\nu}$ as well.

Note that $lb_\chi(\nu)$ also appears in the formula for the t-density (dt) and distribution (tail) functions.

Usage

b_chi     (nu, one.minus = FALSE, c1 = 341, c2 = 1000)
b_chiAsymp(nu, order = 2, one.minus = FALSE)
#lb_chi    (nu, ......) # not yet
lb_chiAsymp(nu, order)

c_dt(nu)       # warning("FIXME: current c_dt() is poor -- base it on lb_chi(nu) !")
c_dtAsymp(nu)  # deprecated in favour of lb_chi(nu)
c_pt(nu)       # warning("use better c_dt()") %---> FIXME deprecate even stronger ?

Arguments

Details

One can see that b_chi() has the properties of a CDF of a continuous positive random variable: It grows monotonely from $b_\chi(0) = 0$ to (asymptotically) one. Specifically, for large nu, b_chi(nu) = b_chiAsymp(nu) and

$1 - b_\chi(\nu) \sim \frac{1}{4\nu}.$

More accurately, derived from Abramowitz and Stegun, 6.1.47 (p.257) for a= 1/2, b=0,

$\Gamma(z + 1/2) / \Gamma(z) \sim \sqrt(z)*(1 - 1/(8z) + 1/(128 z^2) + O(1/z^3)),$

and applied for $b_\chi(\nu)$ with $z = \nu/2$, we get

$b_\chi(\nu) \sim 1 - (1/(4\nu) * (1 - 1/(8\nu)) + O(\nu^{-3})),$

which has been implemented in b_chiAsymp(*, order=2) in 1999.

Even more accurately, Martin Maechler, used Maple to derive an asymptotic expansion up to order 15, here reported up to order 5, namely with $r := \frac{1}{4\nu}$,

$b_\chi(\nu) = c_\chi(r) = 1 - r + \frac{1}{2}r^2 + \frac{5}{2}r^3 - \frac{21}{8}r^4 - \frac{399}{8}r^5 + O(r^6).$

Value

a numeric vector of the same length as nu.

Author(s)

Martin Maechler

References

Johnson, Kotz, Balakrishnan (1995) Continuous Univariate Distributions, Vol 2, 2nd Edition; Wiley; Formula on page 520, after (31.26a)

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain.

See Also

The t-distribution (base R) page pt; our pntJW39().

Examples

curve(b_chi, 0, 20); abline(h=0:1, v=0, lty=3)
r <- curve(b_chi, 1e-10, 1e5, log="x")
with(r, lines(x, b_chi(x, one.minus=TRUE), col = 2))

## Zoom in to c1-region
rc1 <- curve(b_chi, 340.5, 341.5, n=1001)# nothing to see
e <- 1e-3; curve(b_chi, 341-e, 341+e, n=1001) # nothing
e <- 1e-5; curve(b_chi, 341-e, 341+e, n=1001) # see noise, but no jump
e <- 1e-7; curve(b_chi, 341-e, 341+e, n=1001) # see float "granularity"+"jump"

## Zoom in to c2-region
rc2 <- curve(b_chi, 999.5, 1001.5, n=1001) # nothing visible
e <- 1e-3; curve(b_chi, 1000-e, 1000+e, n=1001) # clear small jump
c2 <- 1500
e <- 1e-3; curve(b_chi(x,c2=c2), c2-e, c2+e, n=1001)# still
## - - - -
c2 <- 3000
e <- 1e-3; curve(b_chi(x,c2=c2), c2-e, c2+e, n=1001)# ok asymp clearly better!!
curve(b_chiAsymp, add=TRUE, col=adjustcolor("red", 1/3), lwd=3)
if(requireNamespace("Rmpfr")) {
 xm <- Rmpfr::seqMpfr(c2-e, c2+e, length.out=1000)

}
## - - - -
c2 <- 4000
e <- 1e-3; curve(b_chi(x,c2=c2), c2-e, c2+e, n=1001)# ok asymp clearly better!!
curve(b_chiAsymp, add=TRUE, col=adjustcolor("red", 1/3), lwd=3)

grCol <- adjustcolor("forest green", 1/2)
curve(b_chi,                    1/2, 1e11, log="x")
curve(b_chiAsymp, add = TRUE, col = grCol, lwd = 3)
## 1-b(nu) ~= 1/(4 nu) a power function <==> linear in log-log scale:
curve(b_chi(x, one.minus=TRUE), 1/2, 1e11, log="xy")
curve(b_chiAsymp(x, one.minus=TRUE), add = TRUE, col = grCol, lwd = 3)

Description

Return the $n$-th Bernoulli number $B_n$, (or $B_n^+$, see the reference), where $B_1 = + \frac 1 2$.

Usage

Bern(n, verbose = getOption("verbose", FALSE))

Arguments

Value

The number $B_n$ of type numeric.

A side effect is the caching of computed Bernoulli numbers in the hidden environment .bernoulliEnv.

Author(s)

Martin Maechler

References

https://en.wikipedia.org/wiki/Bernoulli_number

See Also

Bernoulli in Rmpfr in arbitrary precision via Riemann's $\zeta$ function.

The next version of package gmp is to contain BernoulliQ(), providing exact Bernoulli numbers as big rationals (class "bigq").

Examples

(B.0.10 <- vapply(0:10, Bern, 1/2))
## [1]  1.00000000 +0.50000000  0.16666667  0.00000000 -0.03333333  0.00000000
## [7]  0.02380952  0.00000000 -0.03333333  0.00000000  0.07575758
if(requireNamespace("MASS")) {
  print( MASS::fractions(B.0.10) )
  ## 1  +1/2   1/6    0  -1/30     0  1/42     0 -1/30     0  5/66
}

Description

Compute the bpser series approximation of pbeta, the incomplete beta function. Note that when b is integer valued, the series is a sum of $b+1$ terms.

Usage

bpser(a, b, x, log.p = FALSE, eps = 1e-15, verbose = FALSE, warn = TRUE)

Arguments

Value

a list with components

Author(s)

Martin Maechler, ported to DPQ; R-Core team for the code in R.

References

TOMS 708, see pbeta

See Also

R's pbeta; DPQ's pbetaRv1(), and Ixpq(); Rmpfr's pbetaI

Examples

with(bpser(100000, 11, (0:64)/64), # all 0 {last one "wrongly"}
     stopifnot(r == c(rep(0, 64), 1), err == 0))
bp1e5.11L <- bpser(100000, 11, (0:64)/64, log.p=TRUE)# -> 2 "underflow to -Inf" warnings!
pbe <- pbeta((0:64)/64, 100000, 11, log.p=TRUE)

## verbose=TRUE showing info on number of terms / iterations
ps11.5 <- bpser(100000, 11.5, (0:64)/64, log.p=TRUE, verbose=TRUE)

Description

Provides (evaluation of) Chebyshev polynomials, given their coefficients vector coef (using $2 c_0$, i.e., 2*coef[1] as the base R mathlib chebyshev*() functions. Specifically, the following sum is evaluated:

$\sum_{j=0}^n c_j T_j(x)$

where $c_0 :=$coef[1] and $c_j :=$coef[j+1] for $j \ge 1$. $n :=$ chebyshev_nc(coef, .) is the maximal degree and hence one less than the number of terms, and $T_j()$ is the Chebyshev polynomial (of the first kind) of degree $j$.

Usage

chebyshevPoly(coef, nc = chebyshev_nc(coef, eta), eta = .Machine$double.eps/20)

chebyshev_nc(coef, eta = .Machine$double.eps/20)
chebyshevEval(x, coef,
              nc = chebyshev_nc(coef, eta), eta = .Machine$double.eps/20)

Arguments

Value

chebyshevPoly() returns function(x) which computes the values of the underlying Chebyshev polynomial at x.

chebyshev_nc() returns an integer, and chebyshevEval(x, coef) returns a numeric “like” x with the values of the polynomial at x.

Author(s)

R Core team, notably Ross Ihaka; Martin Maechler provided the R interface.

References

https://en.wikipedia.org/wiki/Chebyshev_polynomials

See Also

polyn.eval() from CRAN package sfsmisc; as one example of many more.

Examples

## The first 5 (base) Chebyshev polynomials:
T0 <- chebyshevPoly(2)  # !! 2, not 1
T1 <- chebyshevPoly(0:1)
T2 <- chebyshevPoly(c(0,0,1))
T3 <- chebyshevPoly(c(0,0,0,1))
T4 <- chebyshevPoly(c(0,0,0,0,1))
curve(T0(x), -1,1, col=1, lwd=2, ylim=c(-1,1))
abline(h=0, lty=2)
curve(T1(x), col=2, lwd=2, add=TRUE)
curve(T2(x), col=3, lwd=2, add=TRUE)
curve(T3(x), col=4, lwd=2, add=TRUE)
curve(T4(x), col=5, lwd=2, add=TRUE)

(Tv <- vapply(c(T0=T0, T1=T1, T2=T2, T3=T3, T4=T4),
              function(Tp) Tp(-1:1), numeric(3)))
x <- seq(-1,1, by = 1/64)
stopifnot(exprs = {
   all.equal(chebyshevPoly(1:5)(x),
             0.5*T0(x) + 2*T1(x) + 3*T2(x) + 4*T3(x) + 5*T4(x))
   all.equal(unname(Tv), rbind(c(1,-1), c(1:-1,0:1), rep(1,5)))# warning on rbind()
})

Description

A pure R implementation of R's C API (‘Mathlib’ specifically) dbinom_raw() function which computes binomial probabilities and is continuous in x, i.e., also “works” for non-integer x.

Usage

dbinom_raw (x, n, p, q = 1-p, log = FALSE,
            version = c("2008", "R4.4"),
            verbose = getOption("verbose"))

Arguments

Value

numeric vector of the same length as x which may have to be thought of recycled along n, p and/or q.

Author(s)

R Core and Martin Maechler

See Also

Note that our CRAN package Rmpfr provides dbinom, an mpfr-accurate function to be used used instead of R's or this pure R version relying bd0() and stirlerr() where the latter currently only provides accurate double precision accuracy.

Examples

for(n in c(3, 10, 27, 100, 500, 2000, 5000, 1e4, 1e7, 1e10)) {
 x <- if(n <= 2000) 0:n else round(seq(0, n, length.out=2000))
 p <- 3/4
 stopifnot(all.equal(dbinom_raw(x, n, p, q=1-p) -> dbin,
                     dbinom    (x, n, p), tolerance = 1e-13))# 1.636e-14 (Apple clang 14.0.3)
 stopifnot(all.equal(dbin, dbinom_raw(x, n, p, q=1-p, version = "R4.4") -> dbin44,
                     tolerance = 1e-13))
 cat("n = ", n, ": ", (aeq <- all.equal(dbin44, dbin, tolerance = 0)), "\n")
 if(n < 3000) stopifnot(is.character(aeq)) # check that dbin44 is "better" ?!
}

n <- 1024 ; x <- 0:n
plot(x, dbinom_raw(x, n, p, q=1-p) - dbinom(x, n, p), type="l", main = "|db_r(x) - db(x)|")
plot(x, dbinom_raw(x, n, p, q=1-p) / dbinom(x, n, p) - 1, type="b", log="y",
     main = "rel.err.  |db_r(x / db(x) - 1)|")

Description

Compute the density function $f(x, *)$ of the (noncentral) chi-squared distribution.

Usage

dnchisqR     (x, df, ncp, log = FALSE,
              eps = 5e-15, termSml = 1e-10, ncpLarge = 1000)
dnchisqBessel(x, df, ncp, log = FALSE)
dchisqAsym   (x, df, ncp, log = FALSE)
dnoncentchisq(x, df, ncp, kmax = floor(ncp/2 + 5 * (ncp/2)^0.5))

Arguments

Details

dnchisqR() is a pure R implementation of R's own C implementation in the sources, ‘R/src/nmath/dnchisq.c’, additionally exposing the three “tuning parameters” eps, termSml, and ncpLarge.

dnchisqBessel() implements Fisher(1928)'s exact closed form formula based on the Bessel function $I_{nu}$, i.e., R's besselI() function; specifically formula (29.4) in Johnson et al. (1995).

dchisqAsym() is the simple asymptotic approximation from Abramowitz and Stegun's formula 26.4.27, p. 942.

dnoncentchisq() uses the (typically defining) infinite series expansion directly, with truncation at kmax, and terms $t_k$ which are products of a Poisson probability and a central chi-square density, i.e., terms t.k := dpois(k, lambda = ncp/2) * dchisq(x, df = 2*k + df) for k = 0, 1, ..., kmax.

Value

numeric vector similar to x, containing the (logged if log=TRUE) values of the density $f(x,*)$.

Note

These functions are mostly of historical interest, notably as R's dchisq() was not always very accurate in the noncentral case, i.e., for ncp > 0.

Note

R's dchisq() is typically more uniformly accurate than the approximations nowadays, apart from dnchisqR() which should behave the same. There may occasionally exist small differences between dnchisqR(x, *) and dchisq(x, *) for the same parameters.

Author(s)

Martin Maechler, April 2008

References

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain.

Johnson, N.L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions Vol~2, 2nd ed.; Wiley; chapter 29, Section 3 Distribution, (29.4), p. 436.

See Also

R's own dchisq().

Examples

x <- sort(outer(c(1,2,5), 2^(-4:5)))
fRR <- dchisq  (x, 10, 2)
f.R <- dnchisqR(x, 10, 2)
all.equal(fRR, f.R, tol = 0) # 64bit Lnx (F 30): 1.723897e-16
stopifnot(all.equal(fRR, f.R, tol = 4e-15))

Description

The “binomial deviance” function bd0(x, M) := $D_0(x,M) := M \cdot d_0(x/M)$, where $d_0(r) := r\log(r) + 1-r$.

Mostly, pure R transcriptions of the C code utility functions for dgamma(), dbinom(), dpois(), dt(), and similar “base” density functions by Catherine Loader.
These have extra arguments with defaults that correspond to R's Mathlib C code hardwired cutoffs and tolerances.

Usage

dpois_raw(x, lambda, log=FALSE,
          version, 
          small.x__lambda = .Machine$double.eps,
          ## the defaults for version will probably change in the future
          bd0.delta = 0.1,
          ## optional arguments of log1pmx() :
          tol_logcf = 1e-14, eps2 = 0.01, minL1 = -0.79149064, trace.lcf = verbose,
          logCF = if (is.numeric(x)) logcf else logcfR,
          verbose = FALSE)

dpois_simpl (x, lambda, log=FALSE)
dpois_simpl0(x, lambda, log=FALSE)

bd0(x, np,
    delta = 0.1, maxit = as.integer(-1100 / log2(delta)),
    s0 = .Machine$double.xmin,
    verbose = getOption("verbose"))
bd0C(x, np, delta = 0.1, maxit = 1000L, version = "R4.0", verbose = getOption("verbose"))
# "simple" log1pmx() based versions :
bd0_p1l1d1(x, M, tol_logcf = 1e-14, ...)
bd0_p1l1d (x, M, tol_logcf = 1e-14, ...)
bd0_l1pm  (x, M, tol_logcf = 1e-14, ...)

ebd0 (x, M, verbose = getOption("verbose"), ...) # experimental, may disappear !!
ebd0C(x, M, verbose = getOption("verbose"))

Arguments

Details

bd0():

Loader's “Binomial Deviance” function; for $x, M > 0$ (where the limit $x \to 0$ is allowed). In the case of dbinom, $x$ are integers (and $M = n p$), but in general x is real.

$bd_0(x,M) := M \cdot D_0\bigl(\frac{x}{M}\bigr),$

where $D_0(u) := u \log(u) + 1-u = u(\log(u) - 1) + 1$. Hence

$bd_0(x,M) = M \cdot \bigl(\frac{x}{M}(\log(\frac{x}{M}) -1) +1 \bigr) = x \log(\frac{x}{M}) - x + M.$

A different way to rewrite this from Martyn Plummer, notably for important situation when $\left|x-M \right| \ll M$, is using $t := (x-M)/M$ (and $\left|t \right| \ll 1$ for that situation), equivalently, $\frac{x}{M} = 1+t$. Using $t$,

$bd_0(x,M) = \log(1+t) - t \cdot M = M \cdot [(t+1)(\log(1+t) - 1) + 1] = M \cdot [(t+1) \log(1+t) - t] = M \cdot p_1l_1(t),$

and

$p_1l_1(t) := (t+1)\log(1+t) - t = \frac{t^2}{2} - \frac{t^3}{6} ...$

where the Taylor series expansion is useful for small $|t|$.

Note that bd0(x, M) now also works when x and/or M are arbitrary-accurate mpfr-numbers (package Rmpfr).

bd0C() interfaces to C code which corresponds to R's C Mathlib (Rmath) bd0().

Value

a numeric vector “like” x; in some cases may also be an (high accuracy) "mpfr"-number vector, using CRAN package Rmpfr.

ebd0() (R code) and ebd0C() (interface to C code) are experimental, meant to be precision-extended version of bd0(), returning (yh, yl) (high- and low-part of y, the numeric result). In order to work for long vectors x, yh, yl need to be list components; hence we return a two-column data.frame with column names "yh" and "yl".

Author(s)

Martin Maechler

References

C. Loader (2000), see dbinom's documentation.

Our package vignette log1pmx, bd0, stirlerr - Probability Computations in R.

See Also

stirlerr for Stirling's error function, complementing bd0() for computation of Gamma, Beta, Binomial and Poisson probabilities. dgamma, dpois.

Examples

x <- 800:1200
bd0x1k <- bd0(x, np = 1000)
plot(x, bd0x1k, type="l", ylab = "bd0(x, np=1000)")
bd0x1kC <- bd0C(x, np = 1000)
lines(x, bd0x1kC, col=2)
bd0.1d1 <- bd0_p1l1d1(x, 1000)
bd0.1d  <- bd0_p1l1d (x, 1000)
bd0.1pm <- bd0_l1pm  (x, 1000)
stopifnot(exprs = {
    all.equal(bd0x1kC, bd0x1k,  tol=1e-14) # even tol=0 currently ..
    all.equal(bd0x1kC, bd0.1d1, tol=1e-14)
    all.equal(bd0x1kC, bd0.1d , tol=1e-14)
    all.equal(bd0x1kC, bd0.1pm, tol=1e-14)
})

str(log1pmx) ##--> play with  { tol_logcf, eps2, minL1, trace.lcf, logCF }

ebd0x1k <- ebd0 (x, 1000)
exC     <- ebd0C(x, 1000)
stopifnot(all.equal(exC, ebd0x1k, tol=4e-16))
lines(x, rowSums(ebd0x1k), col=adjustcolor(4, 1/2), lwd=4)

x <- 0:250
dp   <- dpois    (x, 48, log=TRUE)# R's 'stats' pkg function
dp.r <- dpois_raw(x, 48, log=TRUE)
all.equal(dp, dp.r, tol = 0) # on Linux 64b, see TRUE
stopifnot(all.equal(dp, dp.r, tol = 1e-14))
## dpois_raw()  versions:
(vers <- eval(formals(dpois_raw)$version))
mv <- sapply(vers, function(v) dpois_raw(x, 48, version=v))
matplot(x, mv, type="h", log="y", main="dpois_raw(x, 48, version=*)") # "fine"

if(all(mv[,"ebd0_C1"] == mv[,"ebd0_v1"])) {
    cat("versions 'ebd0_C1' and 'ebd0_v1' are identical for lambda=48\n")
    mv <- mv[, vers != "ebd0_C1"]
}
## now look at *relative* errors -- need "Rmpfr" for "truth"
if(requireNamespace("Rmpfr")) {

    dM <- Rmpfr::dpois(Rmpfr::mpfr(x, 256), 48)
    asN <- Rmpfr::asNumeric
    relE <- asN(mv / dM - 1)
    cols <- adjustcolor(1:ncol(mv), 1/2)

    mtit <- "relative Errors of dpois_raw(x, 48, version = * )"
    matplot(x, relE, type="l", col=cols, lwd=3, lty=1, main=mtit)
    legend("topleft", colnames(mv), col=cols, lwd=3, bty="n")

    matplot(x, abs(relE), ylim=pmax(1e-18, range(abs(relE))), type="l", log="y",
            main=mtit, col=cols, lwd=2, lty=1, yaxt="n")
    sfsmisc::eaxis(2)
    legend("bottomright", colnames(mv), col=cols, lwd=2, bty="n", ncol=3)
    ee <- c(.5, 1, 2)* 2^-52; eC <- quote(epsilon[C])
    abline(h=ee, lty=2, col="gray", lwd=c(1,2,1))
    axis(4, at=ee[2:3], expression(epsilon[C], 2 * epsilon[C]), col="gray", las=1)
    par(new=TRUE)
    plot(x, asN(dM), type="h", col=adjustcolor("darkgreen", 1/3), axes=FALSE, ann=FALSE)
    stopifnot(abs(relE) < 8e-13) # seen 2.57e-13
}# Rmpfr

Description

dgamma.R() is aimed to be an R level “clone” of R's C level implementation dgamma (from package stats).

Usage

dgamma.R(x, shape, scale = 1, log,
         dpois_r_args = list())

Arguments

Value

numeric vector of the same length as x (which may have to be thought of recycled along shape and/or scale.

Author(s)

Martin Maechler

See Also

(As R's C code) this depends crucially on the “workhorse” function dpois_raw().

Examples

xy  <- curve(dgamma  (x, 12), 0,30) # R's dgamma()
xyR <- curve(dgamma.R(x, 12, dpois_r_args = list(verbose=TRUE)), add=TRUE,
             col = adjustcolor(2, 1/3), lwd=3)
stopifnot(all.equal(xy, xyR, tolerance = 4e-15)) # seen 7.12e-16
## TODO: check *vectorization* in x --> add tests/*.R				___ TODO ___


## From R's  <R>/tests/d-p-q-r-tst-2.R -- replacing dgamma() w/ dgamma.R()
## PR#17577 - dgamma(x, shape)  for shape < 1 (=> +Inf at x=0) and very small x
stopifnot(exprs = {
    all.equal(dgamma.R(2^-1027, shape = .99 , log=TRUE), 7.1127667376, tol=1e-10)
    all.equal(dgamma.R(2^-1031, shape = 1e-2, log=TRUE), 702.8889158,  tol=1e-10)
    all.equal(dgamma.R(2^-1048, shape = 1e-7, log=TRUE), 710.30007699, tol=1e-10)
    all.equal(dgamma.R(2^-1048, shape = 1e-7, scale = 1e-315, log=TRUE),
              709.96858768, tol=1e-10)
})
## R's dgamma() gave all Inf in R <= 3.6.1 [and still there in 32-bit Windows !]

Description

Compute hypergeometric (point) probabilities via Molenaar's binomial approximation, hyper2binomP().

Usage

dhyperBinMolenaar(x, m, n, k, log = FALSE)

Arguments

Value

a numeric vector, with the length the maximum of the lengths of x, m, n, k.

Author(s)

Martin Maechler

References

See those in phyperBinMolenaar.

See Also

hyper2binomP(); R's own dhyper() which uses more sophisticated computations.

Examples

## The function is simply defined as
function (x, m, n, k, log = FALSE)
  dbinom(x, size = k, prob = hyper2binomP(x, m, n, k), log = log)

Description

Compute pure R implementations of R's C Mathlib (Rmath) dnbinom() binomial probabilities, allowing to see the effect of the cutoff eps.

Usage

dnbinomR  (x, size, prob, log = FALSE, eps = 1e-10)
dnbinom.mu(x, size, mu,   log = FALSE, eps = 1e-10)

Arguments

Value

numeric vector of the same length as x which may have to be thought of recycled along size and prob or mu.

Author(s)

R Core and Martin Maechler

References

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain.

See Also

dbinom_raw; Note that our CRAN package Rmpfr provides dnbinom, dbinom and more, where mpfr-accurate functions are used instead of R's (and our pure R version of) bd0() and stirlerr().

Examples

stopifnot( dnbinomR(0, 1, 1) == 1 )
 size <- 1000 ; x <- 0:size
 dnb <- dnbinomR(x, size, prob = 5/8, log = FALSE, eps = 1e-10)
 plot(x, dnb, type="b")
 all.equal(dnb, dnbinom(x, size, prob = 5/8)) ## mean rel. diff: 0.00017...

 dnbm <- dnbinom.mu(x, size, mu = 123, eps = 1e-10)
 all.equal(dnbm, dnbinom(x, size, mu = 123)) #  Mean relative diff: 0.00069...

Description

dntJKBf1 implements the summation formulas of Johnson, Kotz and Balakrishnan (1995), (31.15) on page 516 and (31.15') on p.519, the latter being typo-corrected for a missing factor $1 / j!$.

dntJKBf() is Vectorize(dntJKBf1, c("x","df","ncp")), i.e., works vectorized in all three main arguments x, df and ncp.

The functions .dntJKBch1() and .dntJKBch() are only there for didactical reasons allowing to check that indeed formula (31.15') in the reference is missing a $j!$ factor in the denominator.

The dntJKBf*() functions are written to also work with arbitrary precise numbers of class "mpfr" (from package Rmpfr) as arguments.

Usage

dntJKBf1(x, df, ncp, log = FALSE, M = 1000)
dntJKBf (x, df, ncp, log = FALSE, M = 1000)

## The "checking" versions, only for proving correctness of formula:
.dntJKBch1(x, df, ncp, log = FALSE, M = 1000, check=FALSE, tol.check = 1e-7)
.dntJKBch (x, df, ncp, log = FALSE, M = 1000, check=FALSE, tol.check = 1e-7)

Arguments

Details

How to choose M optimally has not been investigated yet and is probably also a function of the precision of the first three arguments (see getPrec from Rmpfr).

Note that relatedly, R's source code ‘R/src/nmath/dnt.c’ has claimed from 2003 till 2014 but wrongly that the noncentral t density $f(x, *)$ was

    f(x, df, ncp) =
 	df^(df/2) * exp(-.5*ncp^2) /
 	(sqrt(pi)*gamma(df/2)*(df+x^2)^((df+1)/2)) *
 	   sum_{k=0}^Inf  gamma((df + k + df)/2)*ncp^k / prod(1:k)*(2*x^2/(df+x^2))^(k/2) .

These functions (and this help page) prove that it was wrong.

Value

a number for dntJKBf1() and .dntJKBch1().

a numeric vector of the same length as the maximum of the lengths of x, df, ncp for dntJKBf() and .dntJKBch().

Author(s)

Martin Maechler

References

Johnson, N.L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions Vol~2, 2nd ed.; Wiley; chapter 31, Section 5 Distribution Function, p.514 ff

See Also

R's dt; (an improved version of) Viechtbauer's proposal: dtWV.

Examples

tt <-  seq(0, 10, length.out = 21)
ncp <- seq(0,  6, length.out = 31)
dt3R   <- outer(tt, ncp, dt,     df = 3)
dt3JKB <- outer(tt, ncp, dntJKBf, df = 3)
all.equal(dt3R, dt3JKB) # Lnx(64-b): 51 NA's in dt3R

x <- seq(-1,12, by=1/16)
fx <- dt(x, df=3, ncp=5)
re1 <- 1 - .dntJKBch(x, df=3, ncp=5) / fx ; summary(warnings()) # slow, with warnings
op <- options(warn = 2) # (=> warning == error, for now)
re2 <- 1 -  dntJKBf (x, df=3, ncp=5) / fx  # faster, no warnings
stopifnot(all.equal(re1[!is.na(re1)], re2[!is.na(re1)], tol=1e-6))
head( cbind(x, fx, re1, re2) , 20)
matplot(x, log10(abs(cbind(re1, re2))), type = "o", cex = 1/4)

## One of the numerical problems in "base R"'s non-central t-density:
options(warn = 0) # (factory def.)
x <- 2^seq(-12, 32, by=1/8) ; df <- 1/10
dtm <- cbind(dt(x, df=df,           log=TRUE),
             dt(x, df=df, ncp=df/2, log=TRUE),
             dt(x, df=df, ncp=df,   log=TRUE),
             dt(x, df=df, ncp=df*2, log=TRUE)) #.. quite a few warnings:
summary(warnings())
matplot(x, dtm, type="l", log = "x", xaxt="n",
        main = "dt(x, df=1/10, log=TRUE) central and noncentral")
sfsmisc::eaxis(1)
legend("right", legend=c("", paste0("ncp = df",c("/2","","*2"))),
       lty=1:4, col=1:4, bty="n")

(doExtras <- DPQ:::doExtras()) # TRUE e.g. if interactive()
(ncp <- seq(0, 12, by = if(doExtras) 3/4 else 2))
names(ncp) <- nnMs <- paste0("ncp=", ncp)
tt  <- seq(0, 5, by =  1)
dt3R <- outer(tt, ncp, dt,   df = 3)
if(requireNamespace("Rmpfr")) withAutoprint({
   mt  <- Rmpfr::mpfr(tt , 128)
   mcp <- Rmpfr::mpfr(ncp, 128)
   system.time(
       dt3M <- outer(mt, mcp, dntJKBf, df = 3,
                     M = if(doExtras) 1024 else 256)) # M=1024: 7 sec [10 sec on Winb]
   relE <- Rmpfr::asNumeric(sfsmisc::relErrV(dt3M, dt3R))
   relE[tt != 0, ncp != 0]
})

## all.equal(dt3R, dt3V, tol=0) # 1.2e-12

 # ---- using MPFR high accuracy arithmetic (too slow for routine testing) ---
## no such kink here:
x. <- if(requireNamespace("Rmpfr")) Rmpfr::mpfr(x, 256) else x
system.time(dtJKB <- dntJKBf(x., df=df, ncp=df, log=TRUE)) # 43s, was 21s and only 7s ???
lines(x, dtJKB, col=adjustcolor(3, 1/2), lwd=3)
options(op) # reset to prev.

## Relative Difference / Approximation errors :
plot(x, 1 - dtJKB / dtm[,3], type="l", log="x")
plot(x, 1 - dtJKB / dtm[,3], type="l", log="x", xaxt="n", ylim=c(-1,1)*1e-3); sfsmisc::eaxis(1)
plot(x, 1 - dtJKB / dtm[,3], type="l", log="x", xaxt="n", ylim=c(-1,1)*1e-7); sfsmisc::eaxis(1)
plot(x, abs(1 - dtJKB / dtm[,3]), type="l", log="xy", axes=FALSE, main =
     "dt(*, 1/10, 1/10, log=TRUE) relative approx. error",
     sub= paste("Copyright (C) 2019  Martin Maechler  --- ", R.version.string))
for(j in 1:2) sfsmisc::eaxis(j)

Description

Utility functions for "dpq"-computations, parelling those in R's own C source ‘<Rsource>/src/nmath/dpq.h’, (“dpq” := density–probability–quantile).

Usage

.D_0(log.p) # prob/density == 0   (for log.p=FALSE)
.D_1(log.p) # prob         == 1     "       "

.DT_0(lower.tail, log.p) # == 0  when (lower.tail=TRUE, log.p=FALSE)
.DT_1(lower.tail, log.p) # == 1  when     "                "

.D_Lval(p, lower.tail) # p    {L}ower
.D_Cval(p, lower.tail) # 1-p  {C}omplementary

.D_val (x, log.p)  #     x  in pF(x,..)
.D_qIv (p, log.p)  #     p  in qF(p,..)
.D_exp (x, log.p)  # exp(x)        unless log.p where it's  x
.D_log (p, log.p)  #     p           "      "     "    "   log(p)
.D_Clog(p, log.p)  #   1-p           "      "     "    "   log(1-p) == log1p(-)

.D_LExp(x, log.p)  ## [log](1 - exp(-x))     == log1p(- .D_qIv(x))) even more stable

.DT_val (x, lower.tail, log.p) # := .D_val(.D_Lval(x, lower.tail), log.p) ==    x  in pF
.DT_Cval(x, lower.tail, log.p) # := .D_val(.D_Cval(x, lower.tail), log.p) ==  1-x  in pF

.DT_qIv (p, lower.tail, log.p) # := .D_Lval(.D_qIv(p))	==    p	 in qF
.DT_CIv (p, lower.tail, log.p) # := .D_Cval(.D_qIv(p))	==  1-p  in qF

.DT_exp (x, lower.tail, log.p) #  exp( x )
.DT_Cexp(x, lower.tail, log.p) #  exp(1-x)

.DT_log (p, lower.tail, log.p) # log ( p )  in qF
.DT_Clog(p, lower.tail, log.p) # log (1-p)  in qF
.DT_Log (p, lower.tail)        # log ( p )  in qF(p,..,log.p=TRUE)

Arguments

Value

Typically a numeric vector “as” x or p, respectively.

Author(s)

Martin Maechler

See Also

log1mexp() which is called from .D_LExp() and .DT_Log().

Examples

FT <- c(FALSE, TRUE)
stopifnot(exprs = {
    .D_0(log.p = FALSE) ==    (0)
    .D_0(log.p = TRUE ) == log(0)
    identical(c(1,0), vapply(FT, .D_1, double(1)))
})

## all such functions in package DPQ:
eDPQ <- as.environment("package:DPQ")
ls.str(envir=eDPQ, pattern = "^[.]D", all.names=TRUE)
(nD <- local({ n <- names(eDPQ); n[startsWith(n, ".D")] }))
trimW <- function(ch) sub(" +$","", sub("^ +","", ch))
writeLines(vapply(sort(nD), function(nm) {
    B <- deparse(eDPQ[[nm]])
    sprintf("%31s := %s", trimW(sub("function ", nm, B[[1]])),
            paste(trimW(B[-1]), collapse=" "))
                  }, ""))

do.lowlog <- function(Fn, ...) {
    stopifnot(is.function(Fn),
              all(c("lower.tail", "log.p") %in% names(formals(Fn)))) 
    FT <- c(FALSE, TRUE) ; cFT <- c("F", "T")
    L <- lapply(FT, function(lo) sapply(FT, function(lg) Fn(..., lower.tail=lo, log.p=lg)))
    r <- simplify2array(L)
    `dimnames<-`(r, c(rep(list(NULL), length(dim(r)) - 2L),
                      list(log.p = cFT, lower.tail = cFT)))
}
do.lowlog(.DT_0)
do.lowlog(.DT_1)
do.lowlog(.DT_exp, x = 1/4) ; do.lowlog(.DT_exp, x = 3/4)
do.lowlog(.DT_val, x = 1/4) ; do.lowlog(.DT_val, x = 3/4)
do.lowlog(.DT_Cexp, x = 1/4) ; do.lowlog(.DT_Cexp, x = 3/4)
do.lowlog(.DT_Cval, x = 1/4) ; do.lowlog(.DT_Cval, x = 3/4)
do.lowlog(.DT_Clog, p = (1:3)/4) # w/ warn
do.lowlog(.DT_log,  p = (1:3)/4) # w/ warn
do.lowlog(.DT_qIv,  p = (1:3)/4)

## unfinished: FIXME, the above is *not* really checking
stopifnot(exprs = {

})

Description

Log Gamma derivatives, Psi Gamma functions. dpsifn() is an R interface to the R API function R_dpsifn().

Usage

dpsifn(x, m, deriv1 = 0L, k2 = FALSE)

Arguments

Details

dpsifn() is the underlying “workhorse” of R's own digamma, trigamma and (generalized) psigamma functions.

It is useful, e.g., when several derivatives of $\log\Gamma=$lgamma are desired. It computes and returns length-m sequence $(-1)^{k+1} / \Gamma(k+1) \cdot \psi^{(k)}(x)$ for $k = n, n+1,\ldots, n+m-1$, where $n=$deriv1, and $\psi^{(k)}(x)$ is the k-th derivative of $\psi(x)$, i.e., psigamma(x,k). For more details, see the comments in ‘src/nmath/polygamma.c’.

Value

A numeric $l_x \times m$ matrix (where $l_x=$length(x)) of scaled $\psi^{(k)}(x)$ values. The matrix has attributes

Author(s)

Martin Maechler (R interface); R Core et al., see digamma.

References

See those in psigamma

See Also

digamma, trigamma, psigamma.

Examples

x <- seq(-3.5, 6, by=1/4)
dpx <- dpsifn(x, m = if(getRversion() >= "4.2") 7 else 5)
dpx # in R <= 4.2.1, see that sometimes the 'nz' (under-over-flow count) was uninitialized !!
j <- -1L+seq_len(nrow(dpx)); (fj <- (-1)^(j+1)*gamma(j+1))
## mdpsi <- cbind(di =   digamma(x),      -dpx[1,],
## 	       tri=  trigamma(x),       dpx[2,],
## 	       tetra=psigamma(x,2),  -2*dpx[3,],
## 	       penta=psigamma(x,3),   6*dpx[4,],
## 	       hexa =psigamma(x,4), -24*dpx[5,],
## 	       hepta=psigamma(x,5), 120*dpx[6,],
## 	       octa =psigamma(x,6),-720*dpx[7,])
## cbind(x, ie=attr(dpx,"errorCode"), round(mdpsi, 4))
str(psig <- outer(x, j, psigamma))
dpsi <- t(fj * (`attributes<-`(dpx, list(dim=dim(dpx)))))
if(getRversion() >= "4.2") {
      print( all.equal(psig, dpsi, tol=0) )# -> see 1.185e-16
  stopifnot( all.equal(psig, dpsi, tol=1e-15) )
} else { # R <= 4.1.x; dpsifn(x, ..) *not* ok for x < 0
  i <- x >= 0
      print( all.equal(psig[i,], dpsi[i,], tol=0) )# -> see 1.95e-16
  stopifnot( all.equal(psig[i,], dpsi[i,], tol=1e-15) )
}

Description

Compute the density function $f(x)$ of the t distribution with df degrees of freedom and non-centrality parameter ncp, according to Wolfgang Viechtbauer's proposal in 2002. This is an asymptotic formula for “large” df$= \nu$, or mathematically $\nu \to \infty$.

Usage

dtWV(x, df, ncp = 0, log = FALSE)

Arguments

Details

The formula used is “asymptotic”: Resnikoff and Lieberman (1957), p.1 and p.25ff, proposed to use recursive polynomials for (integer !) degrees of freedom $f = 1,2,\dots, 20$, and then, for df$= f > 20$, use the asymptotic approximation which Wolfgang Viechtbauer proposed as a first version of a non-central t density for R (when dt() did not yet have an ncp argument).

Value

numeric vector of density values, properly recycled in (x, df, ncp).

Author(s)

Wolfgang Viechtbauer (2002) post to R-help (https://stat.ethz.ch/pipermail/r-help/2002-October/026044.html), and Martin Maechler (log argument; tweaks, notably recycling).

References

Resnikoff, George J. and Lieberman, Gerald J. (1957) Tables of the non-central t-distribution; Technical report no. 32 (LIE ONR 32), April 1, 1957; Applied Math. and Stat. Lab., Stanford University. https://statistics.stanford.edu/technical-reports/tables-non-central-t-distribution-density-function-cumulative-distribution

See Also

dt, R's (C level) implementation of the (non-central) t density; dntJKBf, for Johnson et al.'s summation formula approximation.

Examples

tt <- seq(0, 10, len = 21)
ncp <- seq(0, 6, len = 31)
dt3R  <- outer(tt, ncp, dt  , df = 3)
dt3WV <- outer(tt, ncp, dtWV, df = 3)
all.equal(dt3R, dt3WV) # rel.err 0.00063
dt25R  <- outer(tt, ncp, dt  , df = 25)
dt25WV <- outer(tt, ncp, dtWV, df = 25)
all.equal(dt25R, dt25WV) # rel.err 1.1e-5

x <- -10:700
fx  <- dt  (x, df = 22, ncp =100)
lfx <- dt  (x, df = 22, ncp =100, log=TRUE)
lfV <- dtWV(x, df = 22, ncp =100, log=TRUE)

head(lfx, 20) # shows that R's dt(*, log=TRUE) implementation is "quite suboptimal"

## graphics
opa <- par(no.readonly=TRUE)
par(mar=.1+c(5,4,4,3), mgp = c(2, .8,0))
plot(fx ~ x, type="l")
par(new=TRUE) ; cc <- c("red", adjustcolor("orange", 0.4))
plot(lfx ~ x, type = "o", pch=".", col=cc[1], cex=2, ann=FALSE, yaxt="n")
sfsmisc::eaxis(4, col=cc[1], col.axis=cc[1], small.args = list(col=cc[1]))
lines(x, lfV, col=cc[2], lwd=3)
dtt1 <- "      dt"; dtt2 <- "(x, df=22, ncp=100"; dttL <- paste0(dtt2,", log=TRUE)")
legend("right", c(paste0(dtt1,dtt2,")"), paste0(c(dtt1,"dtWV"), dttL)),
       lty=1, lwd=c(1,1,3), col=c("black", cc), bty = "n")
par(opa) # reset

Description

Compute $e^x - 1 - x =$ exp(x) - 1 - x accurately, notably for small $|x|$.

The last two entries in cutx[] denote boundaries where expm1x(x) uses direct formulas. For nC <- length(cutx), exp(x) - 1 - x is used for abs(x) >= cutx[nC], and when abs(x) < cutx[nC] expm1(x) - x is used for abs(x) >= cutx[nC-1].

Usage

expm1x(x, cutx = c( 4.4e-8, 0.1, 0.385, 1.1, 2),
             k = c(2,      9,  12,    17))

expm1xTser(x, k)

Arguments

Value

a vector like x containing (approximations to) $e^x - x - 1$.

Author(s)

Martin Maechler

See Also

expm1(x) for computing $e^x - 1$ is much more widely known, and part of the ISO C standards now.

Examples

## a symmetric set of negative and positive
x <- unique(c(2^-seq(-3/8, 54, by = 1/8), seq(7/8, 3, by = 1/128)))
x <- x0 <- sort(c(-x, 0, x)) # negative *and* positive

## Mathematically,  expm1x() = exp(x) - 1 - x  >= 0  (and == 0 only at x=0):
em1x <- expm1x(x)
stopifnot(em1x >= 0, identical(x == 0, em1x == 0))

plot (x, em1x, type='b', log="y")
lines(x, expm1(x)-x, col = adjustcolor(2, 1/2), lwd = 3) ## should nicely cover ..
lines(x, exp(x)-1-x, col = adjustcolor(4, 1/4), lwd = 5) ## should nicely cover ..
cuts <- c(4.4e-8, 0.10, 0.385, 1.1, 2)[-1] # *not* drawing 4.4e-8
v <- c(-rev(cuts), 0, cuts); stopifnot(!is.unsorted(v))
abline(v = v, lty = 3, col=adjustcolor("gray20", 1/2))

stopifnot(diff(em1x[x <= 0]) <= 0)
stopifnot(diff(em1x[x >= 0]) >= 0)

## direct formula - may be really "bad" :
expm1x.0 <- function(x) exp(x) -1 - x
## less direct formula - improved (but still not universally ok):
expm1x.1 <- function(x) expm1(x)  - x

ax <- abs(x) # ==> show negative and positive x on top of each other
plot (ax, em1x, type='l', log="xy", xlab = "|x|  (for negative and positive x)")
lines(ax, expm1(x)-x, col = adjustcolor(2, 1/2), lwd = 3) ## see problem at very left
lines(ax, exp(x)-1-x, col = adjustcolor(4, 1/4), lwd = 5) ## see huge problems for |x| < ~10^{-7}
legend("topleft", c("expm1x(x)", "expm1(x) - x", "exp(x) - 1 - x"), bty="n",
       col = c(1,2,4), lwd = c(1,3,5))

## -------------------- Relative error of Taylor series approximations :
twoP <- seq(-0.75, 54, by = 1/8)
x <- 2^-twoP
x <- sort(c(-x,x)) # negative *and* positive
e1xAll <- cbind(expm1x.0 = expm1x.0(x),
                expm1x.1 = expm1x.1(x),
                vapply(1:15, \(k) expm1xTser(x, k=k), x))
colnames(e1xAll)[-(1:2)] <- paste0("k=",1:15)
head(e1xAll)
## TODO  plot !!

Description

Format numbers in [0,1] with “precise” result, notably using "1-.." if needed.

Usage

format01prec(x, digits = getOption("digits"), width = digits + 2,
             eps = 1e-06, ...,
             FUN = function(x, ...) formatC(x, flag = "-", ...))

Arguments

Value

a character vector of the same length as x.

Author(s)

Martin Maechler, 14 May 1997

See Also

formatC, format.pval.

Examples

## Show that format01prec()  does reveal more precision :
cbind(format      (1 - 2^-(16:24)),
      format01prec(1 - 2^-(16:24)))

## a bit more variety
e <- c(2^seq(-24,0, by=2), 10^-(7:1))
ee <- sort(unique(c(e, 1-e)))
noquote(ff <- format01prec(ee))
data.frame(ee, format01prec = ff)

Description

Both are R versions of C99 (and POSIX) standard C (and C++) mathlib functions of the same name.

frexp(x) computes base-2 exponent e and “mantissa”, or fraction r, such that $x = r * 2^e$, where $r \in [0.5, 1)$ (unless when x is in c(0, -Inf, Inf, NaN) where r == x and e is 0), and $e$ is integer valued.

ldexp(f, E) is the inverse of frexp(): Given fraction or mantissa f and integer exponent E, it returns $x = f * 2^E$. Viewed differently, it's the fastest way to multiply or divide (double precision) numbers with $2^E$.

Usage

frexp(x)
ldexp(f, E)