pnchi1sq {DPQ} | R Documentation |
Computes probabilities for the non-central chi-squared distribution, in
special cases, currently for df = 1
and df = 3
, using
‘exact’ formulas only involving the standard normal (Gaussian)
cdf Φ() and its derivative φ(), i.e., R's
pnorm()
and dnorm()
.
pnchi1sq(q, ncp = 0, lower.tail = TRUE, log.p = FALSE, epsS = .01) pnchi3sq(q, ncp = 0, lower.tail = TRUE, log.p = FALSE, epsS = .04)
q |
number ( ‘quantile’, i.e., abscissa value.) |
ncp |
non-centrality parameter delta; .... |
lower.tail, log.p |
logical, see, e.g., |
epsS |
small number, determining where to switch from the
“small case” to the regular case, namely by defining
|
In the “small case” (epsS
above), the direct formulas
suffer from cancellation, and we use Taylor series expansions in
s := √{q}, which in turn use
“probabilists'” Hermite polynomials He_n(x).
The default values epsS
have currently been determined by
experiments as those in the ‘Examples’ below.
a numeric vector “like” q+ncp
, i.e., recycled to common length.
Martin Maechler, notably the Taylor approximations in the “small” cases.
Johnson et al.(1995), see ‘References’ in
pnchisqPearson
.
https://en.wikipedia.org/wiki/Hermite_polynomials
pchisq
, the (simple and R-like) approximations, such as
pnchisqPearson
and the wienergerm approximations,
pchisqW()
etc.
qq <- seq(9500, 10500, length=1000) m1 <- cbind(pch = pchisq (qq, df=1, ncp = 10000), p1 = pnchi1sq(qq, ncp = 10000)) matplot(qq, m1, type = "l"); abline(h=0:1, v=10000+1, lty=3) all.equal(m1[,"p1"], m1[,"pch"], tol=0) # for now, 2.37e-12 m3 <- cbind(pch = pchisq (qq, df=3, ncp = 10000), p3 = pnchi3sq(qq, ncp = 10000)) matplot(qq, m3, type = "l"); abline(h=0:1, v=10000+3, lty=3) all.equal(m3[,"p3"], m3[,"pch"], tol=0) # for now, 1.88e-12 stopifnot(exprs = { all.equal(m1[,"p1"], m1[,"pch"], tol=1e-10) all.equal(m3[,"p3"], m3[,"pch"], tol=1e-10) }) ### Very small 'x' i.e., 'q' would lead to cancellation: ----------- ## df = 1 -------------- qS <- c(0, 2^seq(-40,4, by=1/16)) m1s <- cbind(pch = pchisq (qS, df=1, ncp = 1) , p1.0= pnchi1sq(qS, ncp = 1, epsS = 0) , p1.4= pnchi1sq(qS, ncp = 1, epsS = 1e-4) , p1.3= pnchi1sq(qS, ncp = 1, epsS = 1e-3) , p1.2= pnchi1sq(qS, ncp = 1, epsS = 1e-2) ) cols <- adjustcolor(1:5, 1/2); lws <- seq(4,2, by = -1/2) abl.leg <- function(x.leg = "topright", epsS = 10^-(4:2), legend = NULL) { abline(h = .Machine$double.eps, v = epsS^2, lty = c(2,3,3,3), col= adjustcolor(1, 1/2)) if(is.null(legend)) legend <- c(quote(epsS == 0), as.expression(lapply(epsS, function(K) substitute(epsS == KK, list(KK = formatC(K, w=1)))))) legend(x.leg, legend, lty=1:4, col=cols, lwd=lws, bty="n") } matplot(qS, m1s, type = "l", log="y" , col=cols, lwd=lws) matplot(qS, m1s, type = "l", log="xy", col=cols, lwd=lws) ; abl.leg("right") ## ==== "Errors" =================================================== ## Absolute: ------------------------- matplot(qS, m1s[,1] - m1s[,-1] , type = "l", log="x" , col=cols, lwd=lws) matplot(qS, abs(m1s[,1] - m1s[,-1]), type = "l", log="xy", col=cols, lwd=lws) abl.leg("bottomright") ## Relative: ------------------------- matplot(qS, 1 - m1s[,-1]/m1s[,1] , type = "l", log="x", col=cols, lwd=lws) abl.leg() matplot(qS, abs(1 - m1s[,-1]/m1s[,1]), type = "l", log="xy", col=cols, lwd=lws) abl.leg() ## df = 3 -------------- %% FIXME: the 'small' case is clearly wrong <<< qS <- c(0, 2^seq(-40,4, by=1/16)) ee <- c(1e-3, 1e-2, .04) m3s <- cbind(pch = pchisq (qS, df=3, ncp = 1) , p1.0= pnchi3sq(qS, ncp = 1, epsS = 0) , p1.3= pnchi3sq(qS, ncp = 1, epsS = ee[1]) , p1.2= pnchi3sq(qS, ncp = 1, epsS = ee[2]) , p1.1= pnchi3sq(qS, ncp = 1, epsS = ee[3]) ) matplot(qS, m3s, type = "l", log="y" , col=cols, lwd=lws) matplot(qS, m3s, type = "l", log="xy", col=cols, lwd=lws); abl.leg("right", ee) ## ==== "Errors" =================================================== ## Absolute: ------------------------- matplot(qS, m3s[,1] - m3s[,-1] , type = "l", log="x" , col=cols, lwd=lws) matplot(qS, abs(m3s[,1] - m3s[,-1]), type = "l", log="xy", col=cols, lwd=lws) abl.leg("right", ee) ## Relative: ------------------------- matplot(qS, 1 - m3s[,-1]/m3s[,1] , type = "l", log="x", col=cols, lwd=lws) abl.leg(, ee) matplot(qS, abs(1 - m3s[,-1]/m3s[,1]), type = "l", log="xy", col=cols, lwd=lws) abl.leg(, ee)