.Rhistory

## visualize true and estimated precision matrices
par(mfrow=c(1,3), pty="s")
image(pracma::flipud(Omega0),main="Original Precision")
image(pracma::flipud(Omega1$C), main="banded3::Bonferroni")
image(pracma::flipud(Omega2$C), main="banded3::Holm")
Omega1$optk
Omega2$optk
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
p = 200; n = 100
k0 = 5; A0min=0.1; A0max=0.2; D0min=2; D0max=5
set.seed(123)
A0 = matrix(0, p,p)
for(i in 2:p){
vals = runif(n=min(k0, i-1), min=A0min, max=A0max) * sample(c(1,-1), size=min(k0, i-1), replace = TRUE)
vals = vals[ order(abs(vals)) ]
A0[i, max(1, i-k0):(i-1)] = vals
}
D0 = diag(runif(n = p, min = D0min, max = D0max))
Omega0 = t(diag(p) - A0)%*%diag(1/diag(D0))%*%(diag(p) - A0)
## data generation (based on AR representation)
## it is same with generating n random samples from N_p(0, Omega0^{-1})
X = matrix(0, nrow=n, ncol=p)
X[,1] = rnorm(n, sd = sqrt(D0[1,1]))
for(j in 2:p){
mean.vec.j = X[, 1:(j-1)]%*%as.matrix(A0[j, 1:(j-1)])
X[,j] = rnorm(n, mean = mean.vec.j, sd = sqrt(D0[j,j]))
}
## banded estimation using two different schemes
Omega1 <- PreEst.banded3(X, upperK=20, algorithm="Bonferroni")
Omega2 <- PreEst.banded3(X, upperK=20, algorithm="Holm")
## visualize true and estimated precision matrices
par(mfrow=c(1,3), pty="s")
image(pracma::flipud(Omega0), main="Original Precision")
image(pracma::flipud(Omega1$C), main="banded3::Bonferroni")
image(pracma::flipud(Omega2$C), main="banded3::Holm")
Omega1$optk
Omega2$optk
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
?CovDist
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
help("CovTest1")
jay1 <- function(X){
n = nrow(X)
p = ncol(X)
X = scale(X, center=TRUE, scale=FALSE)
hXiXj = 0
for(j in 2:n){
for(i in 1:(j-1)){
hXiXj = hXiXj + (sum(X[i,]*X[j,]))^2 - ( sum(X[i,]^2) + sum(X[j,]^2) ) + p
}
}
Tn = hXiXj*2 / (n*(n-1)) # test statistic
return(Tn)
}
ks1 <- function(X, Sigma0){
if (missing(Sigma0)){
Sigma0 = diag(ncol(X))
}
n = nrow(X)
p = ncol(X)
X = scale(X, center=TRUE, scale=FALSE)
X = X%*%tmpinvroot(Sigma0)
hXiXj = 0
for(j in 2:n){
for(i in 1:(j-1)){
hXiXj = hXiXj + (sum(X[i,]*X[j,]))^2 - ( sum(X[i,]^2) + sum(X[j,]^2) ) + p
}
}
Tn = hXiXj*2 / (n*(n-1)) # test statistic
return(Tn)
}
tmpinvroot <- function(S){
eigS = eigen(S)
output = eigS$vectors %*% (diag(eigS$values^(-0.5))) %*% t(eigS$vectors)
return(output)
}
library(MASS)
n = 100
p = 200
Sigma0 = diag(p)
X = mvrnorm(n=n, mu=rep(0,p), Sigma=Sigma0)
CovTest1(X, Sigma0=diag(p), alpha=0.01)
CovTest1(X, Sigma0=diag(p), alpha=0.05)
jay1(X)
library(CovTools)
jay1(X)
CovTest1(X, Sigma0=diag(p), alpha=0.05)
library(CovTools)
library(CovTools)
CovTest1(X, Sigma0=diag(p), alpha=0.05)
jay1(X)
library(CovTools)
library(CovTools)
library(CovTools)
q()
install.packages(c("shapes","mvtnorm","pracma"))
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
A = matrix(rnorm(10*3),nrow=10)
X = A
p = ncol(X)
n = nrow(X)
log.BF.mat = matrix(0, ncol=p,nrow=p) # log Bayes factors
for(i in 1:p){
Xi = matrix(X[,i], ncol=1)
for(j in (1:p)[-i]){
Xj = matrix(X[,j], ncol=1)
log.BF.mat[i,j] = a0*log(b0) - lgamma(a0) +
1/2 * log(gamma/(1+gamma)) + lgamma(n/2 + a0) +
1/2 * sum((Xi)^2) - (n/2 + a0) * log(1/2 * ( sum((Xi)^2) - sum(Xi*Xj)^2/sum((Xj)^2) /(1+gamma) ) + b0)
}
}
a0 = 1
b0 = 1
p = ncol(X)
n = nrow(X)
log.BF.mat = matrix(0, ncol=p,nrow=p) # log Bayes factors
for(i in 1:p){
Xi = matrix(X[,i], ncol=1)
for(j in (1:p)[-i]){
Xj = matrix(X[,j], ncol=1)
log.BF.mat[i,j] = a0*log(b0) - lgamma(a0) +
1/2 * log(gamma/(1+gamma)) + lgamma(n/2 + a0) +
1/2 * sum((Xi)^2) - (n/2 + a0) * log(1/2 * ( sum((Xi)^2) - sum(Xi*Xj)^2/sum((Xj)^2) /(1+gamma) ) + b0)
}
}
gamma = 1
p = ncol(X)
n = nrow(X)
log.BF.mat = matrix(0, ncol=p,nrow=p) # log Bayes factors
for(i in 1:p){
Xi = matrix(X[,i], ncol=1)
for(j in (1:p)[-i]){
Xj = matrix(X[,j], ncol=1)
log.BF.mat[i,j] = a0*log(b0) - lgamma(a0) +
1/2 * log(gamma/(1+gamma)) + lgamma(n/2 + a0) +
1/2 * sum((Xi)^2) - (n/2 + a0) * log(1/2 * ( sum((Xi)^2) - sum(Xi*Xj)^2/sum((Xj)^2) /(1+gamma) ) + b0)
}
}
log.BF.mat
diag(log.BF.mat) = log.BF.mat[1,2]
log.BF.mat
install.packages(c("MCMCpack","doMC","MASS"))
install.packages(c("MCMCpack", "doMC", "MASS"))
######################################################################
# R code for one-sample covariance test in Lee, Lin and Dunson (2018).
#
# References
# - Kyoungjae Lee, Lizhen Lin and David Dunson. (2018) Maximum Pairwise Bayes Factors for Covariance Structure Testing. [https://arxiv.org/abs/1809.03105]
# - Tony Cai and Zongming Ma. (2013) Optimal hypothesis testing for high dimensional covariance matrices.
# - Muni Srivastava, Hirokazu Yanagihara and Tatsuya Kubokawa. (2014) Tests for covariance matrices in high dimension with less sample size.
#
# made by Kyoungjae Lee.
# version 1 (Sep. 09, 2018)
######################################################################
##################################################################
# Functions
##################################################################
# mxPBF for one-sample covariance test
mxPBF <- function(X, a0, b0, gamma){
###############################################################
# Input
#     X: n times p data matrix
#     a0, b0, gamma: hyperparameters for mxPBF
#
# Output
#     log.BF.mat: p times p matrix of pairwise log Bayes factors
#                 The (i,j)th entry of log.BF.mat is the pairwise log Bayes factor,
#                 log BF_{10}(X_i, X_j).
#
# Description:
#     Let X_i be the i-th column of X.
#     [Model]
#        X_i | X_j ~ N_n( a_{ij}X_j, \tau_{ij}^2 I_n )
#     [Hypothesis testing]
#        H_0: a_{ij}=0 & \tau_{ij}=1  versus  H_1: not H_0
#     [Prior]
#        under H_1: a_{ij}|\tau_{ij}^2 ~ N(0, \tau_{ij}^2/[gamma*||X_j||^2])
#                   \tau_{ij}^2 ~ IG(a0, b0)
###############################################################
p = ncol(X)
n = nrow(X)
log.BF.mat = matrix(0, ncol=p,nrow=p) # log Bayes factors
for(i in 1:p){
Xi = matrix(X[,i], ncol=1)
for(j in (1:p)[-i]){
Xj = matrix(X[,j], ncol=1)
log.BF.mat[i,j] = a0*log(b0) - lgamma(a0) +
1/2 * log(gamma/(1+gamma)) + lgamma(n/2 + a0) +
1/2 * sum((Xi)^2) - (n/2 + a0) * log(1/2 * ( sum((Xi)^2) - sum(Xi*Xj)^2/sum((Xj)^2) /(1+gamma) ) + b0)
}
}
diag(log.BF.mat) = log.BF.mat[1,2] # just to fill out the diagonal parts
return(log.BF.mat)
}
# Cai and Ma (2013): one-sample covariance test
CovTest1 <- function(X, Sigma0 = diag(p), alpha = 0.05){
p = ncol(X)
n = nrow(X)
z.val = qnorm(1 - alpha)
hXiXj = 0
for(j in 2:n){
for(i in 1:(j-1)){
hXiXj = hXiXj + (sum(X[i,]*X[j,]))^2 - ( sum(X[i,]^2) + sum(X[j,]^2) ) + p
}
}
Tn = hXiXj*2 / (n*(n-1)) # test statistic
res = list()
res$Tn = Tn
res$test = ( Tn > z.val * 2 * sqrt(p*(p+1) /(n*(n-1))) )
return( res )
}
# Srivastava et al. (2014): one-sample covariance test
CovTest2 <- function(X, Sigma0 = diag(p), alpha = 0.05){
p = ncol(X)
n = nrow(X)
z.val = qnorm(1 - alpha)
bar.X = matrix(colMeans(X), nrow = p, ncol = 1)
Sn = cov(X)
Yn = X - matrix(c(bar.X), nrow=n, ncol=p, byrow=T)
a1.hat = sum(diag(Sn)) / p
a2.hat = ( (n-1)*(n-2)*sum(diag( (n-1)^2*Sn%*%Sn )) - n*(n-1)*sum((rowSums(Yn^2))^2) + (sum(diag( (n-1)*Sn )))^2 ) /( p*n*(n-1)*(n-2)*(n-3) )
T2 = (n-1)/2 * ( a2.hat - 2*a1.hat + 1 )
res = list()
res$T2 = T2
res$test = ( T2 > z.val )
return( res )
}
# auxiliary function for tridiagonal matrix
tridiag <- function(upper.val, lower.val, diag.val){
p = length(diag.val)
res = matrix(0, p, p)
ind = 1:(p-1)
res[ cbind(ind+1, ind) ] = lower.val
res[ cbind(ind, ind+1) ] = upper.val
diag(res) = diag.val
return(res)
}
library(MCMCpack)
library(MASS)
library(doMC)
##################################################################
# Simulations
##################################################################
n.sim = 100 # number of simulations
p = 200 # number of variables
n = 100 # number of observations
registerDoMC(1) # number of cores to use
# registerDoMC(20)
##########################################
# Under H0
##########################################
tot.res0 <- foreach (sim = 1:n.sim) %dopar% {
tot.res.mat = matrix(0, nrow = 1, ncol = 3)
set.seed(sim)
Sigma0 = diag(p)
X = mvrnorm(n = n, mu = rep(0, p), Sigma = Sigma0) # data generation
# hyperparameters for mxPBF
a0 = 2; b0 = 2; c_gam = 10; alpha = 8.01*(1-1/log(n))
gamma = c_gam*max(n,p)^{-alpha}
log.BF.mat = mxPBF(X, a0, b0, gamma)
tot.res.mat[, 1] = max(log.BF.mat) # mxPBF
tot.res.mat[, 2] = CovTest1(X)$Tn # Cai and Ma's test
tot.res.mat[, 3] = CovTest2(X)$T2 # Srivastava et al.'s test
return(tot.res.mat)
}# data.ind for loop
##########################################
# Under H1
##########################################
# magnitude of signals in H1
rho.vec = seq(from = 0.3, to = 0.8, by = 0.025) # (1) tridiagonal
# rho.vec = seq(from = 0.1, to = 0.5, by = 0.025) # (2) sparse signals
rho.length = length(rho.vec)
tot.res1 <- foreach (sim = 1:n.sim) %dopar% {
test.mat = matrix(0, nrow = rho.length, ncol = 3)
for(rho.ind in 1:rho.length){
set.seed(sim)
# under H1
rho = rho.vec[rho.ind]
Sigma0 = diag(p)
Sigma0 = tridiag(rep(rho, p-1), rep(rho, p-1), rep(1, p)) # (1) tridiagonal
# Sigma0[1,2] = rho; Sigma0[2,1] = rho # (2) Sparse signals
X = mvrnorm(n = n, mu = rep(0, p), Sigma = Sigma0) # data generation
# hyperparameters for mxPBF
a0 = 2; b0 = 2; c_gam = 10; alpha = 8.01*(1-1/log(n))
gamma = c_gam*max(n,p)^{-alpha}
log.BF.mat = mxPBF(X, a0, b0, gamma)
test.mat[rho.ind, 1] = max(log.BF.mat) # mxPBF
test.mat[rho.ind, 2] = CovTest1(X)$Tn # Cai and Ma's test
test.mat[rho.ind, 3] = CovTest2(X)$T2 # Srivastava et al.'s test
} # rho.ind for loop
return(test.mat)
}# data.ind for loop
library(CovTools)
data = mvtnorm::rmvnorm(100, sigma=diag(5))
## run test
CovTest1(data, Sigma0=diag(5))
out = CovTest1(data, Sigma0=diag(5))
out
library(CovTools)
out = CovTest1(data, Sigma0=diag(5))
lout
out
test1.Cai13 <- function(X, Sigma0, alpha){
# 1. setup
scaler = get_invroot(Sigma0)
# 2. mean centering
X.centered = scale(X, center=TRUE, scale=FALSE)
X.centered = matrix(X.centered,nrow(X)) %*% scaler
# 3. parameters and hXiXj
n = nrow(X)
p = ncol(X)
hXiXj = rcpptest1_cai11(X.centered)
# 4. post-adjusting
statistic = (hXiXj*2)/(n*(n-1)) # test statistic
threshold = (qnorm(1-alpha))*2*sqrt((p*(p+1))/(n*(n-1)))
reject = (statistic > threshold)
# 5. return
return(list(statistic=statistic, threshold=threshold, reject=reject))
}
test1.Cai13(data, diag(ncol(data)), 0.05)
rm(list=ls())
data = mvtnorm::rmvnorm(100, sigma=diag(5))
## run test
CovTest1(data, Sigma0=diag(5))
BayesTest1(data)
rm(list=ls())
library(CovTools)
data = matrix(rnorm(100*4), nrow=100)
output = BayesTest1(data)
output
library(CovTools)
library(CovTools)
output = BayesTest1(data)
output
library(CovTools)
library(CovTools)
library(CovTools)
library(CovTools)
usethis::use_news_md()
library(CovTools)
library(CovTools)
citation("DAAG")
install.packages(c("DAAG","Hotelling"))
library(Hotelling)
citation(Hotelling)
citation("Hotelling")
library(DAAG)
citation("DAAG")
library(RcppArmadillo)
citation("RcppArmadillo")
library(Rcpp)
citation("Rcpp")
print("Rcpp",bibtex=TRUE)
library(RcppXPtrUtils)
help("RcppXPtrUtils")
citation("RcppXPtrUtils")
citation("RcppArmadillo")
help("DAAG")
citation("DAAG")
citation("CovTools")
help(t.test)
install.packages(c("geigen","shapes","expm","pracma","Rdpack","SHT"), dependencies = TRUE)
library(CovTools)
library(CovTools)
library(CovTools)
# set a seed for reproducibility
set.seed(11)
#  small data with identity covariance
dat.small <- matrix(rnorm(20*5), ncol=5)
#  run the code with highly structured estimator
out.small <- CovEst.2010RBLW(dat.small)
#  visualize
par(mfrow=c(1,3), pty="s")
install.packages(c("devtools","roxygen2"), dependencies = TRUE)
library(CovTools)
# set a seed for reproducibility
set.seed(11)
#  small data with identity covariance
dat.small <- matrix(rnorm(20*5), ncol=5)
#  run the code with highly structured estimator
out.small <- CovEst.2010RBLW(dat.small)
#  visualize
par(mfrow=c(1,3), pty="s")
image(diag(5), main="true cov")
image(cov(dat.small), main="sample cov")
image(out.small$S, main="estimated cov")
library(CovTools)
# set a seed for reproducibility
set.seed(11)
#  small data with identity covariance
dat.small <- matrix(rnorm(20*5), ncol=5)
#  run the code
out.small <- CovEst.2010OAS(dat.small)
#  visualize
par(mfrow=c(1,3), pty="s")
image(diag(5), main="true cov")
image(cov(dat.small), main="sample cov")
image(out.small$S, main="estimated cov")
\dontrun{
## want to see how delta is determined according to
#  the number of observations we have.
nsamples = seq(from=5, to=200, by=5)
nnsample = length(nsamples)
#  we will record two values; rho and norm difference
vec.rho   = rep(0, nnsample)
vec.normd = rep(0, nnsample)
for (i in 1:nnsample){
dat.norun <- matrix(rnorm(nsamples[i]*5), ncol=5) # sample in R^5
out.norun <- CovEst.2010OAS(dat.norun)            # run with default
vec.rho[i]   = out.norun$delta
vec.normd[i] = norm(out.norun$S - diag(5),"f")    # Frobenius norm
}
## want to see how delta is determined according to
#  the number of observations we have.
nsamples = seq(from=5, to=200, by=5)
nnsample = length(nsamples)
#  we will record two values; rho and norm difference
vec.rho   = rep(0, nnsample)
vec.normd = rep(0, nnsample)
for (i in 1:nnsample){
dat.norun <- matrix(rnorm(nsamples[i]*5), ncol=5) # sample in R^5
out.norun <- CovEst.2010RBLW(dat.norun)           # run with default
vec.rho[i]   = out.norun$rho
vec.normd[i] = norm(out.norun$S - diag(5),"f")    # Frobenius norm
}
# let's visualize the results
par(mfrow=c(1,2))
plot(nsamples, vec.rho,   lwd=2, type="b", col="red", main="estimated rhos")
plot(nsamples, vec.normd, lwd=2, type="b", col="blue",main="Frobenius error")
library(CovTools)
## CRAN-purpose small computation
# set a seed for reproducibility
set.seed(11)
#  small data with identity covariance
dat.small <- matrix(rnorm(20*5), ncol=5)
#  run the code
out.small <- CovEst.2010OAS(dat.small)
#  visualize
par(mfrow=c(1,3), pty="s")
image(diag(5), main="true cov")
image(cov(dat.small), main="sample cov")
image(out.small$S, main="estimated cov")
out.small$S
dat.small <- matrix(rnorm(5*10), ncol=10)
#  run the code
out.small <- CovEst.2010OAS(dat.small)
#  visualize
par(mfrow=c(1,3), pty="s")
image(diag(10), main="true cov")
image(cov(dat.small), main="sample cov")
image(out.small$S, main="estimated cov")
dat.small <- matrix(rnorm(10*5), ncol=10)
#  run the code
out.small <- CovEst.2010RBLW(dat.small)
#  visualize
par(mfrow=c(1,3), pty="s")
image(diag(10), main="true cov")
image(cov(dat.small), main="sample cov")
image(out.small$S, main="estimated cov")
out.small$rho
# set a seed for reproducibility
set.seed(11)
#  small data with identity covariance
dat.small <- matrix(rnorm(5*10), ncol=10)
#  run the code
out.small <- CovEst.2010OAS(dat.small)
#  visualize
par(mfrow=c(1,3), pty="s")
image(diag(10), main="true cov")
image(cov(dat.small), main="sample cov")
image(out.small$S, main="estimated cov")
out.small$S
out.small$rho
library(CovTools)
library(CovTools)
library(CovTools)
A = matrix(rnorm(5*10),ncol=10)
aA = abs(A)
rowSums(aA)-diag(aA)
B = aA
diag(B) = 0
rowSums(B)
diag(aA)
library(CovTools)
help("BDiagTest1.mxPBF")
library(CovTools)
help("BDiagTest1.mxPBF")
library(CovTools)
help("BDiagTest1.mxPBF")
library(CovTools)
library(CovTools)
library(CovTools)