R Under development (unstable) (2024-08-15 r87022 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require("stabledist")
Loading required package: stabledist
> 
> ###--- Tail approximations etc  -- both for pstable() and dstable()
> dPareto <- stabledist:::dPareto
> 
> source(system.file("test-tools-1.R", package = "Matrix"), keep.source=interactive())
Loading required package: tools
> 					#-> identical3(), showProc.time(),...
> (doExtras <- stabledist:::doExtras())
[1] FALSE
> nc <- if(doExtras) 512 else 64 # number of points for curve()
> 
> pdf("stable-tails.pdf")
> 
> pstab.tailratio <- function(alpha, beta, n = nc, prob = 1/4096,
+                             xmin = qstable(0.95, alpha,beta, tol = 0.01),
+                             xmax = qstable(1 - prob, alpha,beta))
+ {
+   ## Purpose: Explore eps(x) where  (1 - pstable(x))/(1- pPareto()) = 1 + eps(x)
+   ##  <==>
+   ## ----------------------------------------------------------------------
+   ## Arguments:
+   ## ----------------------------------------------------------------------
+   ## Author: Martin Maechler, Date: 21 Mar 2011, 10:09
+     cl <- match.call()
+     stopifnot(0 < prob, prob < .5,
+               0 < xmin, xmin < xmax)
+     x <- exp(seq(log(xmin), log(xmax), length.out = n))
+     iF <-               pstable(x, alpha,beta, lower.tail=FALSE)
+     ok <- iF > 0
+     iF <- iF[ok]
+     x <- x[ok]
+     iFp <- stabledist:::pPareto(x, alpha,beta, lower.tail=FALSE)
+     eps <- (iF - iFp)/iFp
+     structure(list(x=x, eps=eps, call = cl, alpha=alpha, beta=beta),
+               class = "pstableTailratio")
+ }
> 
> plot.pstableTailratio <- function(x, type="l", col="blue3",
+                                   lin.col = adjustcolor("red2",.6), ...)
+ {
+   ## Purpose:
+   ## ----------------------------------------------------------------------
+   ## Arguments:
+   ## ----------------------------------------------------------------------
+     stopifnot(is.list(x), is.numeric(x$eps))
+     tit <- substitute("tail ratio approx. for"~~ pstable(alpha==A, beta==B),
+                       list(A=x[["alpha"]], B=x[["beta"]]))
+     dat <- as.data.frame(x[c("x","eps")])
+     ## NOTA BENE: Empirically,  I see that  eps > 0 <==> alpha > 1
+     ##                                      eps < 0 <==> alpha < 1
+     dat <- dat[dat[,"eps"] > 0, ] ## drop invalid  eps[]
+     plot(eps ~ x, log = "xy", data = dat, col=col,
+          ylab = expression(epsilon ~~~ "('eps')"),
+          main = tit, type=type, ...)
+     mtext( expression(epsilon(x) == (bar(F)(x,.) - bar(F)[P](x,.)) / bar(F)[P](x,.)) )
+     fm <- lm(log(eps) ~ log(x), weights = x^2, data = dat)
+     lines(dat[["x"]], exp(predict(fm)), col=lin.col)
+     Form <- function(x) formatC(x, digits=4, wid=1)
+     leg.line <-
+         substitute(log(epsilon) == A + B * log(x),
+                    list(A = Form(coef(fm)[["(Intercept)"]]),
+                         B = Form(coef(fm)[["log(x)"]])))
+     legend("topright", legend=leg.line, bty = "n", lty=1, col=lin.col)
+ }
> 
> 
> plot(tr0  <- pstab.tailratio(1, 0.5))
> plot(tr1  <- pstab.tailratio(1.1, 0.25))
> plot(tr2 <- pstab.tailratio(0.99, +0.992))
> 
> showProc.time()
Time (user system elapsed): 0.34 0 0.34 
> 
> plot(tr   <- pstab.tailratio(1.2, 0.5))
> 
> plot(tr3 <- pstab.tailratio(0.7, +0.9))
> 
> plot(tr4 <- pstab.tailratio(1.7, +0.6))# not really useful: pstable(.) = 1 too early
> 
> showProc.time()
Time (user system elapsed): 0.17 0 0.17 
> 
> ##---------------- Now the density
> 
> ##' @title Explore eps(x) where   dstable(x)/dPareto(x) = 1 + eps(x)
> ##' @param alpha
> ##' @param beta
> ##' @param n
> ##' @param prob
> ##' @param xmin
> ##' @param xmax
> ##' @return an object of \code{\link{class} "dstableTailratio"}, ...
> ##' @author Martin Maechler, 21 Mar 2011
> dstab.tailratio <- function(alpha, beta, n = nc, prob = 1/4096,
+                             xmin = qstable(0.95, alpha,beta, tol = 0.01),
+                             xmax = qstable(1 - prob, alpha,beta))
+ {
+     cl <- match.call()
+     stopifnot(0 < prob, prob < .5,
+               0 < xmin, xmin < xmax)
+     x <- exp(seq(log(xmin), log(xmax), length.out = n))
+     f <-               dstable(x, alpha,beta)
+     ok <- f > 0
+     f <- f[ok]
+     x <- x[ok]
+     fp <- stabledist:::dPareto(x, alpha,beta)
+     eps <- (f - fp)/fp
+     structure(list(x=x, eps=eps, call = cl, alpha=alpha, beta=beta),
+               class = "dstableTailratio")
+ }
> 
> ##' @title plot() method for  dstableTailratio() results
> ##' @param x object of \code{\link{class} "dstableTailratio"}.
> ##' @param type plot type; default simple lines
> ##' @param col
> ##' @param lin.col
> ##' @param ... optional further arguments passed to \code{\link{plot.formula}()}.
> ##' @return
> ##' @author Martin Maechler
> plot.dstableTailratio <- function(x, type="l", col="blue3",
+                                   lin.col = adjustcolor("red2",.6), ...)
+ {
+   ## Purpose:
+   ## ----------------------------------------------------------------------
+   ## Arguments:
+   ## ----------------------------------------------------------------------
+     stopifnot(is.list(x), is.numeric(x$eps))
+     tit <- substitute("tail ratio approx. for"~~ dstable(alpha==A, beta==B),
+                       list(A=x[["alpha"]], B=x[["beta"]]))
+     dat <- as.data.frame(x[c("x","eps")])
+     ## NOTA BENE: Empirically,  I see that  eps > 0 <==> alpha > 1
+     ##                                      eps < 0 <==> alpha < 1
+     dat <- dat[dat[,"eps"] > 0, ] ## drop invalid  eps[]
+     plot(eps ~ x, log = "xy", data = dat, col=col,
+          ylab = expression(epsilon ~~~ "('eps')"),
+          main = tit, type=type, ...)
+     mtext( expression(epsilon(x) == (f(x,.) - f[P](x,.)) / f[P](x,.)) )
+     fm <- lm(log(eps) ~ log(x), weights = x^2, data = dat)
+     lines(dat[["x"]], exp(predict(fm)), col=lin.col)
+     Form <- function(x) formatC(x, digits=4, wid=1)
+     leg.line <-
+         substitute(log(epsilon) == A + B * log(x),
+                    list(A = Form(coef(fm)[["(Intercept)"]]),
+                         B = Form(coef(fm)[["log(x)"]])))
+     legend("topright", legend=leg.line, bty = "n", lty=1, col=lin.col)
+ }
> 
> plot(fr   <- dstab.tailratio(1.01, 0.8))
> plot(fr   <- dstab.tailratio(1.05, 0.4))
> plot(fr   <- dstab.tailratio(1.1,  0.4))
> plot(fr   <- dstab.tailratio(1.2,  0.5))
> plot(fr   <- dstab.tailratio(1.3,  0.6))
> 
> showProc.time()
Time (user system elapsed): 0.53 0 0.53 
> 
> plot(fr   <- dstab.tailratio(1.4, 0.7))
> plot(fr   <- dstab.tailratio(1.5, 0.8))
> 
> plot(fr   <- dstab.tailratio(1.5, 0.8, xmax= 1000))
> plot(fr   <- dstab.tailratio(1.5, 0.8, xmax= 1e4));abline(v=1000, lty=2)
> plot(fr   <- dstab.tailratio(1.5, 0.8, xmax= 1e5));abline(v=1e4, lty=2)
> 
> showProc.time()
Time (user system elapsed): 0.41 0 0.41 
> 
> plot(fr   <- dstab.tailratio(1.6, 0.9))
> plot(fr   <- dstab.tailratio(1.7, 0.1))
> plot(fr   <- dstab.tailratio(1.8, 0.2))
> 
> showProc.time()
Time (user system elapsed): 0.31 0.02 0.33 
> 
> ##------ Some explicit tail problems visualized:
> 
> I <- integrate(dstable, 0,Inf, alpha=0.998, beta=0, subdivisions=1000)
There were 24 warnings (use warnings() to see them)
> str(I)
List of 5
 $ value       : num 0.5
 $ abs.error   : num 5.14e-05
 $ subdivisions: int 1
 $ message     : chr "OK"
 $ call        : language integrate(f = dstable, lower = 0, upper = Inf, alpha = 0.998, beta = 0,      subdivisions = 1000)
 - attr(*, "class")= chr "integrate"
> stopifnot(abs(I$value - 0.5) < 1e-4)
> 
> curve(dstable(x, alpha=.998, beta=0,    log=TRUE), 10, 1e17, log="x")
There were 22 warnings (use warnings() to see them)
> curve(dstable(x, alpha=.999, beta=0.1,  log=TRUE), 10, 1e17, log="x")
There were 20 warnings (use warnings() to see them)
> curve(dstable(x, alpha=.999, beta=0.9,  log=TRUE), 10, 1e17, log="x")
There were 22 warnings (use warnings() to see them)
> curve(dstable(x, alpha=.999, beta=0.99, log=TRUE), 10, 1e17, log="x")
There were 22 warnings (use warnings() to see them)
> curve(dstable(x, alpha=.999, beta=0.99, log=TRUE), 10, 1e170, log="x")
Warning message:
In .integrate2(g1, lower = a, upper = b, subdivisions = subdivisions,  :
  roundoff error is detected in the extrapolation table
> showProc.time()
Time (user system elapsed): 0.52 0.01 0.53 
> 
> ## less problems when alpha > ~= 1 (but it's  __S..L..O..W__ !)
> curve(dstable(x, alpha=1.001,beta=0.99, log=TRUE), 10,  1e7, log="x")
There were 50 or more warnings (use warnings() to see the first 50)
> curve(dstable(x, alpha=1.001,beta=0.99, log=TRUE), 10, 1e17, log="x")
There were 50 or more warnings (use warnings() to see the first 50)
> ## -> problem --> zoom in:
> curve(dstable(x, alpha=1.001,beta=0.99, log=TRUE), 1e12, 160e12)
There were 50 or more warnings (use warnings() to see the first 50)
> curve(dPareto(x, alpha=1.001,beta=0.99, log=TRUE), add=TRUE, lty=3, col=4)
> 
> curve(dstable(x, alpha=1.001,beta=0.99, log=TRUE), 10, 1e40, log="x")
There were 50 or more warnings (use warnings() to see the first 50)
> showProc.time()
Time (user system elapsed): 1.22 0.03 1.25 
> 
> ## NB: alpha == 1   also has problems in tail  --- only as long as "old R"s wrong uniroot is used:
> curve(dstable(x, alpha=1.   ,beta=0.99, log=TRUE), 1, 20)
> curve(dstable(x, alpha=1.   ,beta=0.99, log=TRUE), 1,100)
> 
> 
> showProc.time()
Time (user system elapsed): 0.25 0 0.25 
> 
> proc.time()
   user  system elapsed 
   3.89    0.10    3.98