richard — Feb 19, 2014, 8:06 AM
## Simulation of the "Christian 2.0" model for the EPR-B correlations,
## S^2 version (aka 3-D: think of S^2 as subset of R^3)
## This version optimized for speed at some slight cost to transparency.
## For other versions, and also some other models, take a look at
## http://rpubs.com/gill1109
## "Christian 2.0" refers to the second edition of Joy Christian's book, see
## http://libertesphilosophica.info/blog/
## The present algorithm has a number of inspirations: first of all
## what I like to call the S^1 version (2-D): Michel Fodje's algorithm:
## https://github.com/minkwe/epr-simple/
## The present S^2 version: exactly as Michel's, but with uniform random points
## on S^2 (subset R^3) instead of uniform random points on S^1 (subset R^2).
## This brings us closer to Christian's description of his own "three sphere" model.
## But it is also 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
## determined by theta_0.
## http://arxiv.org/abs/quant-ph/0210150
## http://freespace.virgin.net/ch.thompson1/
## http://freespace.virgin.net/ch.thompson1/Papers/The%20Record/TheRecord.htm
## One can think of this as a kind of "un-sharp" membership of the circular caps,
## an avenue Caroline herself proposed to explore.
## For another kind of un-sharp membership, see Chantal Roth's simple model
## http://rpubs.com/gill1109/13329 and
## http://rpubs.com/chenopodium (various versions EPR-B experiment, by Chantal)
## http://rpubs.com/gill1109 (yet more variations, in particular, Gisin-Gisin, by RDG)
set.seed(9875)
## 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 + 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 my measurement directions are all in the equatorial plane,
## I'll only generate z and x and treat them as x and y
## Measurement angles for setting "a": directions in the equatorial plane
angles <- seq(from = 0, to = 360, by = 7.5) * 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 ...
## I 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
theta <- runif(M, 0, pi/2) ## Joy and Minkwe's "theta_0"
s <- (sin(theta)^2) / 2
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.5, legend = c("Christian", "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.0),
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("Christian", "cosine"), text.col=c("blue", "black"), lty=1, col = c("blue", "black"))
angles[4]*180/pi # one particular angle
[1] 22.5
round(corrs[4], 4) # simulation
[1] 0.9091
round(cos(angles[4]), 4) # cosine
[1] 0.9239
round(1/sqrt(Ns[4]), 4) # standard error
[1] 0.0012