DPQ (version 0.5-9)

pt_Witkovsky_Tab1: Viktor Witosky's Table_1 pt() Examples

Description

A data frame with 17 pt() examples from Witosky (2013)'s ‘Table 1’. We provide the results for the FOSS Softwares, additionally including octave's, running the original 2013 matlab code, and the corrected one from 2022.

Usage

data(pt_Witkovsky_Tab1)

Format

A data frame with 17 observations on the following numeric variables.

x

the abscissa, called q in pt().

nu

the positive degrees of freedom, called df in pt().

delta

the noncentrality parameter, called ncp in pt().

true_pnt

“true” values (computed via higher precision, see Witkovsky(2013)).

NCTCDFVW

the pt() values computed with Witkovsky's matlab implementation. Confirmed by using octave (on Fedora 40 Linux). These correspond to our R (package DPQ) pntVW13() values.

Boost

computed via the Boost C++ library; reported by Witkovsky.

R_3.3.0

computed by R version 3.3.0; confirmed to be identical using R 4.4.1

NCT2013_octave_7.3.0

values computed using Witkovsky's original matlab code, by octave 7.3.0

NCT2022_octave_8.4.0

values computed using Witkovsky's 2022 corrected matlab code, by octave 8.4.0

References

Witkovský, Viktor (2013) A Note on Computing Extreme Tail Probabilities of the Noncentral T Distribution with Large Noncentrality Parameter, Acta Universitatis Palackianae Olomucensis, Facultas Rerum Naturalium, Mathematica 52(2), 131--143.

Examples

data(pt_Witkovsky_Tab1)
stopifnot(is.data.frame(d.W <- pt_Witkovsky_Tab1), # shorter
          nrow(d.W) >= 17)
mW <- as.matrix(d.W); row.names(mW) <- NULL # more efficient
colnames(mW)[1:3] #  "x" "nu" "delta"
## use 'R pt() - compatible' names:
(n3 <- names(formals(pt)[1:3]))# "q" "df" "ncp"
colnames(mW)[1:3] <- n3
ptR <- apply(mW[, 1:3], 1, \(a3) unname(do.call(pt, as.list(a3))))
cNm <- paste0("R_", with(R.version, paste(major, minor, sep=".")))
mW <- cbind(mW, `colnames<-`(cbind(ptR), cNm),
            relErr = sfsmisc::relErrV(mW[,"true_pnt"], ptR))
mW
## is current R better than R 3.3.0?  -- or even "the same"?
all.equal(ptR, mW[,"R_3.3.0"])                    # still true in R 4.4.1
all.equal(ptR, mW[,"R_3.3.0"], tolerance = 1e-14) # (ditto)
table(ptR == mW[,"R_3.3.0"]) # {see only 4 (out of 17) *exactly* equal ??}


## How close to published NCTCDFVW is octave's run of the 2022 code?
with(d.W, all.equal(NCTCDFVW, NCT2022_octave_8.4.0, tolerance = 0)) # 3.977e-16

pVW <- apply(unname(mW[, 1:3]), 1, \(a3) unname(do.call(pntVW13, as.list(a3))))
all.equal(pVW, d.W$NCT2013_oct, tolerance = 0)# 2013-based pntVW13() --> 5.6443e-16
all.equal(pVW, d.W$NCT2022_oct, tolerance = 0)