# Plots the sampling distribution of r for specifiec value of rho

rho <- scan(); rho
N = 50
x <- rnorm(10000)    
y <- rnorm(10000)
zx <- (x - mean(x))/sd(x)
zy <- (y-mean(y))/sd(y)
a <- rho/sqrt(1-rho^2)
z <- a*zx + zy
cor(x,z)                                                          
# This will produce variables x and y with a correaltion very close to rho.



# Now draw 10,000 samples (shuffled) and draw histogram
nreps <- 1000
rsamp <- numeric(nreps)
for (i in 1:nreps) {
  xnew <- sample(x,N, replace = T)
  ynew <- sample(y,N, replace = TRUE)
  ### I need to look up how to sample xy pairs from a data frame or matrix
  ###  See an earlier program
  
      #ynew <- y  #Not necessary
  rsamp[i] <- cor(xnew, ynew)
  }
a <- sort(rsamp)
lower <- a[.025*nreps]
upper <- a[.975*nreps]
cat("Using resampling the lower and upper confidence limits are ","\n")
limits <- list(c("lower" = lower, "upper" = upper))
limits
hist(rsamp, breaks = 50)