1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
|
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)
|