sample=20
repp=10000
dime=2
risk1=rep(0,repp)
for (j in 1:repp){
  # dimention of complex normal distribution
  ecov=matrix(rep(0,dime**2),dime,dime)
  for (i in 1:sample){
    z=rnorm(dime,0,1)
    ecov=ecov+z%*%t(z)
  }
  eigen.w=eigen(ecov)
  latent=c(eigen.w$value)
  #hmat<-matrix(eigen.w$vectors,dime,dime)
  risk1[j]=sum(latent)/sample-log(latent[1]/sample)-log(latent[2]/sample)-dime
}
mrisk1=mean(risk1)
exga<-function(x){log(x)*dgamma(x,sample)}
m1<-integrate(exga,0.0001,20000)$value
exga<-function(x){log(x)*dgamma(x,sample-1)}
m2<-integrate(exga,0.0001,20000)$value
minimax<-log(sample+1)+log(sample-1)-m1-m2
unbiased<-2*log(sample)-m1-m2
print(c(unbiased,mrisk1,minimax))