1require("stabledist") 2dPareto <- stabledist:::dPareto 3 4source(system.file("test-tools-1.R", package = "Matrix"), keep.source=interactive()) 5 #-> identical3(), showProc.time(),... 6(doExtras <- stabledist:::doExtras()) 7if(!require("sfsmisc")) eaxis <- axis # use sfsmisc::eaxis if available 8 9stopifnot(0 < print(dstable(4000., alpha=1.00001, beta=0.6))) 10## gave error in fBasics::dstable() 11## now 18 warnings from uniroot() 12 13pdf("dstab-ex.pdf") 14 15x <- 2^seq(0, 20, length= if(doExtras) 200 else 64) 16## Regression check for alpha=2: <==> *norm() : 17x. <- x/1024 18fx <- dstable(x., alpha = 2, beta = 0, gamma = 1.1, delta=2, pm=2) 19lf <- dstable(x., alpha = 2, beta = 0, gamma = 1.1, delta=2, pm=2, log=TRUE) 20stopifnot( 21 local({i <- is.finite(log(fx)); all.equal(log(fx[i]), lf[i])}), 22 all.equal(fx, dnorm(x., m=2, s=1.1)), 23 all.equal(lf, dnorm(x., m=2, s=1.1, log=TRUE))) 24 25fx <- dstable(x, alpha = 1.0001, beta = 0.6) 26 27plot(x,fx, log="x", type="l")# looks good 28plot(x,fx, log="xy", type="l")# --- perfect now 29stopifnot((dlx <- diff(log(fx))) < 0, # decreasing 30 { dl <- dlx[x[-1] > 1000] # the negative slope: 31 relErr(if(doExtras) -0.13934 else -0.44016, dl) < 4e-4}, 32 diff(dl) > -1e-6)# > 0 "in theory"; > -6e-7, ok on 64-bit Linux 33 34nc <- if(doExtras) 512 else 101 # number of points for curve() 35 36zeta <- function(alpha,beta) if(alpha==1) 0 else -beta*tan(pi/2*alpha) 37 38## negative beta: 39cx <- curve(dstable(x, 0.75, -.5), -.5, 1.5, n=nc)# ok, now 40m <- stableMode(0.75, -.5, tol=1e-14) 41stopifnot(all.equal(m, 0.35810298366, tol = 1e-7), 42 all.equal(dstable(m, 0.75, -.5), 0.3554664043, tol=1e-6)) 43 44showProc.time() 45 46###-------- "small" alpha ----------------- 47## alpha --> 0 --- very heavy tailed -- and numerically challenging. 48 49## symmetric (beta = 0) 50(x0 <- (-16:16)/256) 51fx0 <- dstable(x0, alpha = 0.1, beta=0, gamma = 1e6) 52plot(x0, fx0, type = "o", 53 main = expression(f(x, alpha== 0.1, beta == 0, gamma == 10^6))) 54stopifnot(all.equal(fx0[17],1.15508291498374), 55 all.equal(fx0[ 1],0.02910420736536), 56 all.equal(range(diff(fx0[1:8])), 57 c(0.0011871409, 0.0025179435), tol=1e-6) 58 ) 59 60## beta > 0 61r3 <- curve(dstable(x, alpha = 0.3, beta = 0.5, tol=1e-7), 62 -1, 1) 63m3 <- stableMode(0.3, 0.5, tol=1e-14)# still with 3 warnings 64stopifnot(all.equal(m3, -0.2505743952946, tol = 1e-10)) 65## zoom into the above 66r3. <- curve(dstable(x, alpha = 0.3, beta = 0.5, tol=1e-7), 67 -.27, -.22) 68abline(v = m3, col="gray40", lty=2) 69 70r1 <- curve(dstable(x, alpha = 0.1, beta = 0.5, tol=1e-7), 71 -.4, .2, ylim = c(0, 20), n = nc) 72m1 <- stableMode(0.1, 0.5, tol=1e-15) 73abline(v=m1, h=0, col="gray40", lty=2) 74stopifnot(all.equal(m1, -0.079193, tol=1e-5)) # -0.079193188, was -0.079192219, and -0.0791921758 75title(main = expression(f(x, alpha== 0.1, beta == 0.5))) 76## check mode *and* unimodality 77i. <- r1$x > m1 78stopifnot(## decreasing to the right: 79 diff(r1$y[ i.]) < 0, 80 ## increasing on the left: 81 diff(r1$y[!i.]) > 0) 82 83## Now *really* small alpha --> 0: 84## -------------------------- 85## Note that 86if(require("FMStable")) { 87 try( FMStable::setParam(alpha = 1e-18, location=0, logscale=0, pm = 0) ) 88## gives 89## Error in .... setParam: S0=M parametrization not suitable for small alpha 90} 91## now if this is true (and I think we can trust Geoff Robinson on this), 92## we currently have a "design bug - problem": 93## as internally, we always translate to 'pm=0' parametrization __FIXME_(how ??)__ 94## --> solution: see below: there's a simple (alpha -> 0) asymptotic 95 96## These now "work": .... well with integration warnings !! 97 98pdstab.alpha <- function(x, beta, alphas = 2^-(40:2), 99 type = "o", add=FALSE, log = "xy", ...) 100{ 101 stopifnot(is.numeric(x), length(x) == 1, length(beta) == 1, 102 is.numeric(beta), -1 <= beta, beta <= 1, 103 is.numeric(alphas), 0 <= alphas, alphas <= 2, 104 is.logical(add)) 105 if(is.unsorted(alphas)) alphas <- sort(alphas) 106 dst.alp <- vapply(alphas, function(aaa) 107 dstable(x, aaa, beta = beta, pm = 0), 1.) ## with warnings 108 if(add) 109 lines(alphas, dst.alp, type=type, ...) 110 else { 111 plot(alphas, dst.alp, type=type, log=log, axes=FALSE, 112 xlab = quote(alpha), 113 ylab = bquote(f(x == .(x), alpha)), 114 main = bquote(dstable(x == .(x), alpha, beta == .(beta), pm == 0) 115 ~~~"for (very) small" ~ alpha)) 116 if(!require("sfsmisc")) eaxis <- axis # use sfsmisc::eaxis if available 117 eaxis(1); eaxis(2) 118 } 119 invisible(cbind(alphas, dstable = dst.alp)) 120} 121 122## vary beta, keep x : 123pdstab.alpha(x = -1, beta = 0.1) 124if(doExtras) { 125pdstab.alpha(x = -1, beta = 0.3, add=TRUE, col=2, type="l") 126pdstab.alpha(x = -1, beta = 0.5, add=TRUE, col=3, type="l") 127pdstab.alpha(x = -1, beta = 0.7, add=TRUE, col=4, type="l") 128pdstab.alpha(x = -1, beta = 0.9, add=TRUE, col=5, type="l") 129 130## vary x, keep beta 131pdstab.alpha(x = -1, beta = 0.3) 132pdstab.alpha(x = +1, beta = 0.3, add=TRUE, col=2, type="l") 133pdstab.alpha(x = +5, beta = 0.3, add=TRUE, col=3, type="l") 134pdstab.alpha(x = +50, beta = 0.3, add=TRUE, col=4, type="l") 135pdstab.alpha(x = -10, beta = 0.3, add=TRUE, col=5, type="l") 136} 137## Plots suggest a simple formula 138## log(f(x, alpha, beta))= c(x,beta) + log(alpha) 139## f(x, alpha, beta) = C(x,beta) * (alpha) -- ??? 140 141## for normal alpha, it looks different {which is reassuring!} 142pdstab.alpha(x = -1, beta = 0.3, alphas = seq(1/128, 2, length=100)) 143 144## Show the formula, we derived: 145## f(x, \alpha, \beta) ~ \alpha * \frac{1 + \beta}{2e \abs{x + \pi/2 \alpha\beta}} 146## ONLY ok, when x > zeta := - beta * tan(pi/2 *alpha) 147## otherwise, "swap sign" of (x, beta) 148dst.smlA <- function(x, alpha, beta, log = FALSE) { 149 pa <- pi/2*alpha 150 i. <- x < -pa*beta 151 if(any(i.)) { 152 beta <- rep(beta, length.out=length(x)) 153 beta[i.] <- -beta[i.] 154 x [i.] <- -x [i.] 155 } 156 ## alpha*(1+beta) / (2*exp(1)*(x+ pa*beta)) 157 r <- log(alpha) + log1p(beta) - (1 + log(2*(x+ pa*beta))) 158 if(log) r else exp(r) 159} 160 161set.seed(17) 162 163alpha <- 1e-15 ## must be larger than cutoff in .fct1() in ../R/dpqr-stable.R : 164stopifnot(alpha > stabledist:::.alpha.small.dstable) 165for(beta in runif(if(doExtras) 20 else 2, -1, 1)) { 166 cat(sprintf("beta =%10g: ", beta)) 167 for(N in 1:(if(doExtras) 10 else 1)) { 168 x <- 10*rnorm(100); cat(".") 169 stopifnot(all.equal(dstable (x, alpha, beta), 170 dst.smlA(x, alpha, beta), tol = 1e-13)) 171 }; cat("\n") 172} 173 174 175curve( dstable (x, alpha=1e-10, beta=.8, log=TRUE) , -10, 10) 176curve( dst.smlA(x, alpha=1e-10, beta=.8, log=TRUE), add=TRUE, col=2) 177 178## "same" picture (self-similar !) 179curve( dstable (x, alpha=1e-10, beta=.8, log=TRUE), -1, 1) 180curve( dst.smlA(x, alpha=1e-10, beta=.8, log=TRUE), add=TRUE, col=2) 181 182 183 184## Testing stableMode here: 185 186 187### beta = 1 (extremely skewed) and small alpha: --------- 188## -------- 189## Problem at *left* ("less problematic") tail, namely very close to where the 190## support of the density becomes mathematically exactly zero : 191## 192## clear inaccuracy / bug --- maybe *seems* curable 193## 194curve(dstable(exp(x), alpha= 0.1, beta=1, pm=1, log=TRUE), -40, 10) 195## ------ 196## --> warnings both from uniroot ("-Inf") and .integrate2() 197## about equivalent to 198curve(dstable(x, alpha= 0.1, beta=1, pm=1, log=TRUE), 1e-15, 4e4, 199 log="x", xaxt="n"); eaxis(1) 200## If we decrease zeta.tol "to 0", we get better here: 201curve(dstable(exp(x), alpha= 0.1, beta=1, pm=1, log=TRUE, zeta.tol=1e-100), -40, 20) 202## or here, ... but still not good enough 203curve(dstable(exp(x), alpha= 0.1, beta=1, pm=1, log=TRUE, zeta.tol=1e-200), -45, 30) 204 205 206 207showProc.time() 208 209##------ NB: Pareto tail behavior -- see more in ./tails.R 210## ======= 211 212## alpha ~= 1 ---- and x ~ zeta(a,b) ----------- 213## ========== 214f1 <- dstable(6366.197, alpha= 1.00001, beta= .1) 215f2 <- dstable(-50929.58, alpha= 1.00001, beta= -.8) 216stopifnot(f1 > 0, f2 > 0) 217 218## these all work (luck): 219zet <- zeta(alpha= 1.00001, beta= -.8)# -50929.58 220## here, we must have larger zeta.tol: = 5e-5 is fine; now using "automatic" default 221r4 <- curve(dstable(zet+x, alpha= 1.00001, beta= -.8), -3, 7, 222 xlab = expression(zeta(alpha,beta) - x), 223 ylim=c(2.1, 2.3)*1e-10, n = nc) 224cc <- "pink3" 225abline(v=0, col=cc) 226z.txt <- bquote(paste(x == zeta(.), phantom() == .(signif(zet,6)))) 227mtext(z.txt, at=0, line = -1, adj = -.1, padj=.2, col=cc) 228## no longer much noise (thanks to zeta.tol = 1e-5): 229curve(dPareto(zet+x, alpha= 1.00001, beta= -.8), add=TRUE, col=2) 230stopifnot({ rr <- range(r4$y) 231 2.1e-10 < rr & rr < 2.3e-10 }) 232 233showProc.time() 234 235### ---- alpha == 1 --------- 236curve(dstable(x, alpha = 1, beta = 0.3), -20, 20, log="y", n=nc) 237curve(dstable(x, alpha = 1, beta = 0.3, log=TRUE), -200, 160, n=nc) 238 239curve(dPareto(x, alpha = 1, beta = 0.3, log=TRUE), add=TRUE, col=4) 240## "works", but discontinuous --- FIXME 241## ditto: 242curve(dstable(x, alpha=1, beta= 0.1, log=TRUE), -70,80, col=2) 243curve(dPareto(x, alpha=1, beta= 0.1, log=TRUE), add=TRUE) 244 245showProc.time() 246 247dstable(-44, alpha=1, beta= .1)# failed 248## large x gave problems at times: 249dstable(-1e20, alpha = 0.9, beta = 0.8) 250 251chkUnimodal <- function(x) { 252 ## x = c(x1, x2) and x1 is *increasing* and x2 is *decreasing* 253 stopifnot((n <- length(x)) %% 2 == 0, 254 (n2 <- n %/% 2) >= 2) 255 if(is.unsorted(x[seq_len(n2)])) stop("first part is *not* increasing") 256 if(is.unsorted(x[n:(n2+1)])) stop("seconds part is *not* decreasing") 257 invisible(x) 258} 259 260showProc.time() 261 262xLrg <- c(10^c(10:100,120, 150, 200, 300), Inf) 263xLrg <- sort(c(-xLrg, xLrg)) 264d <- dstable(xLrg, alpha = 1.8, beta = 0.3 ); chkUnimodal(d) 265d <- dstable(xLrg, alpha = 1.01, beta = 0.3 ); chkUnimodal(d) # (was slow!) 266## look at the problem: 267dstCurve <- function(alpha, beta, log=TRUE, NEG=FALSE, 268 from, to, n=nc, cLog=NULL, ...) 269{ 270 if(NEG) { 271 r <- curve(dstable(-x, alpha=alpha, beta=beta, log=log), 272 from=from, to=to, n=n, log = cLog, ...) 273 curve(dPareto(-x, alpha=alpha, beta=beta, log=log), add=TRUE, 274 col=2, lwd=2, lty=2) 275 } else { 276 r <- curve(dstable(x, alpha=alpha, beta=beta, log=log), 277 from=from, to=to, n=n, log = cLog, ...) 278 curve(dPareto(x, alpha=alpha, beta=beta, log=log), add=TRUE, 279 col=2, lwd=2, lty=2) 280 } 281 leg.ab <- paste0("(", if(NEG) "-x" else "x", 282 ", a=",formatC(alpha, digits=3), 283 ", b=",formatC(beta, digits=3),")") 284 legend("topright", paste0(c("dstable ", "dPareto"), leg.ab), 285 col=1:2, lty=1:2, lwd=1:2, bty="n") 286 invisible(r) 287} 288 289## (was *S.L.O.W.* on [2010-03-28] !) 290r <- dstCurve(alpha = 1.01, beta = 0.3, NEG=TRUE, 291 from=1e10, to=1e20, cLog="x", ylim = c(-100, -45)) 292## zoom in: 293r <- dstCurve(alpha = 1.01, beta = 0.3, , , .1e13, 9e13, ylim = c(-80, -55)) 294showProc.time() 295 296d <- dstable(xLrg, alpha = 1.001, beta = -0.9) # >= 50 warnings 297try( chkUnimodal(d) ) # FIXME 298## look at the problem: 299dstCurve(alpha = 1.001, beta = -0.9, log=FALSE, NEG=TRUE, 1e10, 1e20, cLog="xy") 300## and at the right tail, too: 301dstCurve(alpha = 1.001, beta = -0.9, log=FALSE, NEG=FALSE, 1000, 1e17, cLog="xy") 302 303d <- dstable(xLrg, alpha = 1. , beta = 0.3 ); chkUnimodal(d) # "ok" now 304d <- dstable(xLrg, alpha = 0.9, beta = 0.3 ) # 10 warnings (had 11) 305try( chkUnimodal(d) ) # FIXME 306d <- dstable(xLrg, alpha = 0.5, beta = 0.3 ) # 19 warnings (had 22) 307chkUnimodal(d) 308d <- dstable(xLrg, alpha = 0.1, beta = 0.3 ) # 26 warnings (had 21) 309chkUnimodal(d) 310 311showProc.time() 312 313##------------- beta = 1 --------------------- 314options(dstable.debug = TRUE) 315dstable(1, alpha=1.2, beta= 1 - 1e-7)#ok 316dstable(1, alpha=1.2, beta= 1)# gave error, because g(pi/2) < 0 317dstable(0, alpha=13/16, beta= 1 -2^-52)# had error as g(-theta0) |-> NaN 318dstable(0, alpha=19/16, beta= 1) # had error as g(pi/2) |-> NaN 319options(dstable.debug = FALSE) 320 321## NB: (beta=1, alpha = 1/2) is 'Levy' ---> dLevy() and some checks 322## -- in ./pstab-ex.R 323## ~~~~~~~~~~ 324 325if(doExtras) { ## actually "very-Extras" (checkLevel == "FULL") 326 ## This needs 65 seconds (nb-mm3: 54*32*11 dstable() values) 327 328 ep <- 2^-(1:54)## beta := 1 - ep ---> 1 {NB: 1 - 2^-54 == 1 numerically} 329 alph.s <- (1:32)/16 # in (0, 2] 330 f.b1 <- sapply(alph.s, function(alf) 331 sapply(1 - ep, function(bet) 332 dstable(0:10, alpha = alf, beta = bet)), 333 simplify = if(getRversion() >= "2.13") "array" else TRUE) 334 print(summary(f.b1)) 335 r.b1 <- range(f.b1) 336 stopifnot(0 < r.b1[1], r.b1[2] < 0.35) 337 ## "FIXME" test more: monotonicity in x {mode is < 0 = min{x_i}}, beta, alpha, ... 338 showProc.time() 339 340} else message("expensive dstable() checks have been skipped") 341 342cat('Time elapsed: ', proc.time(),'\n') # "stats" 343