richard — Feb 17, 2014, 8:21 PM
## Simultaneous simulation of the "Christian 2.0" model for the EPR-B correlations,
## in two versions.
## "Christian 2.0" refers to the second edition of Joy Christian's book, see
## http://libertesphilosophica.info/blog/
## S^1 version (2-D): implementation of Michel Fodje's algorithm:
## https://github.com/minkwe/epr-simple/
## S^2 version (3-D): same, but with uniform random points on S^2 (subset R^3)
## instead of uniform random points on S^1 (subset R^2).
## This is now identical to an extension of Caroline Thompson's
## "chaotic rotating ball model" - instead of circular caps (on the sphere)
## of a fixed radius, we now have circular caps with a common, random radius.
## http://arxiv.org/abs/quant-ph/0210150
## http://freespace.virgin.net/ch.thompson1/
## One can also think of this as a kind of "fuzzy" membership of the circular caps.
## For another kind of fuzzy membership, see Chantal's model
## http://rpubs.com/gill1109/13329 and
## http://rpubs.com/chenopodium (various versions EPR experiment)
## I only look at measurement directions in the x-y plane and therefore
## don't bother with generating the z component of the points on S^2
## (irrelevant to inner product of state vector with measurement vector).
set.seed(5678)
## For reproducibility. Replace integer seed by
## your own, or delete this line and let your computer
## dream up one for you in its own mysterious way (system time?).
## 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
## (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).
## I simultaneously generate "M" uniform points on the circle using t as above,
## then defining x and y as in (d) and (e) with r = 1
## The measurement angles for setting "a"
angles <- seq(from = 0, to = 360, by = 7.5) * 2 * pi / 360
K <- length(angles)
corrs2 <- numeric(K) ## Container for correlations for S^1 model ("2" refers to R^2)
Ns2 <- numeric(K) ## Container for number of states for S^1 model
corrs3 <- numeric(K) ## Container for correlations for S^2 model ("3" refers to R^3)
Ns3 <- numeric(K) ## Container for number of states for S^2 model
beta <- 0 * 2 * pi / 360 ## Measurement direction "b" is fixed
M <- 10^6 ## Sample size.
## I use the same, single sample of "M" realizations of hidden states
## for all measurement directions and for both models.
## This saves a lot of time, and reduces variance when we look at *differences*
t <- runif(M, 0, 2*pi)
x2 <- cos(t)
y2 <- sin(t)
e2 <- rbind(y2, x2) ## 2 x M matrix
## The M columns of e2 represent M uniform random points on the circle S^1
z <- runif(M, -1, 1)
r <- sqrt(1 - z^2)
x3 <- r * x2
## y3 <- r * y2 # (we don't need the third coordinate here)
## My 3-vector (point on S^2) is Seaman's (z, x, y).
## I discard the third component since I only look at
## measurement vectors "a", "b" which lie in the plane
## generated by the first two components
e3 <- rbind(z, x3) ## 2 x M matrix
## e3alt <- rbind(z, x3, y3) #(this would give us the full 3 x M matrix)
## The M columns of e3 represent the orthogonal projections
## of M uniform random points on the sphere S^2,
## projected onto the horizontal plane
theta <- runif(M, 0, pi/2) ## Joy and Minkwe's "theta_0"
s <- (sin(theta)^2) / 2
## The measurement vectors all live in the plane - in fact, on the circle S^1
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"
ca2 <- colSums(e2*a) ## Inner products of cols of "e2" with "a"
cb2 <- colSums(e2*b) ## Inner products of cols of "e2" with "b"
ca3 <- colSums(e3*a) ## Inner products of cols of "e3" with "a"
cb3 <- colSums(e3*b) ## Inner products of cols of "e3" with "b"
good2 <- abs(ca2) > s & abs(cb2) > s ## Select the "states" for S^1 model
good3 <- abs(ca3) > s & abs(cb3) > s ## Select the "states" for S^2 model
N2 <- sum(good2)
N3 <- sum(good3)
corrs2[i] <- sum(sign(ca2[good2])*sign(cb2[good2]))/N2
Ns2[i] <- N2
corrs3[i] <- sum(sign(ca3[good3])*sign(cb3[good3]))/N3
Ns3[i] <- N3
}
corrs2[K] <- corrs2[1]
Ns2[K] <- Ns2[1]
corrs3[K] <- corrs3[1]
Ns3[K] <- Ns3[1]
plot(angles*180/pi, corrs2, type="l", col="blue", pch = ".", cex = 2,
main = "Three correlation functions",
xlab = "Angle (degrees)",
ylab = "Correlation")
lines(angles*180/pi, corrs3, col="red", pch = ".", cex = 2)
lines(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2)
points(angles*180/pi, corrs3, col="red", pch = ".", cex = 2)
points(angles*180/pi, corrs2, col="blue", pch = ".", cex = 2)
points(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2)
legend(x=0, y=-0.5, legend = c("S1", "cosine", "S2"), text.col=c("blue", "black", "red"), lty=1, col = c("blue", "black", "red"))
plot(angles*180/pi, corrs2, type="b", col="blue", pch = ".", xlim = c(95, 135), ylim = c(-0.75, -0.1), cex = 2,
main = "Three correlation functions",
xlab = "Angle (degrees)",
ylab = "Correlation")
lines(angles*180/pi, corrs3, type="b", col="red", pch = ".", cex = 2)
lines(angles*180/pi, cos(angles), type = "b", col = "black", pch = ".", cex = 2)
legend(x=95, y=-0.5, legend = c("S1", "cosine", "S2"), text.col=c("blue", "black", "red"), lty=1, col = c("blue", "black", "red"))
angles[16]*180/pi
[1] 112.5
round(corrs2[16], 4)
[1] -0.3728
round(cos(angles[16]), 4)
[1] -0.3827
round(corrs3[16], 4)
[1] -0.3906
round(1/sqrt(Ns2[16]), 4)
[1] 0.0012
round(1/sqrt(Ns3[16]), 4)
[1] 0.0013