Nothing
tests/groupTsts.R
In hdi: High-Dimensional Inference
#### Testing groupBound() , clusterGroupBound() , etc
#### ========== =================
## Examples in ../man/groupBound.rd and ../man/clusterGroupBound.rd
## should *not*
## be duplicated here. OTOH, want strict testing (*and* be speedy !)
stopifnot(require(hdi))
## From Matrix/.../tests-tools-1.R :
showSys.time <- function(expr, ...) {
## prepend 'Time' for R CMD Rdiff
st <- system.time(expr, ...)
writeLines(paste("Time", capture.output(print(st))))
invisible(st)
}
doExtras <- function ()
{
interactive() || nzchar(Sys.getenv("R_hdi_check_extra")) ||
identical("true", unname(Sys.getenv("R_PKG_CHECKING_doExtras")))
}
do.Extras <- doExtras()
b64 <- .Machine$sizeof.pointer == 8
cat(if(b64) "64" else "32", "bit platform; do.Extras = ", do.Extras,"\n")
##' @title Generate "Reference" Datasets where variables are in *groups*
##' @param pvec integer vector of number of variables p[k], in k-th group, k = 1..np.
##' @param s0.v integer vector of non-zero coefficient s0[k], in k-th group, k = 1..np.
##' @param n number of observations
##' @param ... further arguments such as \code{xtype}, \code{btype}, passed to
##' \code{\link{rXb}}
##' @param iteration see \code{\link{rXb}}
##' @param do2S ditto
##' @return a \code{\link{list}} with components \code{x} and \code{beta} similar to those
##' of \code{rXb()} and \code{p.vec} and \code{s0.vec} which are the
##' corresponding arguments.
##' @author Martin Maechler
rXbBlock <- function(p.vec, s0.vec, n, ..., iteration = NA, do2S = FALSE) {
stopifnot(is.numeric(p.vec), p.vec == round(p.vec),
is.numeric(s0.vec), s0.vec == round(s0.vec),
(np <- length(p.vec)) >= 1, length(s0.vec) == np,
0 <= s0.vec, s0.vec <= p.vec)
grL <-
lapply(1:np, function(k)
rXb(n = n, p = p.vec[k], s0 = s0.vec[k],
..., iteration=iteration, do2S=do2S))
## now concatenate the X matrices and the beta vectors:
list(x = do.call(cbind, lapply(grL, `[[`, "x")),
beta = do.call(c, lapply(grL, `[[`, "beta")),
p.vec = p.vec, s0.vec = s0.vec)
}
##' @title From (X, beta, sigma), generate y = X b + eps
##' @param rXb a list, e.g. resulting from rXb (== rXb) or rXbBlock()
##' @param sigma a non-negative number
##' @return a list with components (x, y)
##' @author Martin Maechler
rXy <- function(rXb, sigma = 1) {
stopifnot(is.list(rXb), is.matrix(X <- rXb$x), is.numeric(b <- rXb$beta),
ncol(X) == length(b),
length(sigma) == 1, is.finite(sigma), 0 <= sigma)
list(x = X, y = drop(X %*% b) + rnorm(nrow(X), sd = sigma))
}
##' from a p-vector, return corresponding list of 'groups', i.e., group indices
pvec2grp <- function(pv) {
stopifnot(is.numeric(pv), is.finite(pv), pv >= 0)
cv <- cumsum(c(1,pv))
mapply(seq, from=cv[-length(cv)], to = cv[-1L]-1L)
}
if(FALSE) { ## not used currently
suppressWarnings(RNGversion("3.5.0"))
set.seed(47) # <-> iteration = NA to use R's seed
d1 <- rXb(n = 64, p = 16, s0 = 7, iteration=NA, do2S = FALSE)
d2 <- rXb(n = 64, p = 256, s0 = 12, iteration=NA, do2S = FALSE)
stopifnot(identical(dim(d2$x), c(64L, 256L)), sum(d2$beta != 0) == 12)
yhat.1 <- with(d1, drop(x %*% beta))
yhat.2 <- with(d2, drop(x %*% beta))
}
suppressWarnings(RNGversion("3.5.0"))
set.seed(20)
d1B <- rXbBlock(p = c(7, 8, 7), s0 = c(2, 4, 1), n = 16)
p2 <- sample(c(32, 16, 16, 16, 12, 12, 5, 5, 5))
s2 <- c(rmultinom(1, size = p2, prob = ifelse(p2 > 15, 1/4, 1/8)))
rbind(p2, s2)
d2B <- rXbBlock(p = p2, s0 = s2, n = 500)
d3B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "toeplitz", x.par = 0.95)
d4B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "equi.corr")
d5B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "exp.decay", x.par = c(0.25, 4))
if(interactive() && exists("db5", mode="list")) ## save when produced with "new hdi"
save(d1B, d2B, d3B, d4B, d5B, file = "rXb_blk.rda")
if(!exists("rXb", mode="function")) ## "old hdi" -- tests should give the same (?!) -- TODO
print(load("rXb_blk.rda"))
if(dev.interactive(TRUE)) {
require(Matrix)
## if you look carefully, you recognize the three blocks of sizes (7, 8, 7):
image(as(zapsmall(cov(d1B$x), digits=1), "sparseMatrix"))
## Here, have 9 blocks :
image(as(zapsmall(cov(d2B$x), digits=1), "sparseMatrix"))
image(as(zapsmall(cov(d3B$x), digits=1), "sparseMatrix")) # easy to see
image(as(zapsmall(cov(d4B$x), digits=1), "sparseMatrix")) # (ditto)
image(as(zapsmall(cov(d5B$x), digits=1), "sparseMatrix"))
# decaying quickly => basically tri-diagonal blocks
}
##- set.seed(107)
##- d1xy <- rXy(d1B, sigma = 0.5)
##- showSys.time(
##- gb1 <- groupBound(x = d1xy$x, y = d1xy$y, group = 1:ncol(d1xy$x))
##- ) # 2.00 [lynne 2015]
##- set.seed(17)# also for random splits
##- d2xy <- rXy(d2B, sigma = 1/16)
##- showSys.time(
##- gb2 <- groupBound(x = d2xy$x, y = d2xy$y, group = 1:ncol(d2xy$x))
##- )
##- showSys.time(
##- gb2g <- groupBound(x = d2xy$x, y = d2xy$y, group = pvec2grp(d2B$p.vec),
##- silent = TRUE)
##- )
##-
##- stopifnot(
##- as.vector(gb1) == 0
##- ,
##- all.equal(as.vector(gb2),
##- if(b64) 6.529515 else 6.598612, tol = 2e-7)
##- ,
##- all.equal(as.vector(gb2g),
##- if(b64) c(0, 0, 1.924578, 0, 1.436938, 0, 0, 0, 0.3264166)
##- else c(0, 0.01877908, 2.0447942, 0, 1.753219, 0,0,0, 0.3779607),
##- tol = 6e-7)
##- )
if(do.Extras) { ## save time in regular 'R CMD check' ---------
##--------------- 3 --
suppressWarnings(RNGversion("3.5.0"))
set.seed(27)# also for random splits
d3xy <- rXy(d3B, sigma = 1/16)
showSys.time(
gb3 <- groupBound(x = d3xy$x, y = d3xy$y, group = 1:ncol(d3xy$x))
)
stopifnot(
all.equal(as.vector(gb3), if(b64) 7.3760622 else 7.7908584, tol = 2e-7)
)
suppressWarnings(RNGversion("3.5.0"))
set.seed(19)# for random splits
showSys.time(
gb3g <- groupBound(x = d3xy$x, y = d3xy$y, group = pvec2grp(d3B$p.vec),
silent = TRUE)
)
stopifnot(
all.equal(as.vector(gb3g),
if(b64) c(0, 0.5237811, 1.354151, 0, 0, 0, 0.8847217, 0, 0)
else c(0, 1.1008886, 1.9249053,0, 0, 0, 0.78661001, 0, 0),
tol = 5e-7) # 9.4e-8
)
##--------------- 4 --
suppressWarnings(RNGversion("3.5.0"))
set.seed(37)# also for random splits
d4xy <- rXy(d4B, sigma = 1/16)
showSys.time(
gb4 <- groupBound(x = d4xy$x, y = d4xy$y, group = 1:ncol(d4xy$x),
silent = TRUE)
)
showSys.time(
gb4g <- groupBound(x = d4xy$x, y = d4xy$y, group = pvec2grp(d4B$p.vec),
silent = TRUE)
)
stopifnot(
all.equal(as.vector(gb4), if(b64) 5.230821 else 5.3752713, tol = 2e-7)
## all.equal(as.vector(gb4), 8.475868, tol = 2e-7)
,
all.equal(as.vector(gb4g),
if(b64) c(0, 0, 3.328455, 0, 0, 0, 0, 0, 0)
else c(0, 0.043785673, 3.6154918, 0, 0, 0, 0, 0, 0),
tol = 8e-7) # 1.4e-7
)
}# only if(do.Extras)
Try the hdi package in your browser
Any scripts or data that you put into this service are public.
hdi documentation built on May 27, 2021, 3 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#### Testing groupBound() , clusterGroupBound() , etc
#### ========== =================
## Examples in ../man/groupBound.rd and ../man/clusterGroupBound.rd
## should *not*
## be duplicated here. OTOH, want strict testing (*and* be speedy !)
stopifnot(require(hdi))
## From Matrix/.../tests-tools-1.R :
showSys.time <- function(expr, ...) {
## prepend 'Time' for R CMD Rdiff
st <- system.time(expr, ...)
writeLines(paste("Time", capture.output(print(st))))
invisible(st)
}
doExtras <- function ()
{
interactive() || nzchar(Sys.getenv("R_hdi_check_extra")) ||
identical("true", unname(Sys.getenv("R_PKG_CHECKING_doExtras")))
}
do.Extras <- doExtras()
b64 <- .Machine$sizeof.pointer == 8
cat(if(b64) "64" else "32", "bit platform; do.Extras = ", do.Extras,"\n")
##' @title Generate "Reference" Datasets where variables are in *groups*
##' @param pvec integer vector of number of variables p[k], in k-th group, k = 1..np.
##' @param s0.v integer vector of non-zero coefficient s0[k], in k-th group, k = 1..np.
##' @param n number of observations
##' @param ... further arguments such as \code{xtype}, \code{btype}, passed to
##' \code{\link{rXb}}
##' @param iteration see \code{\link{rXb}}
##' @param do2S ditto
##' @return a \code{\link{list}} with components \code{x} and \code{beta} similar to those
##' of \code{rXb()} and \code{p.vec} and \code{s0.vec} which are the
##' corresponding arguments.
##' @author Martin Maechler
rXbBlock <- function(p.vec, s0.vec, n, ..., iteration = NA, do2S = FALSE) {
stopifnot(is.numeric(p.vec), p.vec == round(p.vec),
is.numeric(s0.vec), s0.vec == round(s0.vec),
(np <- length(p.vec)) >= 1, length(s0.vec) == np,
0 <= s0.vec, s0.vec <= p.vec)
grL <-
lapply(1:np, function(k)
rXb(n = n, p = p.vec[k], s0 = s0.vec[k],
..., iteration=iteration, do2S=do2S))
## now concatenate the X matrices and the beta vectors:
list(x = do.call(cbind, lapply(grL, `[[`, "x")),
beta = do.call(c, lapply(grL, `[[`, "beta")),
p.vec = p.vec, s0.vec = s0.vec)
}
##' @title From (X, beta, sigma), generate y = X b + eps
##' @param rXb a list, e.g. resulting from rXb (== rXb) or rXbBlock()
##' @param sigma a non-negative number
##' @return a list with components (x, y)
##' @author Martin Maechler
rXy <- function(rXb, sigma = 1) {
stopifnot(is.list(rXb), is.matrix(X <- rXb$x), is.numeric(b <- rXb$beta),
ncol(X) == length(b),
length(sigma) == 1, is.finite(sigma), 0 <= sigma)
list(x = X, y = drop(X %*% b) + rnorm(nrow(X), sd = sigma))
}
##' from a p-vector, return corresponding list of 'groups', i.e., group indices
pvec2grp <- function(pv) {
stopifnot(is.numeric(pv), is.finite(pv), pv >= 0)
cv <- cumsum(c(1,pv))
mapply(seq, from=cv[-length(cv)], to = cv[-1L]-1L)
}
if(FALSE) { ## not used currently
suppressWarnings(RNGversion("3.5.0"))
set.seed(47) # <-> iteration = NA to use R's seed
d1 <- rXb(n = 64, p = 16, s0 = 7, iteration=NA, do2S = FALSE)
d2 <- rXb(n = 64, p = 256, s0 = 12, iteration=NA, do2S = FALSE)
stopifnot(identical(dim(d2$x), c(64L, 256L)), sum(d2$beta != 0) == 12)
yhat.1 <- with(d1, drop(x %*% beta))
yhat.2 <- with(d2, drop(x %*% beta))
}
suppressWarnings(RNGversion("3.5.0"))
set.seed(20)
d1B <- rXbBlock(p = c(7, 8, 7), s0 = c(2, 4, 1), n = 16)
p2 <- sample(c(32, 16, 16, 16, 12, 12, 5, 5, 5))
s2 <- c(rmultinom(1, size = p2, prob = ifelse(p2 > 15, 1/4, 1/8)))
rbind(p2, s2)
d2B <- rXbBlock(p = p2, s0 = s2, n = 500)
d3B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "toeplitz", x.par = 0.95)
d4B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "equi.corr")
d5B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "exp.decay", x.par = c(0.25, 4))
if(interactive() && exists("db5", mode="list")) ## save when produced with "new hdi"
save(d1B, d2B, d3B, d4B, d5B, file = "rXb_blk.rda")
if(!exists("rXb", mode="function")) ## "old hdi" -- tests should give the same (?!) -- TODO
print(load("rXb_blk.rda"))
if(dev.interactive(TRUE)) {
require(Matrix)
## if you look carefully, you recognize the three blocks of sizes (7, 8, 7):
image(as(zapsmall(cov(d1B$x), digits=1), "sparseMatrix"))
## Here, have 9 blocks :
image(as(zapsmall(cov(d2B$x), digits=1), "sparseMatrix"))
image(as(zapsmall(cov(d3B$x), digits=1), "sparseMatrix")) # easy to see
image(as(zapsmall(cov(d4B$x), digits=1), "sparseMatrix")) # (ditto)
image(as(zapsmall(cov(d5B$x), digits=1), "sparseMatrix"))
# decaying quickly => basically tri-diagonal blocks
}
##- set.seed(107)
##- d1xy <- rXy(d1B, sigma = 0.5)
##- showSys.time(
##- gb1 <- groupBound(x = d1xy$x, y = d1xy$y, group = 1:ncol(d1xy$x))
##- ) # 2.00 [lynne 2015]
##- set.seed(17)# also for random splits
##- d2xy <- rXy(d2B, sigma = 1/16)
##- showSys.time(
##- gb2 <- groupBound(x = d2xy$x, y = d2xy$y, group = 1:ncol(d2xy$x))
##- )
##- showSys.time(
##- gb2g <- groupBound(x = d2xy$x, y = d2xy$y, group = pvec2grp(d2B$p.vec),
##- silent = TRUE)
##- )
##-
##- stopifnot(
##- as.vector(gb1) == 0
##- ,
##- all.equal(as.vector(gb2),
##- if(b64) 6.529515 else 6.598612, tol = 2e-7)
##- ,
##- all.equal(as.vector(gb2g),
##- if(b64) c(0, 0, 1.924578, 0, 1.436938, 0, 0, 0, 0.3264166)
##- else c(0, 0.01877908, 2.0447942, 0, 1.753219, 0,0,0, 0.3779607),
##- tol = 6e-7)
##- )
if(do.Extras) { ## save time in regular 'R CMD check' ---------
##--------------- 3 --
suppressWarnings(RNGversion("3.5.0"))
set.seed(27)# also for random splits
d3xy <- rXy(d3B, sigma = 1/16)
showSys.time(
gb3 <- groupBound(x = d3xy$x, y = d3xy$y, group = 1:ncol(d3xy$x))
)
stopifnot(
all.equal(as.vector(gb3), if(b64) 7.3760622 else 7.7908584, tol = 2e-7)
)
suppressWarnings(RNGversion("3.5.0"))
set.seed(19)# for random splits
showSys.time(
gb3g <- groupBound(x = d3xy$x, y = d3xy$y, group = pvec2grp(d3B$p.vec),
silent = TRUE)
)
stopifnot(
all.equal(as.vector(gb3g),
if(b64) c(0, 0.5237811, 1.354151, 0, 0, 0, 0.8847217, 0, 0)
else c(0, 1.1008886, 1.9249053,0, 0, 0, 0.78661001, 0, 0),
tol = 5e-7) # 9.4e-8
)
##--------------- 4 --
suppressWarnings(RNGversion("3.5.0"))
set.seed(37)# also for random splits
d4xy <- rXy(d4B, sigma = 1/16)
showSys.time(
gb4 <- groupBound(x = d4xy$x, y = d4xy$y, group = 1:ncol(d4xy$x),
silent = TRUE)
)
showSys.time(
gb4g <- groupBound(x = d4xy$x, y = d4xy$y, group = pvec2grp(d4B$p.vec),
silent = TRUE)
)
stopifnot(
all.equal(as.vector(gb4), if(b64) 5.230821 else 5.3752713, tol = 2e-7)
## all.equal(as.vector(gb4), 8.475868, tol = 2e-7)
,
all.equal(as.vector(gb4g),
if(b64) c(0, 0, 3.328455, 0, 0, 0, 0, 0, 0)
else c(0, 0.043785673, 3.6154918, 0, 0, 0, 0, 0, 0),
tol = 8e-7) # 1.4e-7
)
}# only if(do.Extras)
Try the hdi package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.