# R/pDist.R In DistributionUtils: Distribution Utilities

#### Documented in pDistqDist

```pDist <- function(densFn = "norm", q, param = NULL,
subdivisions = 100, lower.tail = TRUE,
intTol = .Machine\$double.eps^0.25,
valueOnly = TRUE, ...)
{
## CALL <- match.call()
dfun <- match.fun(paste("d", densFn, sep = ""))
mode <- distMode(densFn, param = param, ...)
## match the density function for different distributions
qLess <- which((q <= mode)&(is.finite(q)))
## when q is less than mode
qGreater <- which((q > mode)&(is.finite(q)))
## when q is greater than mode
prob <- rep(NA, length(q))
err <- rep(NA, length(q))
prob[q == -Inf] <- 0
prob[q == Inf] <- 0
err[q %in% c(-Inf, Inf)] <- 0

## Integrate density from Inf / -Inf to the mode
for (i in qLess){
if(is.null(param)){
intRes <- integrate(dfun, -Inf, q[i],
subdivisions = subdivisions,
rel.tol = intTol, ...)
} else {
intRes <- integrate(dfun, -Inf, q[i], param = param,
subdivisions = subdivisions,
rel.tol = intTol, ...)
}
prob[i] <- intRes\$value
err[i] <- intRes\$abs.error
}

for (i in qGreater){
if(is.null(param)){
intRes <- integrate(dfun, q[i], Inf,
subdivisions = subdivisions,
rel.tol = intTol, ...)
} else {
intRes <- integrate(dfun, q[i], Inf, param = param,
subdivisions = subdivisions,
rel.tol = intTol, ...)
}
prob[i] <- intRes\$value
err[i] <- intRes\$abs.error
}

if (lower.tail)
prob[q > mode] <- 1 - prob[q > mode]
else
prob[q <= mode] <- 1 - prob[q <= mode]
# Return:
if(valueOnly)
prob
else
list(value = prob, error = err)
}

qDist <- function(densFn = "norm", p, param = NULL,
lower.tail = TRUE, method = "spline",  nInterpol = 501,
uniTol = .Machine\$double.eps^0.25, subdivisions = 100,
intTol = uniTol, ...)
{
## CALL <- match.call()
stopifnot(0 <= p, p <= 1)
if (!lower.tail) {
p <- 1 - p
lower.tail <- TRUE
}
mode <- distMode(densFn, param = param, ...)
pMode <- pDist(densFn, q = mode, param = param, intTol = intTol, ...)
quant <- rep(NA, length(p))
delta <- if(densFn == "skewhyp") {
if(is.null(param)) list(...)\$delta else param[2]
} else 0
xRange <- distCalcRange(densFn, param = param, tol = 10^(-5), ...)

if (method == "integrate"){
less <- which((p <= pMode) & (p >.Machine\$double.eps^7.5))
quant <- ifelse(p <= .Machine\$double.eps^5, -Inf, quant)
if (length(less) > 0){
pLow <- min(p[less])
step <- distStepSize(densFn, param = param,
dist = delta, side = "left", ...)[1]
xLow <- mode - step
if (densFn == "skewhyp"){
while(pDist(densFn, xLow,
param = param, intTol = intTol, ...) >= pLow)
{
xLow <- xLow -
distStepSize(densFn, param = param,
dist = delta, side = "left", ...)[1]
}
} else {
while(pDist(densFn, xLow,
param = param, intTol = intTol, ...) >= pLow)
{
xLow <- xLow - step
}
}
xRange <- c(xLow, mode)
zeroFn <- function(x, p)
pDist(densFn, x, param = param, intTol = intTol, ...) - p
for (i in less) {
quant[i] <- uniroot(zeroFn, p = p[i],
interval = xRange, tol = uniTol)\$root
}
}

greater <- which ((p > pMode) & (p < (1 - .Machine\$double.eps^7.5)))
p[greater] <- 1 - p[greater]
quant <- ifelse(p >= (1 - .Machine\$double.eps^5), Inf, quant)
if (length(greater) > 0){
pHigh <- min(p[greater])
step <-  distStepSize(densFn, param = param,
dist = delta, side = "right", ...)[1]
xHigh <- mode + step
if (densFn == "skewhyp"){
while (pDist(densFn, xHigh, param = param, intTol = intTol,
lower.tail = FALSE, ...) >= pHigh)
{
xHigh <- xHigh +
distStepSize(densFn, param = param,
dist = delta, side = "left", ...)[1]
}
} else {
while (pDist(densFn, xHigh, param = param, intTol = intTol,
lower.tail = FALSE, ...) >= pHigh)
{
xHigh <- xHigh + step
}
}
xRange <- c(mode, xHigh)
zeroFn <- function(x, p){
return(pDist(densFn, x, param = param, intTol = intTol,
lower.tail = FALSE, ...) - p)
}
for (i in greater)
{
quant[i] <- uniroot(zeroFn, p = p[i],
interval = xRange, tol = uniTol)\$root
}
}
} else if (method == "spline"){
inRange <- which((p > pDist(densFn, xRange[1],
param = param, intTol = intTol, ...)) &
(p < pDist(densFn, xRange[2],
param = param, intTol = intTol, ...)))
small <- which((p <= pDist(densFn, xRange[1],
param = param, intTol = intTol, ...)) & (p >= 0))
large <- which((p >= pDist(densFn, xRange[2],
param = param, intTol = intTol, ...)) & (p <= 1))
extreme <- c(small, large)
xVals <- seq(xRange[1], xRange[2], length.out = nInterpol)
yVals <- pDist(densFn, xVals, param = param,
subdivisions = subdivisions, intTol = intTol, ...)
splineFit <- splinefun(xVals, yVals)
for (i in inRange){
quant[i] <- uniroot(function(x) splineFit(x) - p[i],
interval = xRange, tol = uniTol)\$root
}

if (length(extreme) > 0){
quant[extreme] <- qDist(densFn, p[extreme], param = param,
lower.tail = lower.tail,
method = "integrate",
nInterpol = nInterpol, uniTol = uniTol,
subdivisions = subdivisions,
intTol = intTol, ...)
}
}
return(quant)
}
```

## Try the DistributionUtils package in your browser

Any scripts or data that you put into this service are public.

DistributionUtils documentation built on May 2, 2019, 4:46 p.m.