Chantal2.R

richard — Feb 15, 2014, 8:30 PM

set.seed(1234) ## For reproducibility, testing!
##  You might like to change the random seed or let your computer
##  seed the pseudo random generator by its own dark means (system clock?)

## This is Chantal Roth's adaptation of my code of her model
## Chantal's additional features: compare different "fudge factors";
## compute CHSH; report rejection rates.

## Based on RDG's understanding of Chantal's java code at
## https://github.com/chenopodium/EPR

## Chantal says: better version at https://github.com/chenopodium/tango
## (but same simple model)

## This is still a "beta testing" version
## Main missing feature: any explanation of what is going on!

## For "Christian 2.0" models (and Caroline Thompson's rotating ball) see
## http://rpubs.com/gill1109/13325

## Chantal
plotcolors = c("black", "red", "green", "orange", "yellow", "gray", "cyan", "magenta")

dangle<-2.5
angles <- seq(from = 0, to = 350, by = dangle) * 2 * pi / 360
K <- length(angles)
corrs <- numeric(K)
Ns <- numeric(K)
beta <- 0 * 2 * pi / 360
M <- 10^4

theta <- runif(M, 0, 2*pi)

computeCHSH <- function(corrs) {
  # |R(A1B1) - R(A1B2) + R(A2B1) + R(A2B2)|, where the angles are A1=0, A2=90, B1=45 and B2=135 
  A1 <- 0
  A2 <- 90
  B1 <- 45
  B2 <- 135
  CHSH <- corrs[abs(B1-A1)/dangle] - corrs[abs(B2-A1)/dangle] + corrs[abs(A2 - B1)/dangle] + corrs[abs(B2-A2)/dangle]
  return (CHSH)
}
computeCorrs <- function(fudgefactor) {
  for (i in 1:(K-1)) {
    alpha <- angles[i]
    sA <- sin(theta + alpha)
    sB <- sin(theta + beta)
    good <- (runif(M) < fudgefactor * abs(sA)) & (runif(M) < fudgefactor * abs(sB))
    N <- sum(good)
    corrs[i] <<- sum(sign(sA[good])*sign(sB[good]))/N
    Ns[i] <<- N
  }
  corrs[K] <<- corrs[1]
  Ns[K] <<- Ns[1]     
}

runExperiment <- function(factors){   
  plot(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2,
       main = "Two correlation functions",
       xlab = "Angle (degrees)",
       ylab = "Correlation")
  lines(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2)        
  i <- 1
  means <- numeric(length(factors))
  CHSH <- numeric(length(factors))
  for (factor in factors) {
    result = computeCorrs(factor)

    means[i] <- round(mean(Ns/M)*100,1)
    CHSH[i] <- round(computeCHSH(corrs),2)
    print (CHSH[i])
    print ( means[i])        
    points(angles*180/pi, corrs, type="l", col=plotcolors[i+1], pch = ".", cex = 2)
    i <- i + 1    
    #  text(6, i*10, c("factor=",factor), cex = .8)
  }
  legend(x=0, y=-0.5, legend = c( "Cos", paste("Chantal, f=",factors, "det=",format(means, nsmall=0),"%", " CHSH=", CHSH)), text.col=plotcolors, lty=1, col = plotcolors)  
}

factors <- seq(from = 2.1, to = 2.7, by=0.2)
runExperiment(factors)

[1] 2.96
[1] 71.3
[1] 2.83
[1] 73.7
[1] 2.75
[1] 75.7
[1] 2.68
[1] 77.5