CovEst.2010RBLW {CovTools}

Rao-Blackwell Ledoit-Wolf Estimator

Description

Authors propose to estimate covariance matrix by minimizing mean squared error with the following formula,

\hat{\Sigma} = \rho \hat{F} + (1-\rho) \hat{S}

where \rho \in (0,1) a control parameter/weight, \hat{S} an empirical covariance matrix, and \hat{F} a target matrix. It is proposed to use a structured estimate \hat{F} = \textrm{Tr} (\hat{S}/p) \cdot I_{p\times p} where I_{p\times p} is an identity matrix of dimension p.

Usage

CovEst.2010RBLW(X)

Arguments

X

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

Value

a named list containing:

S

a (p\times p) covariance matrix estimate.

rho

an estimate for convex combination weight.

References

Chen Y, Wiesel A, Eldar YC, Hero AO (2010). “Shrinkage Algorithms for MMSE Covariance Estimation.” IEEE Transactions on Signal Processing, 58(10), 5016–5029. ISSN 1053-587X, 1941-0476.

Examples

## CRAN-purpose small computation
# set a seed for reproducibility
set.seed(11)

#  small data with identity covariance
pdim      <- 10
dat.small <- matrix(rnorm(5*pdim), ncol=pdim)

#  run the code
out.small <- CovEst.2010RBLW(dat.small)

#  visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
image(diag(pdim)[,pdim:1],     main="true cov")
image(cov(dat.small)[,pdim:1], main="sample cov")
image(out.small$S[,pdim:1],    main="estimated cov")
par(opar)

## Not run: 
## 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]*pdim), ncol=pdim) # 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
opar <- 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")
par(opar)

## End(Not run)