Skip to content

The Stoppa distribution

3 messages · William Dunlap, Andrew Hoerner

William Dunlap
# 
Try the following, which have the usual lower.tail and log.p arguments to make
it easier to get accurate results in the tails.  logspace_sub() is an R version of
the R C-API function Rf_logspace_sub().  I haven't tested the [pq]stoppa functions
much.

logspace_sub <- function (logx, logy)
{
    # log(exp(logx) - exp(logy)), avoiding unnecessary floating point error
    dxy <- logx - logy
    # log(2) looks like best breakpoint
    logx + ifelse(dxy < 0.693147180559945, log(-expm1(-dxy)), log1p(-exp(-dxy)))
}

pstoppa <- function(q, y0, alpha, theta = 1, lower.tail = TRUE, log.p = FALSE)
{
    logp <- theta * logspace_sub(0, -alpha * log(q/y0))
    if (!lower.tail) {
        logp <- logspace_sub(0, logp)
    }
    if (log.p) logp else exp(logp)
}

qstoppa <- function(p, y0, alpha, theta = 1, lower.tail = TRUE, log.p = FALSE)
{
    logp <- if (log.p) {
                if (lower.tail) p else logspace_sub(0, p)
            } else {
                if (lower.tail) log(p) else log1p(-p)
            }
    logq <- log(y0) - 1/alpha * logspace_sub(0, logp/theta)
    exp(logq)
}

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
2 days later