library(bootstrap)
set.seed(10102)

postscript("Boot-varstab-data.ps")
plot(law,pch=19,xlab="LSAT",ylab="GPA")
dev.off()

law.mat<-as.matrix(law)
rho.hat<-cor(law.mat)[1,2]

n<-nrow(law.mat)

index<-c(1:15)
B1<-10000
B2<-500
r.boot.stab<-se.boot.stab<-rep(0,B1)
r.boot.stab.inner<-rep(0,B2)
for(b1 in 1:B1){

	ind.boot.stab<-sample(index,size=n,rep=T)
	law.boot.stab<-law.mat[ind.boot.stab,]
	r.boot.stab[b1]<-cor(law.boot.stab)[1,2]

	for(b2 in 1:B2){
		inner.ind<-sample(index,size=n,rep=T)
		boot.inner<-law.boot.stab[inner.ind,]
		r.boot.stab.inner[b2]<-cor(boot.inner)[1,2]
	}
	se.boot.stab[b1]<-sqrt(sum((r.boot.stab.inner-mean(r.boot.stab.inner))^2)/B2)

	print(c(b1,r.boot.stab[b1],se.boot.stab[b1]))

}

dx<-0.0001
rx<-seq(0.2,0.9999,by=dx)
ry<-predict(loess(se.boot.stab~r.boot.stab),rx)

postscript("Boot-varstab-s.ps")
plot(r.boot.stab,se.boot.stab,pch=19,xlab=expression(hat(theta)^{(b)}),ylab=expression(hat(se)^{(b)}))
lines(rx,ry,col="red",lwd=2)
dev.off()

sx<-rx
sy<-ry

Sy<-cumsum(1/sy)*dx

postscript("Boot-varstab-ghat.ps")
plot(sx,Sy,type="l",ylab=expression(g(hat(theta))),xlab=expression(hat(theta)),col="red",lwd=2)
lines(sx,0.5*log((1+sx)/(1-sx)),lty=2,lwd=2)
dev.off()

tx<-sx
gx<-Sy

B3<-2000
r.boot.stab3<-se.boot.stab3<-phi.boot<-rep(0,B3)

for(b3 in 1:B3){

	ind.boot.stab3<-sample(index,size=n,rep=T)
	law.boot.stab3<-law.mat[ind.boot.stab3,]
	r.boot.stab3[b3]<-cor(law.boot.stab3)[1,2]

	i1<-which.min(abs(tx-r.boot.stab3[b3]))

	phi.boot[b3]<-Sy[i1]
}

boot.ord<-order(phi.boot)


postscript("Boot-varstab-phi.ps")
par(mfrow=c(1,2))
hist(r.boot.stab3,br=seq(0.0,1,by=0.05),main=" ",xlab=expression(hat(rho)),col="black",border="white",ylim=range(0,400))
hist(phi.boot,br=c(0:20)/2,main=" ",xlab=expression(hat(phi)),col="black",border="white",ylim=range(0,400))
dev.off()


Z.int<-0.5*log((1+rho.hat)/(1-rho.hat))+qnorm(c(0.05,0.95))/sqrt(n-3)
Fisher.int<-(exp(2*Z.int)-1)/(exp(2*Z.int)+1)
g.int<-r.boot.stab3[boot.ord][c(100,1900)]
boot.int<-quantile(r.boot.stab,prob=c(0.05,0.95))