Gisin Experiment, adapted from Richard Gill's code

Adding a plot

set.seed(4567)
N <- 10^4
nrsettings <- 10^3
z <- runif(N, -1, 1)
t <- runif(N, 0, 2 * pi)
r <- sqrt(1 - z^2)
x <- r * cos(t)
y <- r * sin(t)
e <- rbind(x, y, z)

whose <- sample(c(TRUE, FALSE), N, replace = TRUE)
unifs <- runif(N)

dangle <- 2.5
angles <- seq(from = 0, to = 180 - dangle, by = dangle) * 2 * pi/360
nrangles <- length(angles)
plotdata <- numeric(nrangles)
detect <- numeric(nrangles)

for (i in 1:nrsettings) {
    a <- rnorm(3)
    a <- a/sqrt(sum(a^2))
    b <- rnorm(3)
    b <- b/sqrt(sum(b^2))

    ca <- colSums(e * a)
    cb <- colSums(e * b)

    good <- (unifs < ifelse(whose, abs(ca), abs(cb)))
    detected <- sum(good)
    theta <- acos(sum(a * b)/(sqrt(sum(a * a)) * sqrt(sum(b * b))))
    angle <- theta * 180/pi
    corr <- sum(sign(ca[good]) * sign(cb[good]))/detected
    angleindex <- (angle + dangle)/dangle
    plotdata[angleindex] <- corr
    detect[angleindex] <- detected/N
    if (i%%nrsettings/100 == 0) 
        print(paste("Experiment ", i, ", angle=", angle, " detected: ", detected/N * 
            100, " actual corr ", round(sum(a * b), 4), ", observed corr ", 
            round(corr, 4)))
}

## [1] "Experiment  1000 , angle= 93.0297918170739  detected:  50.17  actual corr  -0.0529 , observed corr  -0.0636"

plotdata[nrangles] <- -plotdata[1]
detect[nrangles] <- detect[1]

Just adding a plot:

plotcolors <- c("black", "blue", "orange")
plot(angles * 180/pi, cos(angles), col = "black", pch = ".", cex = 2, main = "Gisin Experiment", 
    xlab = "Angle (degrees)", ylab = "Correlation")
lines(angles * 180/pi, cos(angles), col = "black", pch = ".", cex = 2)

points(angles * 180/pi, plotdata, type = "l", col = "blue", pch = ".", cex = 2)
points(angles * 180/pi, detect, type = "l", col = "orange", pch = ".", cex = 2)

legend(x = 0, y = -0.5, legend = c("cos", "correlation", "detected"), text.col = plotcolors, 
    lty = 1, col = plotcolors)