# 2009 年 03 月 19 日作成
# 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 matrix1
logsd<-rep(x3,p)
diag.elements<-ifelse(logsd>0,rlnorm(p,meanlog=-logsd[1]**2/2,meansd<-logsd[1]),rep(1,p))
diag.elements<-rev(sort(diag.elements))
#diag.elements<-rev(sort(rlnorm(p,meanlog=-logsd[1]**2/2,meansd<-logsd[1])))
meaneig<-mean(diag.elements)
diag.elements<-diag.elements/meaneig
diag.elements<-ifelse(diag.elements>1.0e-1,diag.elements,1.0e-1)
alpha<-(p-(sum(1/diag.elements))**2/sum(1/diag.elements**2))/p
#mmeaneig<-mean(diag.elements)
#malpha<-(p-(sum(1/diag.elements))**2/sum(1/diag.elements**2))/p
#alphabeta<-(p+1)/(n*malpha+p+1)
#alphabeta1<-1-(n+p+1)*malpha/(n*malpha+p+1)
halfsig<-diag(sqrt(diag.elements))
insig<-diag(1/diag.elements)
################################################
# number of samples
n<-round(p*ratio)
################################################
# number of repitation
rep.no<-1000
np<-n*p
srisk<-rep(0,rep.no)
sriskmle<-rep(0,rep.no)
haffrisk1<-rep(0,rep.no)
haffrisk2<-rep(0,rep.no)
haffrisk12<-rep(0,rep.no)
haffrisk22<-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]
}
cons2<-2*(n-1)*(p-n-1)/((p-n+1)*(p-n+3))
cons<-(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)
haff12<-(wishart+(cons2/invt)*semiorth%*%t(semiorth))/(p+n+1)
quadh1<-(haff1%*%insig-id)%*%(haff1%*%insig-id)
quadh12<-(haff12%*%insig-id)%*%(haff12%*%insig-id)
for (j in 1:p){
haffrisk1[i]<-haffrisk1[i]+quadh1[j,j]
haffrisk12[i]<-haffrisk12[i]+quadh12[j,j]
}
# calculation of estimate for the nonsingular Haff
haff2<-(wishart+(cons/invt)*id)/(p+n+1)
haff22<-(wishart+(cons2/invt)*id)/(p+n+1)
quadh2<-(haff2%*%insig-id)%*%(haff2%*%insig-id)
quadh22<-(haff22%*%insig-id)%*%(haff22%*%insig-id)
for (j in 1:p){
haffrisk2[i]<-haffrisk2[i]+quadh2[j,j]
haffrisk22[i]<-haffrisk22[i]+quadh22[j,j]
}
}
###########################################################
# Output of results
be<-mean(srisk); mle<-mean(sriskmle)
besd<-sd(srisk)/sqrt(rep.no); mlesd<-sd(sriskmle)/sqrt(rep.no)
hf1<-mean(haffrisk1); hf2<-mean(haffrisk2)
hf12<-mean(haffrisk12); hf22<-mean(haffrisk22)
hf1sd<-sd(haffrisk1)/sqrt(rep.no); hf2sd<-sd(haffrisk2)/sqrt(rep.no)
hf12sd<-sd(haffrisk12)/sqrt(rep.no); hf22sd<-sd(haffrisk22)/sqrt(rep.no)
print("************** settings of simulation **************************")
print("dimension"); print(p)
print("sample size"); print(n)
print("standard deviation of lognormal");print(logsd[1]) 
#print("mean of eigenvalues"); print(meaneig)
#print("modified mean of eigenvalues"); print(mmeaneig)
print("value of alpha"); print(alpha)
#print("value of alphabeta"); print(100*alphabeta)
#print("value of alphabeta1"); print(100*alphabeta1)
print("maximum of eigenvalues"); print(diag.elements[1])
print("minimun of eigenvalues"); print(diag.elements[p])
print("mean of eigenvalues"); print(mean(diag.elements))
print("standard deviation of eigenvalues"); print(sd(diag.elements))
#print("modified value of alpha"); print(malpha)
print("************** results of simulation **************************")
print("exact risk of mle"); print(eriskmle)
print("simulated risk  of mle"); print(mle)
print("standard error(simulated  risk  of mle)"); print(mlesd)
print("exact risk of the best multiplier"); print(erisk)
print("simulated risk  of the best multiple"); print(be)
print("standard error"); print(besd)
print("best multiplier's risk reduction over mle"); print(100*(mle-be)/mle)
print("####### singular #########")
print("simulated risk  of the singular Haff estimator"); print(hf1)
print("standard error"); print(hf1sd)
print("singular Haff's risk reduction over mle"); print(100*(mle-hf1)/mle)
print("singular Haff's risk reduction over the best multiplier"); print(100*(be-hf1)/be)
print("### singular2 ###")
print("simulated risk  of the singular 2Haff estimator"); print(hf12)
print("standard error"); print(hf12sd)
print("singular 2Haff's risk reduction over mle"); print(100*(mle-hf12)/mle)
print("singular 2Haff's risk reduction over the best multiplier"); print(100*(be-hf12)/be)
print("######## nonsingular ########")
print("simulated risk  of the nonsingular Haff estimator"); print(hf2)
print("standard error"); print(hf2sd)
print("nonsingular Haff's risk reduction over mle"); print(100*(mle-hf2)/mle)
print("nonsingular Haff's risk reduction over the best multiplier"); print(100*(be-hf2)/be)
print("### nonsingular2 ####")
print("simulated risk  of the nonsingular 2Haff estimator"); print(hf22)
print("standard error"); print(hf22sd)
print("nonsingular 2Haff's risk reduction over mle"); print(100*(mle-hf22)/mle)
print("nonsingular 2Haff's risk reduction over the best multiplier"); print(100*(be-hf22)/be)
print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print("%%%%%%%%%%%%%%%%%%%%% End of simulation     %%%%%%%%%%%%%%%%%%%%%%%%%%%")
print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
return(c(n,p,alpha,hf12,100*(be-hf12)/be,hf22,100*(be-hf22)/be))
}
######################################################################
######################################################################
reppp<-100
a<-matrix(0,1,7)
for (i in 1:reppp) {
a<-append(a,simtest(20,4/5,3))
}

print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print("lnorm=3")
print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
ab<-matrix(a,,7,byrow=TRUE)
reppp1<-reppp+1
ab1<-ab[2:reppp1,]
print(ab1)