### Generate Bivariate Normal Data

#######################################################################

# help(rmvnorm)
# rmvnorm(n, mean=rep(0,d), cov=diag(d), sd, rho, d=2)

n <- 30

rho <- 0.90

X <- rmvnorm(n, rho=0.90, d=2)

z1 <- X[,1]
z2 <- X[,2]

r.samp <- cor(z1,z2)		# rho.hat
r.samp

plot(z1,z2)

# Fisher's z-transformation

zrho <- c(-1,0,1)*1.96*(1/sqrt(n-3)) + 0.50*log((1+r.samp)/(1-r.samp))
zrho

# 95% Fisher's CI centered at the sample correlation, r.samp

(exp(2*zrho)-1)/(exp(2*zrho)+1)

# Nonparametric Bootstrap

B <- 1000

r.boot <- numeric(B)

for(i in 1:B){
	index <- seq(1:n)	
	index.boot <- sample(index, n, replace = T)
	r.boot[i] <- cor(z1[index.boot],z2[index.boot])
}

summary(r.boot)
r.boot.ci <- quantile(r.boot,c(0.025,0.975))
r.boot.ci
hist(r.boot)

r.boot.ztran <- 0.50*log((1+r.boot)/(1-r.boot))
hist(r.boot.ztran)

#######################################################################

# Nonparametric Bootstrap using bootstrap(), emp

# function to compute the correlation and Fisher's z-transformation

cor2 <- function(X){
	z1 <- X[,1]
	z2 <- X[,2]
	a <- cor(z1,z2)
	b <- 0.50*log((1+a)/(1-a))
	return(c(a,b))
}

r.boot <- bootstrap(X, cor2)

summary(r.boot)
plot(r.boot)

#######################################################################

# Nonparametric Bootstrap using bootstrap(), emp, bt, bca

# function to compute the correlation and Fisher's z-transformation

cor2 <- function(X){
	z1 <- X[,1]
	z2 <- X[,2]
	a <- cor(z1,z2)
	b <- 0.50*(log((1+a)/(1-a)) - log((1+r.samp)/(1-r.samp)))/sqrt(1/(n-3))
	return(c(a,b))
}

r.boot <- bootstrap(X, cor2)

plot(r.boot)

print(r.boot)

secor <- r.boot$estimate[1,3]	# location of the SE in the bootstrap data.frame
secor

r.boot.emp <- limits.emp(r.boot)
r.boot.emp

finv025 <- r.boot.emp[2,1]
finv025

finv975 <- r.boot.emp[2,4]
finv975

r.btci <- c(r.samp - (secor*finv95), r.samp - (secor*finv05))
r.btci

r.boot.bca <- limits.bca(r.boot, details = T)
r.boot.bca