pRangeStat <-
function(q, sigma=1, nPerGroup=5, nGroups=1, stat=c("ES", "FoM", "D"),
lower.tail=TRUE, loUp) {
stat <- match.arg(toupper(stat), choices=c("ES", "FOM", "D"))
stat <- c(ES="ES", FOM="FoM", D="D")[stat]
nPerGroup <- as.integer(nPerGroup[1])
nGroups <- as.integer(nGroups[1])
stopifnot(nPerGroup > 1L, nPerGroup <= max(shotGroups::DFdistr[["n"]]),
nGroups > 0L, nGroups <= max(shotGroups::DFdistr[["nGroups"]]))
sigma <- sigma[1]
stopifnot(is.numeric(sigma), sigma > 0)
## available probabilities -> range is search interval for uniroot()
allPV <- regmatches(names(shotGroups::DFdistr),
regexpr(paste0("^", stat, "_Q[[:digit:]]{3}$"),
names(shotGroups::DFdistr)))
if(missing(loUp)) {
## use whole probability range from lookup table for search
loUp <- range(as.numeric(sub(paste0(stat, "_Q([[:digit:]]{3})$"), "\\1" , allPV)) / 1000)
}
# pp <- numeric(length(q)) # initialize probabilities to 0
pp <- rep(NA_real_, length(q)) # initialize probabilities to NA
keep <- which((q >= 0) | !is.finite(q)) # keep non-negative q, NA, NaN, -Inf, Inf
qdf <- function(p, q, sigma, nPerGroup, nGroups, stat, lower.tail) {
qobs <- qRangeStat(p=p, sigma=sigma, nPerGroup=nPerGroup, nGroups=nGroups,
stat=stat, lower.tail=lower.tail)
# print(c(p, q, qobs, qobs-q))
qobs - q
}
getP <- function(q, sigma, nPerGroup, nGroups, stat, loUp, lower.tail) {
tryCatch(uniroot(qdf, interval=loUp, q=q, sigma=sigma,
nPerGroup=nPerGroup, nGroups=nGroups,
stat=stat, lower.tail=lower.tail)$root,
error=function(e) { return(NA_real_) })
}
pp[keep] <- unlist(Map(getP, q=q[keep], sigma=list(sigma),
nPerGroup=nPerGroup, nGroups=nGroups,
stat=stat, loUp=list(loUp),
lower.tail=lower.tail[1]))
## special cases that are distribution independent
if(lower.tail) {
pp[q == -Inf] <- 0
pp[q == Inf] <- 1
} else {
pp[q < 0] <- 1
pp[q == Inf] <- 0
}
return(pp)
}
qRangeStat <-
function(p, sigma=1, nPerGroup=5, nGroups=1, stat=c("ES", "FoM", "D"),
method=c("linear", "spline"),
lower.tail=TRUE) {
stat <- match.arg(toupper(stat), choices=c("ES", "FOM", "D"))
stat <- c(ES="ES", FOM="FoM", D="D")[stat]
method <- match.arg(toupper(method), choices=c("LINEAR", "SPLINE"))
nPerGroup <- as.integer(nPerGroup[1])
nGroups <- as.integer(nGroups[1])
stopifnot(nPerGroup > 1L, nPerGroup <= max(shotGroups::DFdistr$n),
nGroups > 0L, nGroups <= max(shotGroups::DFdistr$nGroups))
sigma <- sigma[1]
stopifnot(is.numeric(sigma), sigma > 0)
if(!lower.tail) { p <- 1-p}
p_var <- paste0(stat, "_Q", sprintf("%.03d", round(1000*p)))
qq <- setNames(rep(NA_real_, length(p)), p_var)
idxGroup <- shotGroups::DFdistr$nGroups == nGroups
idxN <- shotGroups::DFdistr$n == nPerGroup
haveN <- unique(shotGroups::DFdistr$n[idxGroup])
## in checking whether probability is in lookup table, use 1 more decimal place
havePV <- paste0(stat, "_Q", sprintf("%.04d", trunc(10000*p))) %in%
paste0(names(shotGroups::DFdistr), "0")
# hasName(shotGroups::DFdistr, p_var)
if(sum(idxGroup) < 1L) {
warning("Lookup table does not have quantile(s) for nGroups=", nGroups)
} else {
if(all(havePV)) {
## all probabilities available
qq <- if(nPerGroup %in% haveN) {
## nPerGroup available
c(data.matrix(shotGroups::DFdistr[idxN & idxGroup, p_var[havePV], drop=FALSE]))
} else {
## interpolate nPerGroup
warning("Quantile(s) based on monotone spline interpolation for nPerGroup")
vapply(p_var[havePV], function(pvo) {
splinefun(haveN, shotGroups::DFdistr[idxGroup, pvo], method="monoH.FC")(nPerGroup) },
FUN.VALUE=numeric(1))
}
} else {
## interpolate p, but only within available min/max prob
## available probabilities
allPV <- regmatches(names(shotGroups::DFdistr),
regexpr(paste0("^", stat, "_Q[[:digit:]]{3}$"),
names(shotGroups::DFdistr)))
allP <- as.numeric(sub(paste0(stat, "_Q([[:digit:]]{3})$"), "\\1" , allPV)) / 1000
## make sure probabilities are sorted
allPV <- allPV[order(allP)]
allP <- allP[ order(allP)]
## check range of p is within limits
p_mm <- range(allP)
keep <- (p >= p_mm[1]) & (p <= p_mm[2])
if(nPerGroup %in% haveN) {
warning("Quantile(s) based on monotone spline interpolation for p")
## nPerGroup available
## get all quantiles in correct order
allQ <- shotGroups::DFdistr[idxN & idxGroup, allPV]
## interpolate
qq[keep] <- splinefun(allP, shotGroups::DFdistr[idxN & idxGroup, allPV],
method="monoH.FC")(p[keep])
} else {
## bivariate interpolation over support grid
grid_npq <- expand.grid(n=haveN, p=allP)
allQ <- vapply(haveN, function(n) {
qRangeStat(allP, sigma=sigma,
nPerGroup=n,
nGroups=nGroups,
stat=stat,
lower.tail=lower.tail) },
FUN.VALUE=numeric(length(allP)))
grid_npq[["q"]] <- c(t(allQ))
if(method == "LINEAR") {
if(requireNamespace("interp", quietly=TRUE)) {
warning("Quantile(s) based on bilinear interpolation for p and nPerGroup")
qq[keep] <- interp::interp(grid_npq[["n"]], grid_npq[["p"]], grid_npq[["q"]],
xo=rep(nPerGroup, length(p[keep])),
yo=p[keep],
linear=TRUE, extrap=FALSE, output="points")[["z"]]
} else {
warning("Install package 'interp' to enable bilinear interpolation\\n
for p and nPerGroup")
}
} else if(method == "SPLINE") {
if(requireNamespace("MBA", quietly=TRUE)) {
warning("Quantile(s) based on bivariate spline interpolation for p and nPerGroup")
surf <- MBA::mba.surf(grid_npq, no.X=500, no.Y=500)
idx_x <- which.min(abs(surf$xyz.est$x - nPerGroup))
idx_y <- vapply(p[keep], function(pp) {
which.min(abs(surf$xyz.est$y - pp)) }, FUN.VALUE=integer(1))
qq[keep] <- surf$xyz.est$z[cbind(rep(idx_x, length(idx_y)), idx_y)]
} else {
warning("Install package 'MBA' to enable bivariate spline interpolation\\n
for p and nPerGroup")
}
}
}
}
}
qq*sigma
}
## simulates range statistics for 1 group of shots with
## circular bivariate normal distribution
simRangeStat <-
function(nPerGroup, sigma) {
xy <- matrix(rnorm(2*nPerGroup, mean=0, sd=sigma), ncol=2)
H <- chull(xy) # convex hull indices (vertices ordered clockwise)
ES <- max(dist(xy[H, ], method="euclidean")) # extreme spread
bbW <- abs(diff(range(xy[ , 1])))
bbH <- abs(diff(range(xy[ , 2])))
FoM <- (bbW + bbH) / 2 # figure of merit
D <- sqrt(bbW^2 + bbH^2) # diagonal
c(ES=ES, FoM=FoM, D=D)
}
## random deviates for range statistics
rRangeStat <-
function(n, sigma=1, nPerGroup=5, nGroups=1, stat=c("ES", "FoM", "D")) {
stat <- match.arg(toupper(stat), choices=c("ES", "FOM", "D"))
stat <- c(ES="ES", FOM="FoM", D="D")[unique(stat)]
nPerGroup <- as.integer(nPerGroup[1])
nGroups <- as.integer(nGroups[1])
sigma <- sigma[1]
n <- if(length(n) > 1L) { length(n) } else { n }
stopifnot(is.numeric(sigma), sigma > 0)
stopifnot(nPerGroup > 1L, nPerGroup <= max(shotGroups::DFdistr$n),
nGroups > 0L, nGroups <= max(shotGroups::DFdistr$nGroups))
simOneRangeStat <- function() {
rs <- vapply(integer(nGroups), function(x) {
simRangeStat(nPerGroup, sigma) }, FUN.VALUE=numeric(3))
rowMeans(rs)[stat]
}
unname(replicate(n, simOneRangeStat()))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.