# 2006 年 10 月 04 日作成
# Minimax risk and risk of improved estimators
sample=10
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
minimax
unbiased
repp=100000
ga=1
dime=2
risk1=rep(0,repp)
risk2=rep(0,repp)
risk3=matrix(0,ga,repp)
for (j in 1:repp){
  # dimention of complex normal distribution
  md2<-rep(0,dime)
  md3<-rep(0,dime)
  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)
  for (i in 1:dime){
    md2[i]=latent[i]/(sample +dime+1-2*i)
  }
  for (k in 1:ga){
    kk=k
    a1<-latent[1]^kk/(latent[1]^kk+latent[2]^kk)
    a2<-latent[2]^kk/(latent[1]^kk+latent[2]^kk)
    md3[1]=(a1/(sample +1 )+ a2/(sample-1))*latent[1]
    md3[2]=(a2/(sample +1 )+ a1/(sample-1))*latent[2]
    risk3[k,j]=sum(md3)-log(det(diag(md3)))-dime
  }
#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(md2)-log(det(diag(md2)))-dime
}
mrisk1=mean(risk1)
mrisk2=mean(risk2)
mrisk3=rep(0,ga)
for (k in 1:ga){
 mrisk3[k]=mean(risk3[k,]) 
}
reduction2=(mrisk1-mrisk2)/mrisk1
reduction3=(mrisk1-mrisk3)/mrisk1
for (k in 1:ga){
print(c(mrisk1,mrisk2,reduction2,round(k),mrisk3[k],reduction3[k]))
}