PreEst.2014An {CovTools}

Banded Precision Matrix Estimation via Bandwidth Test

Description

PreEst.2014An returns an estimator of the banded precision matrix using the modified Cholesky decomposition. It uses the estimator defined in Bickel and Levina (2008). The bandwidth is determined by the bandwidth test suggested by An, Guo and Liu (2014).

Usage

PreEst.2014An(
  X,
  upperK = floor(ncol(X)/2),
  algorithm = c("Bonferroni", "Holm"),
  alpha = 0.01
)

Arguments

X

an (n\times p) data matrix where each row is an observation.

upperK

upper bound of bandwidth k.

algorithm

bandwidth test algorithm to be used.

alpha

significance level for the bandwidth test.

Value

a named list containing:

C

a (p\times p) estimated banded precision matrix.

optk

an estimated optimal bandwidth acquired from the test procedure.

References

An B, Guo J, Liu Y (2014). “Hypothesis testing for band size detection of high-dimensional banded precision matrices.” Biometrika, 101(2), 477–483. ISSN 0006-3444, 1464-3510.

Bickel PJ, Levina E (2008). “Regularized estimation of large covariance matrices.” The Annals of Statistics, 36(1), 199–227. ISSN 0090-5364.

Examples

## Not run: 
## parameter setting
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){
  term1 = runif(n=min(k0,i-1),min=A0min, max=A0max)
  term2 = sample(c(1,-1),size=min(k0,i-1),replace=TRUE)
  vals  = term1*term2
  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.2014An(X, upperK=20, algorithm="Bonferroni")
Omega2 <- PreEst.2014An(X, upperK=20, algorithm="Holm")

## visualize true and estimated precision matrices
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
image(Omega0[,p:1],   main="Original Precision")
image(Omega1$C[,p:1], main="banded3::Bonferroni")
image(Omega2$C[,p:1], main="banded3::Holm")
par(opar)

## End(Not run)