[R] The Stoppa distribution

[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


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of Andrew Hoerner
> Sent: Friday, December 13, 2013 4:13 PM
> To: r-help
> Subject: [R] The Stoppa distribution
> 
> The Stoppa distribution is a 3-parameter distribution that generalizes the
> Pareto distribution, adding a second shape parameter but no location term.
> The CDF is
> 
>                                     F(x) =  [1-(x/x0)-á]è    0 < x0 < x
> 
> Kleiber & Kotz (2003). *Statistical Size Distributions in Economics and
> Actuarial Sciences.*    I do not believe that the Stoppa distribution has
> been implemented in R under that name.
> 
> 
> Does anyone know if the Stoppa distribution has been implimented in R under
> a different name, or if there is a more generalized distribution containing
> the Stoppa as a special case which has been implemented in R? (The
> generalized beta of the second kind, maybe?) And if so, how you need to set
> the parameters of the more general distribution to collapse it to the
> Stoppa?
> 
> Any help anyone could offer would be greatly appreciated.
> 
> Sincerely, andrewH
> 
> 	[[alternative HTML version deleted]]