library(bootstrap)

y<-patch$y
z<-patch$z
n<-length(y)

plug.in<-mean(y)/mean(z)

Bvec<-c(seq(50,950,by=50),1000,5000,10000,20000,100000)
bias.est<-matrix(0,nrow=length(Bvec),ncol=2)
ind<-1:n

for(ib in 1:length(Bvec)){

	B<-Bvec[ib]
	boot.rat<-rep(0,B)
	Pj.tot<-0
	for(b in 1:B){
		iboot<-sort(sample(ind,size=n,rep=T))
		Pj.tab<-table(iboot)
		Pj<-rep(0,n)
		Pj[as.numeric(names(Pj.tab))]<-Pj.tab
		Pj.tot<-Pj.tot+Pj
		yboot<-y[iboot]
		zboot<-z[iboot]
		boot.rat[b]<-mean(yboot)/mean(zboot)
	}
	Pj.mean<-Pj.tot/B

	hist(boot.rat)

	bias.est[ib,1]<-mean(boot.rat)-plug.in
	bias.est[ib,2]<-mean(boot.rat)-sum(Pj.mean*y)/sum(Pj.mean*z)

}
bs<-20
postscript("Boot-bias.ps")
plot(Bvec[1:bs],bias.est[1:bs,1],pch=19,xlab="B",ylab="Est. Bias")
lines(Bvec[1:bs],bias.est[1:bs,1])
points(Bvec[1:bs],bias.est[1:bs,2],col="red",pch=19)
lines(Bvec[1:bs],bias.est[1:bs,2],col="red",lty=2)
abline(h=bias.est[nrow(bias.est),1])
abline(h=bias.est[nrow(bias.est),2],col="red",lty=2)
legend(600,0.03,c(expression(hat(b))," ",expression(tilde(b))),lty=c(1,0,2),col=c("black","white","red"),lwd=2)
dev.off()

postscript("Boot-bias-Hormone.ps")
hist(boot.rat,col="black",border="white",xlab=expression(hat(theta)^"*"),main=" ")
abline(v=c(-0.2,0.2),col="red",lwd=2)
dev.off()