theta<-1
X<-matrix(rpois(250000,lambda=theta),1000,250)
Y<-X
Y[X==0]<-1
Y[X>0]<-0
# 
biasmle<-mean(exp(-X[,1]))-exp(-theta)
biasmvue<-mean(Y[,1])-exp(-theta)
print(c("biasmle=", biasmle))
print(c("biasmvue=",biasmvue))
# try for different values of theta
#
#
umvue<-function(datasum,n) {
  umvue<-(1-1/n)^datasum }
#
XX<-t(apply(X,1,cumsum))
# applies sum operator cumulatively to each of the rows so that
# XX[2,10] is the sum of the 10th first element in row 2 etc.
#
XXm<-XX
for (kk in (1:dim(X)[1])) {
 XXm[kk,]<-XX[kk,]/seq(1,dim(X)[2]) }
# creates the mean vector for X
YYm<-XXm
for (jj in (1:dim(X)[2])) {
 YYm[,jj]<-(1-1/jj)^(jj*XXm[,jj]) }
# creates the UMVUEs
#
mle<-exp(-XXm)
mvue<-YYm
plot(seq(1,dim(X)[2]),apply(mle,2,mean)-exp(-theta),type='l',xlab="n", ylab="bias")
lines(seq(1,dim(X)[2]),apply(mvue,2,mean)-exp(-theta),lty=2)
title(main="BIAS, mle solid line, mvue dashed line")
abline(h=0,lty=3)
# run the above a couple of times to see that on average bias cancels out
# for umvue but not for mle.
#
plot(seq(1,dim(X)[2]),apply(mle,2,var),type='l',xlab="n", ylab="bias")
lines(seq(1,dim(X)[2]),apply(mvue,2,var),lty=2)
title(main="VARIANCE, mle solid line, mvue dashed line")
# variance not much difference
#
#
# bias correction 
# XX is the data matrix
biascorr<-function(dataset,n) {
  bhat<-(n-1)*(exp(-(n/(n-1))*mean(dataset))/n*sum(exp(dataset/(n-1)))-exp(-mean(dataset)))
  biascorr<-exp(-mean(dataset))-bhat }
# note, biascorrection for n=1 is 0.
#
biascorrmle<-mle
for (jj in (2:dim(X)[2])) {
  biascorrmle[,jj]<-apply(X[,1:jj],1,biascorr,n=jj) }
#
plot(seq(1,dim(X)[2]),apply(mle,2,mean)-exp(-theta),type='l',xlab="n", ylab="bias")
lines(seq(1,dim(X)[2]),apply(mvue,2,mean)-exp(-theta),lty=2)
lines(seq(1,dim(X)[2]),apply(biascorrmle,2,mean)-exp(-theta),lty=5,col=3)
title(main="BIAS, mle solid line, mvue dashed line, biascorr mle green dotted")
abline(h=0,lty=3)
# bias essentially gone!
plot(seq(1,dim(X)[2]),apply(mle,2,var),type='l',xlab="n", ylab="bias")
lines(seq(1,dim(X)[2]),apply(mvue,2,var),lty=2)
lines(seq(1,dim(X)[2]),apply(biascorrmle,2,var),lty=5,col=3)
title(main="VARIANCE, mle solid line, mvue dashed line, biascorr mle green dotted")
#
#
#  bootstrap
# X contain true data sets
# I will only do this for nvec
nvec<-c(2,5,10,15,20,25,30,40,50,75,100,250)
B<-250
# 250 bootstrap samples
mlebootparam<-mle
mlebootnonparam<-mle
for (kk in (1:length(nvec))) {
  #different n I will examine
 for (jj in (1:dim(X)[1])) {
   # the number of datarepetions I do this for to examine average bias etc
  paramest<-rep(0,B)
  nparamest<-rep(0,B)
  #here I will store the bootstrap estimates
  for (bb in (1:B)) {
      xnewnp<-sample(X[jj,1:nvec[kk]],size=nvec[kk],replace=T)
      xnewp<-rpois(nvec[kk],lambda=mean(X[jj,1:nvec[kk]]))
      #
      paramest[bb]<-exp(-mean(xnewp))
      nparamest[bb]<-exp(-mean(xnewnp)) }
  biaspar<-mean(paramest)-exp(-mean(X[jj,1:nvec[kk]]))
  biasnpar<-mean(nparamest)-exp(-mean(X[jj,1:nvec[kk]]))
  #estimating bias using the mean of the boostrap estimates
  mlebootparam[jj,nvec[kk]]<-mle[jj,nvec[kk]]-biaspar
  mlebootnonparam[jj,nvec[kk]]<-mle[jj,nvec[kk]]-biasnpar
  #biascorrection
 }}
#
plot(seq(1,dim(X)[2]),apply(mle,2,mean)-exp(-theta),type='l',xlab="n", ylab="bias")
lines(seq(1,dim(X)[2]),apply(mvue,2,mean)-exp(-theta),lty=2)
lines(seq(1,dim(X)[2]),apply(biascorrmle,2,mean)-exp(-theta),lty=5,col=3)
lines(nvec,apply(mlebootparam[,nvec],2,mean)-exp(-theta),lty=4,col=2)
lines(nvec,apply(mlebootnonparam[,nvec],2,mean)-exp(-theta),lty=1,col=5)
title(main="BIAS, mle solid line, mvue dashed line, biascorr mle green dotted")
abline(h=0,lty=3)