qbBaha2017 | R Documentation |
Compuate "accurate" qbeta()
values from
Baharev et al (2017)'s Program.
data("qbBaha2017")
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" ...
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")
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.
\Sexpr[results=rd]{tools:::Rd_expr_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/
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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.