# 2006 年 01 月 10 日作成
# Simulated risk of shrinkage estimators in possion estimation problem
# p; dimension of Possion models
# rep: reputation of simulation

riskpo<-function(p,lambda)
{
rep<-5000
risk1<-rep(0,rep)
risk2<-rep(0,rep)
risk3<-rep(0,rep)
x<-rep(0,p)
# simulation process starts here
for (j in 1:rep){
for (i in 1:p){
x[i]<-rpois(1,lambda[i])
}
shrink1<-(1-p/(sum(x)+p))*x
shrink1p<-(shrink1+abs(shrink1))/2
shrink2<-(1-2*(p-1)/(sum(x)+p))*x
if(sum(x)<=p-1) shrink2<-x
shrink2p<-(shrink2+abs(shrink2))/2
#print(shrink1p)
#print(x)
risk1[j]<-sum((x-lambda)**2/lambda)
risk2[j]<-sum((shrink1p-lambda)**2/lambda)
risk3[j]<-sum((shrink2p-lambda)**2/lambda)
}
# Calculate simulated risk here
mrisk1<-mean(risk1)
sdrisk1<-sd(risk1)/sqrt(rep)
mrisk2<-mean(risk2)
sdrisk2<-sd(risk2)/sqrt(rep)
mrisk3<-mean(risk3)
sdrisk3<-sd(risk3)/sqrt(rep)
print("risk of mle ")
print(mrisk1)
print(sdrisk1)
print("risk of shrinkage 1 ")
print(mrisk2)
print(sdrisk2)
print("improvement of shrinkage 1 ")
print((mrisk1-mrisk2)/mrisk1)
print("risk of shrinkage 3 ")
print(mrisk3)
print(sdrisk3)
print("improvement of shrinkage 2 ")
print((mrisk1-mrisk3)/mrisk1)
}