Nothing
require("stabledist")
dPareto <- stabledist:::dPareto
source(system.file("test-tools-1.R", package = "Matrix"), keep.source=interactive())
#-> identical3(), showProc.time(),...
(doExtras <- stabledist:::doExtras())
if(!require("sfsmisc")) eaxis <- axis # use sfsmisc::eaxis if available
stopifnot(0 < print(dstable(4000., alpha=1.00001, beta=0.6)))
## gave error in fBasics::dstable()
## now 18 warnings from uniroot()
pdf("dstab-ex.pdf")
x <- 2^seq(0, 20, length= if(doExtras) 200 else 64)
## Regression check for alpha=2: <==> *norm() :
x. <- x/1024
fx <- dstable(x., alpha = 2, beta = 0, gamma = 1.1, delta=2, pm=2)
lf <- dstable(x., alpha = 2, beta = 0, gamma = 1.1, delta=2, pm=2, log=TRUE)
stopifnot(
local({i <- is.finite(log(fx)); all.equal(log(fx[i]), lf[i])}),
all.equal(fx, dnorm(x., m=2, s=1.1)),
all.equal(lf, dnorm(x., m=2, s=1.1, log=TRUE)))
fx <- dstable(x, alpha = 1.0001, beta = 0.6)
plot(x,fx, log="x", type="l")# looks good
plot(x,fx, log="xy", type="l")# --- perfect now
stopifnot((dlx <- diff(log(fx))) < 0, # decreasing
{ dl <- dlx[x[-1] > 1000] # the negative slope:
relErr(if(doExtras) -0.13934 else -0.44016, dl) < 4e-4},
diff(dl) > -1e-6)# > 0 "in theory"; > -6e-7, ok on 64-bit Linux
nc <- if(doExtras) 512 else 101 # number of points for curve()
zeta <- function(alpha,beta) if(alpha==1) 0 else -beta*tan(pi/2*alpha)
## negative beta:
cx <- curve(dstable(x, 0.75, -.5), -.5, 1.5, n=nc)# ok, now
m <- stableMode(0.75, -.5, tol=1e-14)
stopifnot(all.equal(m, 0.35810298366, tol = 1e-7),
all.equal(dstable(m, 0.75, -.5), 0.3554664043, tol=1e-6))
showProc.time()
###-------- "small" alpha -----------------
## alpha --> 0 --- very heavy tailed -- and numerically challenging.
## symmetric (beta = 0)
(x0 <- (-16:16)/256)
fx0 <- dstable(x0, alpha = 0.1, beta=0, gamma = 1e6)
plot(x0, fx0, type = "o",
main = expression(f(x, alpha== 0.1, beta == 0, gamma == 10^6)))
stopifnot(all.equal(fx0[17],1.15508291498374),
all.equal(fx0[ 1],0.02910420736536),
all.equal(range(diff(fx0[1:8])),
c(0.0011871409, 0.0025179435), tol=1e-6)
)
## beta > 0
r3 <- curve(dstable(x, alpha = 0.3, beta = 0.5, tol=1e-7),
-1, 1)
m3 <- stableMode(0.3, 0.5, tol=1e-14)# still with 3 warnings
stopifnot(all.equal(m3, -0.2505743952946, tol = 1e-10))
## zoom into the above
r3. <- curve(dstable(x, alpha = 0.3, beta = 0.5, tol=1e-7),
-.27, -.22)
abline(v = m3, col="gray40", lty=2)
r1 <- curve(dstable(x, alpha = 0.1, beta = 0.5, tol=1e-7),
-.4, .2, ylim = c(0, 20), n = nc)
m1 <- stableMode(0.1, 0.5, tol=1e-15)
abline(v=m1, h=0, col="gray40", lty=2)
stopifnot(all.equal(m1, -0.079193, tol=1e-5)) # -0.079193188, was -0.079192219, and -0.0791921758
title(main = expression(f(x, alpha== 0.1, beta == 0.5)))
## check mode *and* unimodality
i. <- r1$x > m1
stopifnot(## decreasing to the right:
diff(r1$y[ i.]) < 0,
## increasing on the left:
diff(r1$y[!i.]) > 0)
## Now *really* small alpha --> 0:
## --------------------------
## Note that
if(require("FMStable")) {
try( FMStable::setParam(alpha = 1e-18, location=0, logscale=0, pm = 0) )
## gives
## Error in .... setParam: S0=M parametrization not suitable for small alpha
}
## now if this is true (and I think we can trust Geoff Robinson on this),
## we currently have a "design bug - problem":
## as internally, we always translate to 'pm=0' parametrization __FIXME_(how ??)__
## --> solution: see below: there's a simple (alpha -> 0) asymptotic
## These now "work": .... well with integration warnings !!
pdstab.alpha <- function(x, beta, alphas = 2^-(40:2),
type = "o", add=FALSE, log = "xy", ...)
{
stopifnot(is.numeric(x), length(x) == 1, length(beta) == 1,
is.numeric(beta), -1 <= beta, beta <= 1,
is.numeric(alphas), 0 <= alphas, alphas <= 2,
is.logical(add))
if(is.unsorted(alphas)) alphas <- sort(alphas)
dst.alp <- vapply(alphas, function(aaa)
dstable(x, aaa, beta = beta, pm = 0), 1.) ## with warnings
if(add)
lines(alphas, dst.alp, type=type, ...)
else {
plot(alphas, dst.alp, type=type, log=log, axes=FALSE,
xlab = quote(alpha),
ylab = bquote(f(x == .(x), alpha)),
main = bquote(dstable(x == .(x), alpha, beta == .(beta), pm == 0)
~~~"for (very) small" ~ alpha))
if(!require("sfsmisc")) eaxis <- axis # use sfsmisc::eaxis if available
eaxis(1); eaxis(2)
}
invisible(cbind(alphas, dstable = dst.alp))
}
## vary beta, keep x :
pdstab.alpha(x = -1, beta = 0.1)
if(doExtras) {
pdstab.alpha(x = -1, beta = 0.3, add=TRUE, col=2, type="l")
pdstab.alpha(x = -1, beta = 0.5, add=TRUE, col=3, type="l")
pdstab.alpha(x = -1, beta = 0.7, add=TRUE, col=4, type="l")
pdstab.alpha(x = -1, beta = 0.9, add=TRUE, col=5, type="l")
## vary x, keep beta
pdstab.alpha(x = -1, beta = 0.3)
pdstab.alpha(x = +1, beta = 0.3, add=TRUE, col=2, type="l")
pdstab.alpha(x = +5, beta = 0.3, add=TRUE, col=3, type="l")
pdstab.alpha(x = +50, beta = 0.3, add=TRUE, col=4, type="l")
pdstab.alpha(x = -10, beta = 0.3, add=TRUE, col=5, type="l")
}
## Plots suggest a simple formula
## log(f(x, alpha, beta))= c(x,beta) + log(alpha)
## f(x, alpha, beta) = C(x,beta) * (alpha) -- ???
## for normal alpha, it looks different {which is reassuring!}
pdstab.alpha(x = -1, beta = 0.3, alphas = seq(1/128, 2, length=100))
## Show the formula, we derived:
## f(x, \alpha, \beta) ~ \alpha * \frac{1 + \beta}{2e \abs{x + \pi/2 \alpha\beta}}
## ONLY ok, when x > zeta := - beta * tan(pi/2 *alpha)
## otherwise, "swap sign" of (x, beta)
dst.smlA <- function(x, alpha, beta, log = FALSE) {
pa <- pi/2*alpha
i. <- x < -pa*beta
if(any(i.)) {
beta <- rep(beta, length.out=length(x))
beta[i.] <- -beta[i.]
x [i.] <- -x [i.]
}
## alpha*(1+beta) / (2*exp(1)*(x+ pa*beta))
r <- log(alpha) + log1p(beta) - (1 + log(2*(x+ pa*beta)))
if(log) r else exp(r)
}
set.seed(17)
alpha <- 1e-15 ## must be larger than cutoff in .fct1() in ../R/dpqr-stable.R :
stopifnot(alpha > stabledist:::.alpha.small.dstable)
for(beta in runif(if(doExtras) 20 else 2, -1, 1)) {
cat(sprintf("beta =%10g: ", beta))
for(N in 1:(if(doExtras) 10 else 1)) {
x <- 10*rnorm(100); cat(".")
stopifnot(all.equal(dstable (x, alpha, beta),
dst.smlA(x, alpha, beta), tol = 1e-13))
}; cat("\n")
}
curve( dstable (x, alpha=1e-10, beta=.8, log=TRUE) , -10, 10)
curve( dst.smlA(x, alpha=1e-10, beta=.8, log=TRUE), add=TRUE, col=2)
## "same" picture (self-similar !)
curve( dstable (x, alpha=1e-10, beta=.8, log=TRUE), -1, 1)
curve( dst.smlA(x, alpha=1e-10, beta=.8, log=TRUE), add=TRUE, col=2)
## Testing stableMode here:
### beta = 1 (extremely skewed) and small alpha: ---------
## --------
## Problem at *left* ("less problematic") tail, namely very close to where the
## support of the density becomes mathematically exactly zero :
##
## clear inaccuracy / bug --- maybe *seems* curable
##
curve(dstable(exp(x), alpha= 0.1, beta=1, pm=1, log=TRUE), -40, 10)
## ------
## --> warnings both from uniroot ("-Inf") and .integrate2()
## about equivalent to
curve(dstable(x, alpha= 0.1, beta=1, pm=1, log=TRUE), 1e-15, 4e4,
log="x", xaxt="n"); eaxis(1)
## If we decrease zeta.tol "to 0", we get better here:
curve(dstable(exp(x), alpha= 0.1, beta=1, pm=1, log=TRUE, zeta.tol=1e-100), -40, 20)
## or here, ... but still not good enough
curve(dstable(exp(x), alpha= 0.1, beta=1, pm=1, log=TRUE, zeta.tol=1e-200), -45, 30)
showProc.time()
##------ NB: Pareto tail behavior -- see more in ./tails.R
## =======
## alpha ~= 1 ---- and x ~ zeta(a,b) -----------
## ==========
f1 <- dstable(6366.197, alpha= 1.00001, beta= .1)
f2 <- dstable(-50929.58, alpha= 1.00001, beta= -.8)
stopifnot(f1 > 0, f2 > 0)
## these all work (luck):
zet <- zeta(alpha= 1.00001, beta= -.8)# -50929.58
## here, we must have larger zeta.tol: = 5e-5 is fine; now using "automatic" default
r4 <- curve(dstable(zet+x, alpha= 1.00001, beta= -.8), -3, 7,
xlab = expression(zeta(alpha,beta) - x),
ylim=c(2.1, 2.3)*1e-10, n = nc)
cc <- "pink3"
abline(v=0, col=cc)
z.txt <- bquote(paste(x == zeta(.), phantom() == .(signif(zet,6))))
mtext(z.txt, at=0, line = -1, adj = -.1, padj=.2, col=cc)
## no longer much noise (thanks to zeta.tol = 1e-5):
curve(dPareto(zet+x, alpha= 1.00001, beta= -.8), add=TRUE, col=2)
stopifnot({ rr <- range(r4$y)
2.1e-10 < rr & rr < 2.3e-10 })
showProc.time()
### ---- alpha == 1 ---------
curve(dstable(x, alpha = 1, beta = 0.3), -20, 20, log="y", n=nc)
curve(dstable(x, alpha = 1, beta = 0.3, log=TRUE), -200, 160, n=nc)
curve(dPareto(x, alpha = 1, beta = 0.3, log=TRUE), add=TRUE, col=4)
## "works", but discontinuous --- FIXME
## ditto:
curve(dstable(x, alpha=1, beta= 0.1, log=TRUE), -70,80, col=2)
curve(dPareto(x, alpha=1, beta= 0.1, log=TRUE), add=TRUE)
showProc.time()
dstable(-44, alpha=1, beta= .1)# failed
## large x gave problems at times:
dstable(-1e20, alpha = 0.9, beta = 0.8)
chkUnimodal <- function(x) {
## x = c(x1, x2) and x1 is *increasing* and x2 is *decreasing*
stopifnot((n <- length(x)) %% 2 == 0,
(n2 <- n %/% 2) >= 2)
if(is.unsorted(x[seq_len(n2)])) stop("first part is *not* increasing")
if(is.unsorted(x[n:(n2+1)])) stop("seconds part is *not* decreasing")
invisible(x)
}
showProc.time()
xLrg <- c(10^c(10:100,120, 150, 200, 300), Inf)
xLrg <- sort(c(-xLrg, xLrg))
d <- dstable(xLrg, alpha = 1.8, beta = 0.3 ); chkUnimodal(d)
d <- dstable(xLrg, alpha = 1.01, beta = 0.3 ); chkUnimodal(d) # (was slow!)
## look at the problem:
dstCurve <- function(alpha, beta, log=TRUE, NEG=FALSE,
from, to, n=nc, cLog=NULL, ...)
{
if(NEG) {
r <- curve(dstable(-x, alpha=alpha, beta=beta, log=log),
from=from, to=to, n=n, log = cLog, ...)
curve(dPareto(-x, alpha=alpha, beta=beta, log=log), add=TRUE,
col=2, lwd=2, lty=2)
} else {
r <- curve(dstable(x, alpha=alpha, beta=beta, log=log),
from=from, to=to, n=n, log = cLog, ...)
curve(dPareto(x, alpha=alpha, beta=beta, log=log), add=TRUE,
col=2, lwd=2, lty=2)
}
leg.ab <- paste0("(", if(NEG) "-x" else "x",
", a=",formatC(alpha, digits=3),
", b=",formatC(beta, digits=3),")")
legend("topright", paste0(c("dstable ", "dPareto"), leg.ab),
col=1:2, lty=1:2, lwd=1:2, bty="n")
invisible(r)
}
## (was *S.L.O.W.* on [2010-03-28] !)
r <- dstCurve(alpha = 1.01, beta = 0.3, NEG=TRUE,
from=1e10, to=1e20, cLog="x", ylim = c(-100, -45))
## zoom in:
r <- dstCurve(alpha = 1.01, beta = 0.3, , , .1e13, 9e13, ylim = c(-80, -55))
showProc.time()
d <- dstable(xLrg, alpha = 1.001, beta = -0.9) # >= 50 warnings
try( chkUnimodal(d) ) # FIXME
## look at the problem:
dstCurve(alpha = 1.001, beta = -0.9, log=FALSE, NEG=TRUE, 1e10, 1e20, cLog="xy")
## and at the right tail, too:
dstCurve(alpha = 1.001, beta = -0.9, log=FALSE, NEG=FALSE, 1000, 1e17, cLog="xy")
d <- dstable(xLrg, alpha = 1. , beta = 0.3 ); chkUnimodal(d) # "ok" now
d <- dstable(xLrg, alpha = 0.9, beta = 0.3 ) # 10 warnings (had 11)
try( chkUnimodal(d) ) # FIXME
d <- dstable(xLrg, alpha = 0.5, beta = 0.3 ) # 19 warnings (had 22)
chkUnimodal(d)
d <- dstable(xLrg, alpha = 0.1, beta = 0.3 ) # 26 warnings (had 21)
chkUnimodal(d)
showProc.time()
##------------- beta = 1 ---------------------
options(dstable.debug = TRUE)
dstable(1, alpha=1.2, beta= 1 - 1e-7)#ok
dstable(1, alpha=1.2, beta= 1)# gave error, because g(pi/2) < 0
dstable(0, alpha=13/16, beta= 1 -2^-52)# had error as g(-theta0) |-> NaN
dstable(0, alpha=19/16, beta= 1) # had error as g(pi/2) |-> NaN
options(dstable.debug = FALSE)
## NB: (beta=1, alpha = 1/2) is 'Levy' ---> dLevy() and some checks
## -- in ./pstab-ex.R
## ~~~~~~~~~~
if(doExtras) { ## actually "very-Extras" (checkLevel == "FULL")
## This needs 65 seconds (nb-mm3: 54*32*11 dstable() values)
ep <- 2^-(1:54)## beta := 1 - ep ---> 1 {NB: 1 - 2^-54 == 1 numerically}
alph.s <- (1:32)/16 # in (0, 2]
f.b1 <- sapply(alph.s, function(alf)
sapply(1 - ep, function(bet)
dstable(0:10, alpha = alf, beta = bet)),
simplify = if(getRversion() >= "2.13") "array" else TRUE)
print(summary(f.b1))
r.b1 <- range(f.b1)
stopifnot(0 < r.b1[1], r.b1[2] < 0.35)
## "FIXME" test more: monotonicity in x {mode is < 0 = min{x_i}}, beta, alpha, ...
showProc.time()
} else message("expensive dstable() checks have been skipped")
cat('Time elapsed: ', proc.time(),'\n') # "stats"
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.