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"