## Population correlation coefficients to consider.
R = seq(-0.9, 0.9, 0.1)

## Sample size.
n = 20

## Number of simulation replications.
nrep = 1e3

PROP=NULL

## Loop over the different population correlation coefficients.
for (j in 1:length(R)) 
{
  r = R[j]

  ## Generate nrep sample correlation coefficients.
  C = array(0, nrep)
  for (i in 1:nrep) 
  {
    X = rnorm(n)
    Y = r*X + sqrt(1-r^2)*rnorm(n)
    C[i] = cor(X,Y)
  }

  ## Genearte nrep approximate 95% CI's

  logterm = 0.5*log(1+C)-0.5*log(1-C) 
  CL = logterm-2/sqrt(n-3)
  CU = logterm+2/sqrt(n-3)
  
  LB = (exp(2*CL)-1)/(exp(2*CL)+1)
  UB = (exp(2*CU)-1)/(exp(2*CU)+1)

  ##  Find the proportion of population pearson-correlation coeffs
  ##  that are contained in the 95 percent CI for this value value 
  ##  of r

  PROP = c(PROP, mean( (LB <= r) & (UB >= r) ))
}

## Look at the output
print(cbind(R, PROP))






## Population correlation coefficients to consider.
R = seq(-0.9, 0.9, 0.1)

## Sample size.
n = 20

## Number of simulation replications.
nrep = 1e3

POWER=NULL

## Loop over the different population correlation coefficients.
for (j in 1:length(R)) 
{
  r = R[j]

  ## Generate nrep sample correlation coefficients.
  C = array(0, nrep)
  for (i in 1:nrep) 
  {
    X = rnorm(n)
    Y = r*X + sqrt(1-r^2)*rnorm(n)
    C[i] = cor(X,Y)
  }

  ## Genearte nrep approximate 95% CI's

  logterm = 0.5*log(1+C)-0.5*log(1-C) 
  CL = logterm-2/sqrt(n-3)
  CU = logterm+2/sqrt(n-3)
  
  LB = (exp(2*CL)-1)/(exp(2*CL)+1)
  UB = (exp(2*CU)-1)/(exp(2*CU)+1)

  ##  Find the proportion of CI's
  ##  that do not contain 0

  POWER = c(POWER, mean( (LB > 0) | (UB < 0) ))
}

## Look at the output
print(cbind(R, POWER))




dichot=function(X,Y,c)
{
  ## Initialize the table
  N=array(0,c(2,2))
  
  ## Dichotomize the data
  A=1*(X>=c)
  B=1*(Y>=c)
  
  ##  Populate the table
  N[1,1]=sum( (A==1) & (B==1) )
  N[1,2]=sum( (A==1) & (B==0) )
  N[2,1]=sum( (A==0) & (B==1) )
  N[2,2]=sum( (A==0) & (B==0) )
  return(N)
}


n=100
r=0.01
## Generate nrep sample correlation coefficients.
C = array(0, nrep)
SLOR = array(0, nrep)
SE = array(0, nrep)
for (i in 1:nrep) 
{
  X = rnorm(n)
  Y = r*X + sqrt(1-r^2)*rnorm(n)
  C[i] = cor(X,Y)
  N=dichot(X,Y,c=0)
  SLOR[i] = log(N[1,1]) + log(N[2,2]) - log(N[1,2]) - log(N[2,1])
  SE[i] = sqrt(1/N[1,1] + 1/N[1,2] + 1/N[2,1] + 1/N[2,2])
}





cor.406 = function(X,Y)
{
  ## X,Y are n by reps matrices
  DY=scale(Y, center=TRUE, scale=FALSE)
  DX=scale(X, center=TRUE, scale=FALSE)  
  r=apply(DX*DY,2,sum)/sqrt(apply(DX^2,2,sum)*apply(DY^2,2,sum))
  return(r)
}