Nothing
require("stabledist")
pPareto <- stabledist:::pPareto
source(system.file("test-tools-1.R", package = "Matrix"), keep.source=interactive())
#-> identical3(), showProc.time(),...
(doExtras <- stabledist:::doExtras())
options(pstable.debug = FALSE)
options(pstable.debug = TRUE)# want to see when uniroot() is called
stopifnot(all.equal(pstable(0.3, 0.75, -.5, tol= 1e-14),
0.6688726496, tol = 1e-8))
## was 0.66887227658457, tol = 1e-10))
pstable(-4.5, alpha = 1, beta = 0.01)## gave integration error (now uniroot..)
## a "outer vectorized" version:
pstabALL <- function(x, alpha, beta, ...)
sapply(alpha, function(alph)
sapply(beta, function(bet)
pstable(x, alph, bet, ...)))
alph.s <- (1:32)/16 # in (0, 2]
beta.s <- (-16:16)/16 # in [-1, 1]
stopifnot(pstabALL( Inf, alph.s, beta.s) == 1,
pstabALL(-Inf, alph.s, beta.s, log.p=TRUE) == -Inf,
pstabALL( 0, alph.s, beta = 0) == 0.5,
TRUE)
pdf("pstab-ex.pdf")
##---- log-scale -------------
r <- curve(pstable(x, alpha=1.8, beta=.9,
lower.tail=FALSE, log.p=TRUE),
5, 150, n=500,
log="x",type="b", cex=.5)
curve(pPareto(x, alpha=1.8, beta=.9,
lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2)
##--> clearly potential for improvement!
## the less extreme part - of that:
r <- curve(pstable(x, alpha=1.8, beta=.9,
lower.tail=FALSE, log.p=TRUE),
1, 50, n=500, log="x")
curve(pPareto(x, alpha=1.8, beta=.9, lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2)
## Check that pstable() is the integral of dstable() --- using simple Simpson's rule
## in it's composite form:
## \int_a^b f(x) dx\approx \frac{h}{3}
## \bigg[ f(x_0) + 2 \sum_{j=1}^{n/2-1}f(x_{2j}) +
## + 4 \sum_{j=1}^{n/2} f(x_{2j-1}) +
## + f(x_n) \bigg],
intSimps <- function(fx, h) {
stopifnot((n <- length(fx)) %% 2 == 0,
n >= 4, length(h) == 1, h > 0)
n2 <- n %/% 2
j2 <- 2L * seq_len(n2-1)
j4 <- 2L * seq_len(n2) - 1L
h/3 * sum(fx[1], 2* fx[j2], 4* fx[j4], fx[n])
}
chk.pd.stable <- function(alpha, beta, xmin=NA, xmax=NA,
n = 256, do.plot=TRUE,
comp.tol = 1e-13, eq.tol = 1e-3)
{
stopifnot(n >= 20)
if(is.na(xmin)) xmin <- qstable(0.01, alpha, beta)
if(is.na(xmax)) xmax <- qstable(0.99, alpha, beta)
dx <- ceiling(1024*grDevices::extendrange(r = c(xmin, xmax), f = 0.01))/1024
h <- diff(dx)/n
x <- seq(dx[1], dx[2], by = h)
fx <- dstable(x, alpha=alpha, beta=beta, tol= comp.tol)
Fx <- pstable(x, alpha=alpha, beta=beta, tol=2*comp.tol)
i.ev <- (i <- seq_along(x))[i %% 2 == 0 & i >= max(n/10, 16)]
## integrate from x[1] up to x[i] (where i is even);
## the exact value will be F(x[i]) - F(x[1]) == Fx[i] - Fx[1]
Fx. <- vapply(lapply(i.ev, seq_len),
function(ii) intSimps(fx[ii], h), 0)
a.eq <- all.equal(Fx., Fx[i.ev] - Fx[1], tol = eq.tol)
if(do.plot) {
## Show the fit
plot(x, Fx - Fx[1], type = "l")
lines(x[i.ev], Fx., col=adjustcolor("red", 0.5), lwd=3)
op <- par(ask=TRUE) ; on.exit(par(op))
## show the "residual", i.e., the relative error
plot(x[i.ev], 1- Fx./(Fx[i.ev] - Fx[1]),
type = "l", xlim = range(x))
abline(h=0, lty=3, lwd = .6)
}
if(!isTRUE(a.eq)) stop(a.eq)
invisible(list(x=x, f=fx, F=Fx, i. = i.ev, F.appr. = Fx.))
}
op <- par(mfrow=2:1, mar = .1+c(3,3,1,1), mgp=c(1.5, 0.6,0))
c1 <- chk.pd.stable(.75, -.5, -1, 1.5, eq.tol = .006)
c2 <- chk.pd.stable(.95, +0.6, -1, 1.5, eq.tol = .006)# with >= 50 warnings
## here are the "values"
with(c1, all.equal(F.appr., F[i.] - F[1], tol = 0)) # (.0041290 on 64-Lnx)
with(c2, all.equal(F.appr., F[i.] - F[1], tol = 0)) # (.0049307 on 64-Lnx)
showProc.time() #
c3 <- chk.pd.stable(.95, +0.9, -3, 15) # >= 50 warnings
curve(dstable(x, .999, -0.9), -20, 5, log="y")
curve(pstable(x, .999, -0.9), -20, 5, log="y")#-> using uniroot
c4 <- chk.pd.stable(.999, -0.9, -20, 5)
showProc.time() #
## alpha == 1 , small beta ---- now perfect
curve(pstable(x, alpha=1, beta= .01), -6, 8, ylim=0:1)
abline(h=0:1, v=0, lty=3, col="gray30")
n <- length(x <- seq(-6,8, by = 1/16))
px <- pstable(x, alpha=1, beta= .01)
## now take approximation derivative by difference ratio:
x. <- (x[-n]+x[-1])/2
plot (x., diff(px)*16, type="l")
## now check convexity/concavity :
f2 <- diff(diff(px))
stopifnot(f2[x[-c(1,n)] < 0] > 0,
f2[x[-c(1,n)] > 0] < 0)
## and compare with dstable() ... which actually shows dstable() problem:
fx. <- dstable(x., alpha=1, beta=.01)
lines(x., fx., col = 2, lwd=3, lty="5111")
if(dev.interactive(orNone=TRUE)) {
curve(dstable(x, 1., 0.99), -6, 50, log="y")# "uneven" (x < 0); 50 warnings
curve(dstable(x, 1.001, 0.95), -10, 30, log="y")# much better
}
showProc.time() #
if(doExtras) {
c5 <- chk.pd.stable(1., 0.99, -6, 50)# -> uniroot
c6 <- chk.pd.stable(1.001, 0.95, -10, 30)# -> uniroot; 2nd plot *clearly* shows problem
with(c5, all.equal(F.appr., F[i.] - F[1], tol = 0)) # .00058 on 64-Lnx
with(c6, all.equal(F.appr., F[i.] - F[1], tol = 0)) # 2.611e-5 on 64-Lnx
## right tail:
try(## FIXME:
c1.0 <- chk.pd.stable(1., 0.8, -6, 500)# uniroot; rel.difference = .030
)
## show it more clearly
curve(pstable(x, alpha=1, beta=0.5), 20, 800, log="x", ylim=c(.97, 1))
curve(pPareto(x, alpha=1, beta=0.5), add=TRUE, col=2, lty=2)
abline(h=1, lty=3,col="gray")
# and similarly (alpha ~= 1, instead of == 1):
curve(pstable(x, alpha=1.001, beta=0.5), 20, 800, log="x", ylim=c(.97, 1))
curve(pPareto(x, alpha=1.001, beta=0.5), add=TRUE, col=2, lty=2)
abline(h=1, lty=3,col="gray")
## zoom
curve(pstable(x, alpha=1.001, beta=0.5), 100, 200, log="x")
curve(pPareto(x, alpha=1.001, beta=0.5), add=TRUE, col=2, lty=2)
## but alpha = 1 is only somewhat better as approximation:
curve(pstable(x, alpha=1 , beta=0.5), add=TRUE, col=3,
lwd=3, lty="5131")
showProc.time() #
}
c7 <- chk.pd.stable(1.2, -0.2, -40, 30)
c8 <- chk.pd.stable(1.5, -0.999, -40, 30)# two warnings
showProc.time() #
### Newly found -- Marius Hofert, 18.Sept. (via qstable):
stopifnot(all.equal(qstable(0.6, alpha = 0.5, beta = 1,
tol=1e-15, integ.tol=1e-15),
2.636426573120147))
##--> which can be traced to the first of
stopifnot(pstable(q= -1.1, alpha=0.5, beta=1) == 0,
pstable(q= -2.1, alpha=0.6, beta=1) == 0)
## Found by Tobias Hudec, 2 May 2015:
stopifnot(
all.equal(1.5, qstable(p=0.5, alpha=1.5, beta=0, gamma=2, delta = 1.5)),
all.equal(1.5, qstable(p=0.5, alpha=0.6, beta=0, gamma=0.2, delta = 1.5))
)
## Stable(alpha = 1/2, beta = 1, gamma, delta, pm = 1) <===> Levy(delta, gamma)
source(system.file("xtraR", "Levy.R", package = "stabledist"), keep.source=interactive())
##-> dLevy(x, mu, c, log) and
##-> pLevy(x, mu, c, log.p, lower.tail)
set.seed(101)
show.Acc <- (interactive() && require("Rmpfr"))
if(show.Acc) { ## want to see accuracies, do not stop "quickly"
format.relErr <- function(tt, cc)
format(as.numeric(relErr(tt, cc)), digits = 4, width = 8)
}
## FIXME: Look why pstable() is so much less accurate than dstable()
## even though the integration in dstable() is more delicate in general
pTOL <- 1e-6 # typically see relErr of 5e-7
dTOL <- 1e-14 # typically see relErr of 1.3...3.9 e-15
showProc.time()
## Note that dstable() is more costly than pstable()
for(ii in 1:(if(doExtras) 32 else 8)) {
Z <- rnorm(2)
mu <- sign(Z[1])*exp(Z[1])
sc <- exp(Z[2])
x <- seq(mu, mu+ sc* 100*rchisq(1, df=3),
length.out= if(doExtras) 512 else 32)
## dLevy() and pLevy() using only pnorm() are "mpfr-aware":
x. <- if(show.Acc) mpfr(x, 256) else x
pL <- pLevy(x., mu, sc)
pS <- pstable(x, alpha=1/2, beta=1, gamma=sc, delta=mu,
pm = 1)
dL <- dLevy(x., mu, sc)
dS <- dstable(x, alpha=1/2, beta=1, gamma=sc, delta=mu,
pm = 1)
if(show.Acc) {
cat("p: ", format.relErr(pL, pS), "\t")
cat("d: ", format.relErr(dL, dS), "\n")
} else {
cat(ii %% 10)
}
stopifnot(all.equal(pL, pS, tol = pTOL),
all.equal(dL, dS, tol = dTOL))
}; cat("\n")
showProc.time()## ~ 75 sec (doExtras=TRUE) on lynne (2012-09)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.