richard — *Feb 28, 2014, 7:41 PM*

```
## A simulation of Joy Christian's 3-sphere model for the EPR-Bohm
## correlations; see also http://libertesphilosophica.info/blog/
## S^2 version (aka 3-D: think of S^2 as subset of R^3)
## This version has been adapted from Richard Gill's optimized version of
## Michel Fodje's original simulation, http://rpubs.com/gill1109/EPRB3opt.
## All these simulations are inspired by the simulation of Joy Christian's
## 3-sphere model by Chantal Roth, https://github.com/chenopodium/JCS2.
## The present S^2 version is similar to Michel's, but with uniform random
## points on S^2 (subset R^3) instead of uniform random points on S^1.
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^7 ## Sample size. Next, try 10^7 or even 10^8 ...
## Follow Gill using 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
s <- rbeta(M, 1.29, 4.56) ## Zen's improved solution
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("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, 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("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), 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"))
```

```
plot(angles * 180/pi, corrs - cos(angles), type = "l")
```

```
max(abs(corrs - cos(angles)))
```

```
[1] 0.01505
```