# 2006 年 02 月 28 日作成
# number of repitation of simulation
repp=100
risk1=rep(0,repp)
risk2=rep(0,repp)
for (j in 1:repp){
sample=50
# dimention of complex normal distribution
dime=40
md<-rep(0,dime)
ecov=matrix(rep(0,dime**2),dime,dime)
for (i in 1:sample){
z=complex(real=rnorm(dime,0,1),imag=rnorm(dime,0,1))
ecov=ecov+Conj(z)%*%t(z)
}
eigen.w=eigen(ecov)
latent=c(eigen.w$value)
hmat<-matrix(eigen.w$vectors,dime,dime)
for (i in 1:dime){
md[i]=latent[i]/(sample +dime+1-2*i)
}
#ecov/sample
#hmat%*%diag(latent)%*%t(Conj(hmat))/sample
mecov=hmat%*%diag(md)%*%t(Conj(hmat))
risk1[j]=sum(latent)/sample-log(det(diag(latent)/sample))-dime
risk2[j]=sum(md)-log(det(diag(md)))-dime
}
mrisk1=mean(risk1)
mrisk2=mean(risk2)
reduction=(mrisk1-mrisk2)/mrisk1
print(c(mrisk1,mrisk2,reduction))