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-boun...@r-project.org [mailto:r-help-boun...@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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
