richard — Mar 3, 2014, 11:57 AM
## A simulation of Pearle's (1970) model for the EPR-Bohm correlations
set.seed(9875)
## For reproducibility. Replace integer seed by your own, or delete this line
## and let your computer dream up one for you (system time + process ID).
## Uniform random points on sphere generated using 'trig method' (method 3)
## of Dave Seaman: http://www.math.niu.edu/~rusin/known-math/96/sph.rand
## http://mathforum.org/kb/message.jspa?messageID=393612 See
## http://rpubs.com/gill1109/13340 for an R illustration
## (a) Choose z uniformly distributed in [-1,1]. (b) Choose t uniformly
## distributed on [0, 2*pi). (c) Let r = sqrt(1-z^2). (d) Let x = r *
## cos(t). (e) Let y = r * sin(t).
## Since the measurement directions are all in the equatorial plane, only z
## and x have been generated and treated as x and y
## Measurement angles for setting 'a': directions in the equatorial plane
angles <- seq(from = 0, to = 360, by = 1) * 2 * pi/360
K <- length(angles)
corrs <- numeric(K) ## Container for correlations
Ns <- numeric(K) ## Container for number of states
beta <- 0 * 2 * pi/360 ## Measurement direction 'b' fixed, in equatorial plane
M <- 10^6 ## Sample size. Next, try 10^7 or even 10^8 ...
## Use the same, single sample of 'M' realizations of hidden
## states for all measurement directions. This saves a lot of time,
## and reduces variance when we look at *differences*.
z <- runif(M, -1, 1)
t <- runif(M, 0, 2 * pi)
r <- sqrt(1 - z^2)
x <- r * cos(t)
e <- rbind(z, x) ## 2 x M matrix
## The M columns of e represent the x and y coordinates of M uniform random
## points on the sphere S^2
U <- runif(M)
s <- (2/sqrt(3*U+1)) - 1 # Pearle's "r" is arc cosine of "s"
b <- c(cos(beta), sin(beta)) ## Measurement vector 'b'
## Loop through measurement vectors 'a' (except last = 360 degrees = first)
for (i in 1:(K - 1)) {
alpha <- angles[i]
a <- c(cos(alpha), sin(alpha)) ## Measurement vector 'a'
ca <- colSums(e * a) ## Inner products of cols of 'e' with 'a'
cb <- colSums(e * b) ## Inner products of cols of 'e' with 'b'
good <- abs(ca) > s & abs(cb) > s ## Select the 'states'
N <- sum(good)
corrs[i] <- sum(sign(ca[good]) * sign(cb[good]))/N
Ns[i] <- N
}
corrs[K] <- corrs[1]
Ns[K] <- Ns[1]
plot(angles * 180/pi, corrs, type = "l", col = "blue", main = "Two correlation functions",
xlab = "Angle (degrees)", ylab = "Correlation")
points(angles * 180/pi, corrs, col = "blue", pch = ".", cex = 2)
lines(angles * 180/pi, cos(angles), col = "black")
points(angles * 180/pi, cos(angles), col = "black", pch = ".", cex = 2)
legend(x = 0, y = 0.2, legend = c("Pearle", "cosine"), text.col = c("blue",
"black"), lty = 1, col = c("blue", "black"))
plot(angles * 180/pi, corrs, type = "l", col = "blue", xlim = c(0, 90), ylim = c(0,
1), main = "Two correlation functions", xlab = "Angle (degrees)", ylab = "Correlation")
points(angles * 180/pi, corrs, col = "blue", pch = ".", cex = 2)
lines(angles * 180/pi, cos(angles), col = "black")
points(angles * 180/pi, cos(angles), col = "black", pch = ".", cex = 2)
legend(x = 0, y = 0.2, legend = c("Pearle", "cosine"), text.col = c("blue",
"black"), lty = 1, col = c("blue", "black"))
plot(angles * 180/pi, corrs, type = "l", col = "blue", xlim = c(0, 50), ylim = c(0.8,
1), main = "Two correlation functions", xlab = "Angle (degrees)", ylab = "Correlation")
points(angles * 180/pi, corrs, type = "b", col = "blue", pch = ".", cex = 2)
lines(angles * 180/pi, cos(angles), col = "black")
points(angles * 180/pi, cos(angles), type = "b", col = "black", pch = ".", cex = 2)
legend(x = 0, y = 0.85, legend = c("Pearle", "cosine"), text.col = c("blue",
"black"), lty = 1, col = c("blue", "black"))
plot(angles * 180/pi, corrs - cos(angles), type = "l", col = "blue",
main = "Difference, and +/- 1 standard error in red")
abline(h = 0, col = "black", lwd = 2)
lines(angles * 180/pi, 1/sqrt(Ns), col = "red")
lines(angles * 180/pi, -1/sqrt(Ns), col = "red")
max(abs(corrs - cos(angles)))
[1] 0.002333