PreEst.2014An: Banded Precision Matrix Estimation via Bandwidth Test

In CovTools: Statistical Tools for Covariance Analysis

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

\insertRef

an_hypothesis_2014CovTools

\insertRef

bickel_regularized_2008CovTools

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)

CovTools documentation built on Aug. 14, 2021, 1:08 a.m.