# R-code for Fisher z-transformation 

# generate bivariate normals, mean 0, variances 1, correlation rho

n<-10000
R<-c(1:n) 
Z<-c(1:n) 
for(i in 1:n){ 
X1<- rnorm(10) 
X2<- rnorm(10) 
rho <- 0.6
Y1<- X1*( sqrt((1+rho)/2))  + X2*( sqrt((1-rho)/2))  
Y2<- X1*( sqrt((1+rho)/2))  - X2*( sqrt((1-rho)/2))   
R[i]<- cor(Y1,Y2)
Z[i]<- 1/2*log( ( 1+ R[i])/(1-R[i])) 
 } 
par(mfrow=c(1,2))
hist(R, xlab="histogram for R; rho=0.6;  n=10; m=10,000 repeats") 
hist(Z, xlab="histogram for Z, n=10; m=10,000 repeats") 

# the histogram shows that the distribution of Z is closer to normal than the 
# distribution of R 


# calculate the approximate 95% confidence limits for one value: 

confint1lower <- R[100]-1.96*(1-R[100] )/ sqrt(10)
confint1upper<-  R[100]+1.96*(1-R[100] )/ sqrt(10)
confint1lower 
confint1upper

c2l<- Z[100] - 1.96/ sqrt(10)
c2u<- Z[100] + 1.96/ sqrt(10)

confint2lower<- (exp(2*c2l)-1)/(exp(2*c2l) + 1) 
confint2upper<- (exp(2*c2u)-1)/(exp(2*c2u) + 1) 

confint2lower
confint2upper 

# the confidence interval using the Fisher transform tends to be wider

# now see how often the true parameter is covered by the confidence interval: 


SR<-c(1:n) # the vector: 1 if rho is in the interval, 0 if not, for R 
SZ<-c(1:n) # the vector: 1 if rho is in the interval, 0 if not, for Z 

confint1lower<-c(1:n) 
confint1upper<-c(1:n)
confint2lower<-c(1:n)
confint2upper<-c(1:n)

for(i in 1:n){ 

confint1lower[i] <- R[i]-1.96*(1-R[i] )/ sqrt(10)
confint1upper[i]<-  R[i]+1.96*(1-R[i] )/ sqrt(10)

SR[i]<-1
if(confint1upper[i] <= rho) SR[i]<-0
if(rho<=confint1lower[i]) SR[i]<-0 

c2l<- Z[i] - 1.96/ sqrt(10)
c2u<- Z[i] + 1.96/ sqrt(10)

confint2lower[i]<- (exp(2*c2l)-1)/(exp(2*c2l) + 1) 
confint2upper[i]<- (exp(2*c2u)-1)/(exp(2*c2u) + 1) 

SZ[i]<-1
if(confint2upper[i] <= rho) SZ[i]<-0
if(rho<=confint2lower[i]) SZ[i]<-0 
}

sum(SR)
sum(SZ)

# we see that the true parameter is contained much more often in the confidence
# interval using the Fisher z-transformation.