qbBaha2017 {DPQmpfr}

Accurate qbeta() values from Baharev et al (2017)'s Program

Description

Compuate "accurate" qbeta() values from Baharev et al (2017)'s Program.

Usage

data("qbBaha2017")

Format

FIXME: Their published table only shows 6 digits, but running their (32-bit statically linked) Linux executable ‘mindiffver’ (from their github repos, see "source") with their own ‘input.txt’ gives 12 digits accuracy, which we should be able to increase even more, see https://github.com/baharev/mindiffver/blob/master/README.md

A numeric matrix, 9 \times 22 with guaranteed accuracy qbeta(0.95, a,b) values, for a = 0.5, 1, 1.5, 2, 2.5, 3, 5, 10, 25 and b = with str()

   num [1:9, 1:22] 0.902 0.95 0.966 0.975 0.98 ...
   - attr(*, "dimnames")=List of 2
    ..$ a: chr [1:9] "0.5" "1" "1.5" "2" ...
    ..$ b: chr [1:22] "1" "2" "3" "4" ...

Details

MM constructed this data as follows (TODO: say more..):

        ff <- "~/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/Baharev_et_al-2017_table3.txt"
        qbB2017 <- t( data.matrix(read.table(ff)) )
        dimnames(qbB2017) <- dimnames(qbet)
        saveRDS(qbB2017, "..../qbBaha2017.rds")
    

Source

This matrix comprises all entries of Table 3, p. 776 of
Baharev, A., Schichl, H. and Rév, E. (2017) Computing the noncentral-F distribution and the power of the F-test with guaranteed accuracy; Comput. Stat. 32(2), 763–779. doi:10.1007/s00180-016-0701-3

The paper mentions the first author's ‘github’ repos where source code and executables are available from: https://github.com/baharev/mindiffver/

Examples

data(qbBaha2017)
str(qbBaha2017)
str(ab <- lapply(dimnames(qbBaha2017), as.numeric))
stopifnot(ab$a == c((1:6)/2, 5, 10, 25),
          ab$b == c(1:15, 10*c(2:5, 10, 25, 50)))
matplot(ab$b, t(qbBaha2017)[,9:1], type="l", log = "x", xlab = "b",
        ylab = "qbeta(.95, a,b)",
        main = "Guaranteed accuracy 95% percentiles of Beta distribution")
legend("right", paste("a = ", format(ab$a)),
       lty=1:5, col=1:6, bty="n")

## Relative error of R's qbeta() -- given that the table only shows 6
## digits, there is *no* relevant error: R's qbeta() is accurate enough:
x.ab <- do.call(expand.grid, ab)
matplot(ab$b, 1 - t(qbeta(0.95, x.ab$a, x.ab$b) / qbBaha2017),
        main = "rel.error of R's qbeta() -- w/ 6 digits, it is negligible",
        ylab = "1 - qbeta()/'true'",
        type = "l", log="x", xlab="b")
abline(h=0, col=adjustcolor("gray", 1/2))