# Some code to look at the bootstrap sampling distribution of
# a correlation coefficient.

# Data (from Larson and Marx on ducks) (used by Owen, 1990)

x <- c(7,13,14,6,14,15,4,8,7,9,14)  # plumage index
y <- c(3,10,11,5,15,15,7,10,4,9,11) # behavioral index

# bootstrap sampling

boot.cor <- function(x,y,nboot)
	{
	 n <- length(x)
	 cboot <- numeric(nboot)
	 for( i in 1:nboot )
	  {
	   id <- sample( 1:n, size=n, replace=T )
	   cboot[i] <- cor( x[id], y[id] )  
	  }
	 return(cboot)
	}

test <- boot.cor(x,y,2000)

#postscript("boot.eps", horizontal=FALSE, onefile=FALSE, width=8, height=8)
par(mfrow=c(2,1), mgp=c(2,.5,0), mar=c(4,4,.1,.1) )

plot(x,y,xlab="plumage index", ylab="behavior index" )

alpha <- .05
ll <- sort(test)[2000*alpha/2]
rr <- sort(test)[2000*(1-alpha/2)]
#hist(test,nclass=100,xaxt="i",yaxt="i",xlab="sample correlation",ylab="bootstrap frequency" )
hist(test,nclass=100,xlab="sample correlation",ylab="bootstrap frequency" )
lines( rep(ll,2), c(-10,50), lwd=2)
lines( rep(rr,2), c(-10,50), lwd=2)
box(bty="l")
text(.15,130,"correlation coefficient = 0.825", adj=0,cex=.9)
text(.15,120,"95% bootstrap percentile interval (0.678,0.947)", adj=0,cex=.9)
#graphics.off()
#dev.off()