richard — *Feb 17, 2014, 8:23 PM*

```
## Simultaneous simulation of the "Christian 2.0" model for the EPR-B correlations,
## in two versions.
## This version: **** sample size 100 million; just one angle ****
## "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(4567)
## 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 <- 105 * 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^8 ## 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) {
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
}
angles*180/pi
```

```
[1] 105
```

```
round(corrs2, 5)
```

```
[1] -0.2479
```

```
round(cos(angles), 5)
```

```
[1] -0.2588
```

```
round(corrs3, 5)
```

```
[1] -0.2625
```

```
round(1/sqrt(Ns2), 5)
```

```
[1] 0.00012
```

```
round(1/sqrt(Ns3), 5)
```

```
[1] 0.00013
```