EPRB23big.R

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