tests/special-fun-dgamma.R

In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

### dgamma(): ----------------------- was part of  ./special-fun-ex.R -------------------
stopifnot(require("Rmpfr"))
require(sfsmisc)# -> eaxis(); relErrV()

(doExtras <- Rmpfr:::doExtras())
options(nwarnings = 50000, width = 99)

##                                      vvvvvvvvvvvvvvvv
## to enhance  |rel.Err| plots:  from ./special-fun-ex.R, also in ~/R/Pkgs/DPQ/tests/pow-tst.R }
drawEps.h <- function(p2 = -(53:51), side = 4, lty=3, lwd=2, col=adjustcolor(2, 1/2)) {
    abline(h = 2^p2, lty=lty, lwd=lwd, col=col)
    axis(side, las=2, line=-1, at = 2^p2,
         labels = as.expression(lapply(p2, function(p) substitute(2^E, list(E=p)))),
         col.axis = col, col=NA, col.ticks=NA)
}

(do.pdf <- !dev.interactive(orNone = TRUE))

if(do.pdf) pdf("special-fun-dgamma.pdf")

xe <- c(-2e5, -1e5, -2e4, -1e4, -2000, -1000, -500, -200, -100, -50, -20, -10)
(xe <- c(xe, -8:8, -rev(xe)))
two <- mpfr(2, 256)
## For centering at E[.], will use xP(x, shp) :
xP <- function(x, d) {
    ## cannot eliminate them, as for <mpfr> they are all finite ..
    ## x <- x[is.finite(x)]
    x - d*(x > d)
}
aEQformat <- function(xy, ...) format(xy, digits = 7, ...)
allEQ_0 <- function (target, current, ...)
    all.equal(target, current, tolerance = 0, formatFUN = aEQformat, ...)
stopIfNot <-
    if("allow.logical0" %in% names(formals(stopifnot))) { # experimental (MM only)
        stopifnot
    } else function(exprs, allow.logical0) stopifnot(exprs=exprs)
abs19 <- function(r) pmax(abs(r), 1e-19) # cut  |err| to positive {for log-plots}

for(shp in 2^c(-20, -3, -1:1, 4, 10, 14, 20, 50)) {
    cat("shape = 2^", log2(shp), ":\n-------------\n")
    d.dg  <- dgamma(xP(2 ^ xe, shp) -> x, shape=shp)
    m.dg  <- dgamma(xP(two^xe, shp), shape=shp)
    m.ldg <- dgamma(xP(two^xe, shp), shape=shp, log=TRUE)
    relE <- asNumeric(relErrV(m.dg, d.dg))
    ## Plots: do *not* observe any problems yet
    plot(x, relE, log="x", type="l",
         main = paste0("rel.Errors dgamma(., shape = 2^", log2(shp),")"))
    abline(h=0, col=adjustcolor("gray10", 1/2), lty=3, lwd=2)
    plot(x, abs19(relE), log="xy", type="l", ylim = pmax(4e-17, range(abs19(relE), finite=TRUE)))
    abline(h = 2^-(52:50), col=adjustcolor("red4",1/2), lty=3)
    ##
    stopIfNot(exprs = {
        !is.unsorted(xe)
        is.finite(m.dg)
        m.dg >= 0
        shp > 1  || all(diff(m.dg) <= 0)
        shp > 100|| all((m.dg > 0) >= (d.dg > 0))
        any(fin.d <- is.finite(d.dg))
        m.dg[!fin.d] > 1e300
        { cat("all.EQ(<mpfr>, <doubl>):", allEQ_0  (m.dg[fin.d], d.dg[fin.d]), "\n")
          shp > 100  ||                   all.equal(m.dg[fin.d], d.dg[fin.d],
                                                    tol = 1e-13) # 2.063241e-14
        }
        ## compare with log scale :
        if(any(pos.d <- m.dg > 0)) {
            cat("For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)):",
                allEQ_0  (log(m.dg[pos.d]), m.ldg[pos.d]),"\n")
            ##
            all.equal(log(m.dg[pos.d]), m.ldg[pos.d], tol = 1e-14)
        } else TRUE
    })#, allow.logical0 = TRUE)
}

## NB:  dgamma(x, sh)  sh >= 1 calls
## --   dpois_raw(sh-1, x)   which then
## adds stirlerr(x.) to bd0(x., lambda) ;  where x. <- sh-1; lambda <- x
## bd0(x,L) ~= 0 iff x ~= L  <==> sh-1 ~= x  <==>  x+1 ~= sh

sh2x_gamma <- function(sh, nx, k = 12, f1 = 0.5, f2 = 1.25) {
    stopifnot(is.numeric(sh), length(sh) == 1, sh >= 0, length(k) == 1, k >= 3,
              f1 < f2, length(f1) == length(f2), length(f2) == 1)
    p2 <- 2^-(k:3)
    1 + sh* unique(sort(c(1-p2, 1, 1+p2, # <- values x very close to sh -- does *not* make any diff (????)
                          seq(f1, f2, length=nx))))
}

relEgamma <- function(sh, nx = 1001, k = 12, precBits = 256,
                      x = sh2x_gamma(sh, nx=nx, k=k)) {
    dg  <- dgamma(x, sh)
    dgM <- dgamma(mpfr(x,  precBits),
                  mpfr(sh, precBits))
    structure(cbind(x, relE = asNumeric(relErrV(dgM, dg))), shape=sh)
}

shs <- 1/32 + seq(6, 16, by = 1/8)
stopifnot(all(!is.whole(shs*2))) # shs are *not* half-integers
system.time(
    LrelE <- lapply(shs, relEgamma)
) # 5.3 sec

m.relE <- sapply(LrelE, function(m) m[,"relE"])
qrelE <- t(apply(abs(m.relE), 2, quantile, probs = c(.05, .25, .50, .75, .90, 1)))
##               ^^^^^^^^^^^

## Heureka! --- this shows quite a difference between R 4.3.3  and R-devel (R 4.4.0) !!

iS <- sort.list(qrelE[,"50%"], decreasing = TRUE)
cbind(shs, qrelE)[iS,]
## For R 4.3.3 :
##      shs             5%          25%          50%          75%          90%         100%
## 14.53125   9.410630e-15 9.815160e-15 1.023178e-14 1.065722e-14 1.092232e-14 1.138372e-14
## 15.03125   8.265317e-15 8.702900e-15 9.086072e-15 9.506915e-15 9.756106e-15 1.007928e-14
## 15.90625   6.799207e-15 7.137733e-15 7.611360e-15 8.057580e-15 8.343992e-15 8.670817e-15
## 13.53125   6.716182e-15 7.103502e-15 7.566360e-15 8.004966e-15 8.276645e-15 8.630780e-15
## 15.65625   6.031124e-15 6.389848e-15 6.803347e-15 7.261310e-15 7.527491e-15 8.031559e-15
## ..........
## ..........

myRversion <- paste(R.version.string, "--", osVersion)
if((mach <- Sys.info()[["machine"]]) != "x86_64")
    myRversion <- paste0(myRversion, "_", mach)
if(!capabilities("long.double"))
    myRversion <- paste0(myRversion, "_no_LDbl")
myRversion

rngP <- function(y, M = 1e-14) { yr <- range(y); if(yr[2] < M) yr[2] <- M; yr }

boxplot(abs19(m.relE), at = shs, log="y", ylim = c(7e-17, rngP(abs(m.relE))[2]), yaxt="n")
eaxis(2); drawEps.h(); mtext(myRversion, adj=1, cex=3/4)

matplot(shs, qrelE, type="l", log="y", yaxt="n", ylim = rngP(qrelE))
title("|relErr( dgamma(x, sh) |   for  x / sh  in [.5, 1.25]")
eaxis(2); drawEps.h(); mtext(myRversion, adj=1, cex=3/4)




cat('Time elapsed: ', proc.time(),'\n') # "stats"
if(!interactive()) warnings()