# 2008 年 03 月 15 日作成
# Simulation of risk for Large Covariance Estimation
# number of dimension
p<-100
id<-diag(rep(1,p))
diag.elements<-c(rep(5,p/2),rep(1,p/2))
halfsig<-diag(sqrt(diag.elements))
insig<-diag(1/diag.elements)
# number of samples
n<-50
# number of repitation
rep.no<-100
np<-n*p
srisk<-rep(0,rep.no)
haffrisk1<-rep(0,rep.no)
haffrisk2<-rep(0,rep.no)
# exact risk of best multiple
erisk<-(p*p+p)/(p+n+1)
print("exact risk"); print(erisk)
# Simulation of risk
for (i in 1:rep.no){
xx<-rnorm(n*p,0,1)
x<-matrix(xx,nrow=p,ncol=n)
w<-x%*%t(x)
# generate Wishart matrix
wishart<-halfsig%*%w%*%halfsig
# calculation of estimate for the best multiple
bestmultiple<-wishart/(p+n+1)
quad<-(bestmultiple%*%insig-id)%*%(bestmultiple%*%insig-id)
for (j in 1:p){
srisk[i]<-srisk[i]+quad[j,j]
}
cons<-2*(n-1)*(p-n-1)/((p-n+1)*(p-n+3))
eigen.w<-eigen(wishart)$values
orth<-eigen(wishart)$vectors
semiorth<-orth[,1:n,drop=FALSE]
invt<-sum(1/eigen.w[1:n])
# calculation of estimate for the singular Haff
haff1<-(wishart+(cons/invt)*semiorth%*%t(semiorth))/(p+n+1)
quadh1<-(haff1%*%insig-id)%*%(haff1%*%insig-id)
for (j in 1:p){
haffrisk1[i]<-haffrisk1[i]+quadh1[j,j]
}
# calculation of estimate for the nonsingular Haff
haff2<-(wishart+(cons/invt)*id)/(p+n+1)
quadh2<-(haff2%*%insig-id)%*%(haff2%*%insig-id)
for (j in 1:p){
haffrisk2[i]<-haffrisk2[i]+quadh2[j,j]
}
}
# Output of results
print("exact risk of the best multiple"); print(erisk)
be<-mean(srisk); hf1<-mean(haffrisk1); hf2<-mean(haffrisk2)
print("simulated risk  of the best multiple"); print(be)
print("simulated risk  of the singular Haff estimator"); print(hf1)
print("risk reduction"); print(100*(be-hf1)/be)
print("simulated risk  of the nonsingular Haff estimator"); print(hf2)
print("risk reduction"); print(100*(be-hf2)/be)