1require("stabledist") 2pPareto <- stabledist:::pPareto 3 4source(system.file("test-tools-1.R", package = "Matrix"), keep.source=interactive()) 5 #-> identical3(), showProc.time(),... 6(doExtras <- stabledist:::doExtras()) 7 8options(pstable.debug = FALSE) 9options(pstable.debug = TRUE)# want to see when uniroot() is called 10 11stopifnot(all.equal(pstable(0.3, 0.75, -.5, tol= 1e-14), 12 0.6688726496, tol = 1e-8)) 13## was 0.66887227658457, tol = 1e-10)) 14pstable(-4.5, alpha = 1, beta = 0.01)## gave integration error (now uniroot..) 15 16## a "outer vectorized" version: 17pstabALL <- function(x, alpha, beta, ...) 18 sapply(alpha, function(alph) 19 sapply(beta, function(bet) 20 pstable(x, alph, bet, ...))) 21 22alph.s <- (1:32)/16 # in (0, 2] 23beta.s <- (-16:16)/16 # in [-1, 1] 24stopifnot(pstabALL( Inf, alph.s, beta.s) == 1, 25 pstabALL(-Inf, alph.s, beta.s, log.p=TRUE) == -Inf, 26 pstabALL( 0, alph.s, beta = 0) == 0.5, 27 TRUE) 28 29pdf("pstab-ex.pdf") 30 31##---- log-scale ------------- 32r <- curve(pstable(x, alpha=1.8, beta=.9, 33 lower.tail=FALSE, log.p=TRUE), 34 5, 150, n=500, 35 log="x",type="b", cex=.5) 36curve(pPareto(x, alpha=1.8, beta=.9, 37 lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2) 38##--> clearly potential for improvement! 39 40## the less extreme part - of that: 41r <- curve(pstable(x, alpha=1.8, beta=.9, 42 lower.tail=FALSE, log.p=TRUE), 43 1, 50, n=500, log="x") 44curve(pPareto(x, alpha=1.8, beta=.9, lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2) 45 46## Check that pstable() is the integral of dstable() --- using simple Simpson's rule 47 48## in it's composite form: 49## \int_a^b f(x) dx\approx \frac{h}{3} 50## \bigg[ f(x_0) + 2 \sum_{j=1}^{n/2-1}f(x_{2j}) + 51## + 4 \sum_{j=1}^{n/2} f(x_{2j-1}) + 52## + f(x_n) \bigg], 53intSimps <- function(fx, h) { 54 stopifnot((n <- length(fx)) %% 2 == 0, 55 n >= 4, length(h) == 1, h > 0) 56 n2 <- n %/% 2 57 j2 <- 2L * seq_len(n2-1) 58 j4 <- 2L * seq_len(n2) - 1L 59 h/3 * sum(fx[1], 2* fx[j2], 4* fx[j4], fx[n]) 60} 61 62chk.pd.stable <- function(alpha, beta, xmin=NA, xmax=NA, 63 n = 256, do.plot=TRUE, 64 comp.tol = 1e-13, eq.tol = 1e-3) 65{ 66 stopifnot(n >= 20) 67 if(is.na(xmin)) xmin <- qstable(0.01, alpha, beta) 68 if(is.na(xmax)) xmax <- qstable(0.99, alpha, beta) 69 dx <- ceiling(1024*grDevices::extendrange(r = c(xmin, xmax), f = 0.01))/1024 70 h <- diff(dx)/n 71 x <- seq(dx[1], dx[2], by = h) 72 fx <- dstable(x, alpha=alpha, beta=beta, tol= comp.tol) 73 Fx <- pstable(x, alpha=alpha, beta=beta, tol=2*comp.tol) 74 i.ev <- (i <- seq_along(x))[i %% 2 == 0 & i >= max(n/10, 16)] 75 ## integrate from x[1] up to x[i] (where i is even); 76 ## the exact value will be F(x[i]) - F(x[1]) == Fx[i] - Fx[1] 77 Fx. <- vapply(lapply(i.ev, seq_len), 78 function(ii) intSimps(fx[ii], h), 0) 79 a.eq <- all.equal(Fx., Fx[i.ev] - Fx[1], tol = eq.tol) 80 if(do.plot) { 81 ## Show the fit 82 plot(x, Fx - Fx[1], type = "l") 83 lines(x[i.ev], Fx., col=adjustcolor("red", 0.5), lwd=3) 84 op <- par(ask=TRUE) ; on.exit(par(op)) 85 ## show the "residual", i.e., the relative error 86 plot(x[i.ev], 1- Fx./(Fx[i.ev] - Fx[1]), 87 type = "l", xlim = range(x)) 88 abline(h=0, lty=3, lwd = .6) 89 } 90 91 if(!isTRUE(a.eq)) stop(a.eq) 92 invisible(list(x=x, f=fx, F=Fx, i. = i.ev, F.appr. = Fx.)) 93} 94 95op <- par(mfrow=2:1, mar = .1+c(3,3,1,1), mgp=c(1.5, 0.6,0)) 96 97c1 <- chk.pd.stable(.75, -.5, -1, 1.5, eq.tol = .006) 98c2 <- chk.pd.stable(.95, +0.6, -1, 1.5, eq.tol = .006)# with >= 50 warnings 99## here are the "values" 100with(c1, all.equal(F.appr., F[i.] - F[1], tol = 0)) # (.0041290 on 64-Lnx) 101with(c2, all.equal(F.appr., F[i.] - F[1], tol = 0)) # (.0049307 on 64-Lnx) 102 103showProc.time() # 104 105c3 <- chk.pd.stable(.95, +0.9, -3, 15) # >= 50 warnings 106 107curve(dstable(x, .999, -0.9), -20, 5, log="y") 108curve(pstable(x, .999, -0.9), -20, 5, log="y")#-> using uniroot 109c4 <- chk.pd.stable(.999, -0.9, -20, 5) 110 111showProc.time() # 112 113## alpha == 1 , small beta ---- now perfect 114curve(pstable(x, alpha=1, beta= .01), -6, 8, ylim=0:1) 115abline(h=0:1, v=0, lty=3, col="gray30") 116n <- length(x <- seq(-6,8, by = 1/16)) 117px <- pstable(x, alpha=1, beta= .01) 118## now take approximation derivative by difference ratio: 119x. <- (x[-n]+x[-1])/2 120plot (x., diff(px)*16, type="l") 121## now check convexity/concavity : 122f2 <- diff(diff(px)) 123stopifnot(f2[x[-c(1,n)] < 0] > 0, 124 f2[x[-c(1,n)] > 0] < 0) 125## and compare with dstable() ... which actually shows dstable() problem: 126fx. <- dstable(x., alpha=1, beta=.01) 127lines(x., fx., col = 2, lwd=3, lty="5111") 128 129if(dev.interactive(orNone=TRUE)) { 130 curve(dstable(x, 1., 0.99), -6, 50, log="y")# "uneven" (x < 0); 50 warnings 131 curve(dstable(x, 1.001, 0.95), -10, 30, log="y")# much better 132} 133showProc.time() # 134 135if(doExtras) { 136c5 <- chk.pd.stable(1., 0.99, -6, 50)# -> uniroot 137c6 <- chk.pd.stable(1.001, 0.95, -10, 30)# -> uniroot; 2nd plot *clearly* shows problem 138with(c5, all.equal(F.appr., F[i.] - F[1], tol = 0)) # .00058 on 64-Lnx 139with(c6, all.equal(F.appr., F[i.] - F[1], tol = 0)) # 2.611e-5 on 64-Lnx 140 141## right tail: 142try(## FIXME: 143c1.0 <- chk.pd.stable(1., 0.8, -6, 500)# uniroot; rel.difference = .030 144) 145## show it more clearly 146curve(pstable(x, alpha=1, beta=0.5), 20, 800, log="x", ylim=c(.97, 1)) 147curve(pPareto(x, alpha=1, beta=0.5), add=TRUE, col=2, lty=2) 148abline(h=1, lty=3,col="gray") 149# and similarly (alpha ~= 1, instead of == 1): 150curve(pstable(x, alpha=1.001, beta=0.5), 20, 800, log="x", ylim=c(.97, 1)) 151curve(pPareto(x, alpha=1.001, beta=0.5), add=TRUE, col=2, lty=2) 152abline(h=1, lty=3,col="gray") 153## zoom 154curve(pstable(x, alpha=1.001, beta=0.5), 100, 200, log="x") 155curve(pPareto(x, alpha=1.001, beta=0.5), add=TRUE, col=2, lty=2) 156## but alpha = 1 is only somewhat better as approximation: 157curve(pstable(x, alpha=1 , beta=0.5), add=TRUE, col=3, 158 lwd=3, lty="5131") 159showProc.time() # 160} 161 162 163c7 <- chk.pd.stable(1.2, -0.2, -40, 30) 164c8 <- chk.pd.stable(1.5, -0.999, -40, 30)# two warnings 165 166showProc.time() # 167 168### Newly found -- Marius Hofert, 18.Sept. (via qstable): 169stopifnot(all.equal(qstable(0.6, alpha = 0.5, beta = 1, 170 tol=1e-15, integ.tol=1e-15), 171 2.636426573120147)) 172##--> which can be traced to the first of 173stopifnot(pstable(q= -1.1, alpha=0.5, beta=1) == 0, 174 pstable(q= -2.1, alpha=0.6, beta=1) == 0) 175 176## Found by Tobias Hudec, 2 May 2015: 177stopifnot( 178 all.equal(1.5, qstable(p=0.5, alpha=1.5, beta=0, gamma=2, delta = 1.5)), 179 all.equal(1.5, qstable(p=0.5, alpha=0.6, beta=0, gamma=0.2, delta = 1.5)) 180) 181 182## Stable(alpha = 1/2, beta = 1, gamma, delta, pm = 1) <===> Levy(delta, gamma) 183source(system.file("xtraR", "Levy.R", package = "stabledist"), keep.source=interactive()) 184##-> dLevy(x, mu, c, log) and 185##-> pLevy(x, mu, c, log.p, lower.tail) 186 187set.seed(101) 188show.Acc <- (interactive() && require("Rmpfr")) 189if(show.Acc) { ## want to see accuracies, do not stop "quickly" 190 format.relErr <- function(tt, cc) 191 format(as.numeric(relErr(tt, cc)), digits = 4, width = 8) 192} 193## FIXME: Look why pstable() is so much less accurate than dstable() 194## even though the integration in dstable() is more delicate in general 195 196pTOL <- 1e-6 # typically see relErr of 5e-7 197dTOL <- 1e-14 # typically see relErr of 1.3...3.9 e-15 198showProc.time() 199 200## Note that dstable() is more costly than pstable() 201for(ii in 1:(if(doExtras) 32 else 8)) { 202 Z <- rnorm(2) 203 mu <- sign(Z[1])*exp(Z[1]) 204 sc <- exp(Z[2]) 205 x <- seq(mu, mu+ sc* 100*rchisq(1, df=3), 206 length.out= if(doExtras) 512 else 32) 207 ## dLevy() and pLevy() using only pnorm() are "mpfr-aware": 208 x. <- if(show.Acc) mpfr(x, 256) else x 209 pL <- pLevy(x., mu, sc) 210 pS <- pstable(x, alpha=1/2, beta=1, gamma=sc, delta=mu, 211 pm = 1) 212 dL <- dLevy(x., mu, sc) 213 dS <- dstable(x, alpha=1/2, beta=1, gamma=sc, delta=mu, 214 pm = 1) 215 if(show.Acc) { 216 cat("p: ", format.relErr(pL, pS), "\t") 217 cat("d: ", format.relErr(dL, dS), "\n") 218 } else { 219 cat(ii %% 10) 220 } 221 stopifnot(all.equal(pL, pS, tol = pTOL), 222 all.equal(dL, dS, tol = dTOL)) 223}; cat("\n") 224showProc.time()## ~ 75 sec (doExtras=TRUE) on lynne (2012-09) 225