# 2009 年 03 月 14 日作成
# Simulation of risk for Large Covariance Estimation
# number of dimension
simtest<-function(x1,x2,x3){
p<-x1
ratio<-x2
id<-diag(rep(1,p))
################################################
# standard deviation of lognormal
logsd<-x3
################################################
# diagonal elements of parameter covarinace matrix
diag.elements<-rev(sort(rlnorm(p,meanlog=-logsd**2/2,meansd<-logsd)))
#diag.elements<-rep(1,p)
meaneig<-mean(diag.elements)
alpha<-(p-(sum(1/diag.elements))**2/sum(1/diag.elements**2))/p
diag.elements<-diag.elements/meaneig
mmeaneig<-mean(diag.elements)
malpha<-(p-(sum(1/diag.elements))**2/sum(1/diag.elements**2))/p
halfsig<-diag(sqrt(diag.elements))
insig<-diag(1/diag.elements)
################################################
# number of samples
n<-round(p*ratio)
################################################
# number of repitation
rep.no<-500
np<-n*p
srisk<-rep(0,rep.no)
sriskmle<-rep(0,rep.no)
haffrisk1<-rep(0,rep.no)
haffrisk2<-rep(0,rep.no)
###############################################
# exact risk of the MLE
eriskmle<-p*(p+1)/n
# exact risk of best multiple
erisk<-(p*p+p)/(p+n+1)
#################################################
# 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
mle<-wishart/n
bestmultiple<-wishart/(p+n+1)
quad<-(bestmultiple%*%insig-id)%*%(bestmultiple%*%insig-id)
quadmle<-(mle%*%insig-id)%*%(mle%*%insig-id)
for (j in 1:p){
srisk[i]<-srisk[i]+quad[j,j]
}
for (j in 1:p){
sriskmle[i]<-sriskmle[i]+quadmle[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
be<-mean(srisk); mle<-mean(sriskmle); hf1<-mean(haffrisk1); hf2<-mean(haffrisk2)
#print("mean of eigenvalues"); print(meaneig)
#print("modified mean of eigenvalues"); print(mmeaneig)
print("value of alpha"); print(alpha)
#print("modified value of alpha"); print(malpha)
print("exact risk of mle"); print(eriskmle)
print("simulated risk  of mle"); print(mle)
print("exact risk of the best multiplier"); print(erisk)
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 over mle"); print(100*(mle-hf2)/mle)
print("risk reduction over the best multiplier"); print(100*(be-hf2)/be)
}