bell_game.R

###In this tiny Bell game I fill a very small data matrix completely at random
## both with settings 1 or 2 and with outcomes -1 or +1
## That is a model in which the four true correlations all equal zero,
## so the one sided CHSH statistic S = sum of three correlations minus fourth is zero too.
## Here, I take the fourth correlation to be the one corresponding to both settings = 1
############################################################################################
set.seed(12345) # Set random seed (for reproducibility)
N <- 10 # Sample size
A <- sample(c(1, 2), N, replace = TRUE) # Alice's settings - numerical vector of completely random choices of 1, 2
B <- sample(c(1, 2), N, replace = TRUE) # Bob's settings - numerical vector of completely random choices of 1, 2
X <- sample(c(-1, +1), N, replace = TRUE) # Alice's outcomes - numerical vector of completely random choices of -1, +1
Y <- sample(c(-1, +1), N, replace = TRUE) # Bob's outcomes - numerical vector of completely random choices of -1, +1
data <- cbind (A, B, X, Y) # Data matrix: N x 4, formed by pasting A, B, X, Y columnwise together
data # Display the data

##       A B  X  Y
##  [1,] 2 2  1 -1
##  [2,] 1 2 -1 -1
##  [3,] 2 2  1  1
##  [4,] 2 1  1 -1
##  [5,] 2 2 -1  1
##  [6,] 2 1  1  1
##  [7,] 2 2  1  1
##  [8,] 1 2  1 -1
##  [9,] 1 1  1  1
## [10,] 2 2  1  1

Z <- sum((X == Y) & (A == 2 | B == 2)) + sum((X != Y) & (A == 1 & B == 1)) # Count number of successes in Bell game
# One wants a big positive correlation for setting pairs (1,2), (2,1), (2,2),
#      a big negative correlation for setting pair  (1,1)
# Note: under LR the success probability can't exceed 0.75, under QM it can be 0.85
Z # Display the score

## [1] 5

# Display the p-value, the number of times a binomially (N, 3/4) variable 
pbinom(Z - 1, N, 0.75, lower.tail = FALSE)#could exceed or equal the observed score

## [1] 0.9802723

# export data set
# write.csv2(data, "bell_game.csv")

###########################################################################
#
# Now I'll simulate 10000 Bell games under both LR and QM,
#          each with 200 trials
#
###########################################################################
Nrep <- 10000
N <- 200
pValuesLR <- rep(0, Nrep)
pValuesQM <- rep(0, Nrep)
LRprobs <- c(0.375, 0.125, 0.125, 0.375)
q1 <- (1 + sqrt(2)/2)/4          # 0.4267767
q2 <- (1 - sqrt(2)/2)/4          # 0.0732233
QMprobs <- c(q1, q2, q2, q1)

pValue <- function(A, B, X, Y){
    Z <- sum((X == Y) & (A == 2 | B == 2)) + sum((X != Y) & (A == 1 & B == 1))
    pbinom(Z - 1, N, 0.75, lower.tail = FALSE)
}

for (r in 1:Nrep){  
    XY_LR <- sample(c("11", "12", "21", "22"), N, replace = TRUE, prob = LRprobs)
    XY_QM <- sample(c("11", "12", "21", "22"), N, replace = TRUE, prob = QMprobs)
    X_LR <- 2 * (XY_LR == "11" | XY_LR == "12") - 1
    Y_LR <- 2 * (XY_LR == "11" | XY_LR == "21") - 1
    X_QM <- 2 * (XY_QM == "11" | XY_QM == "12") - 1
    Y_QM <- 2 * (XY_QM == "11" | XY_QM == "21") - 1
    A <- sample(c(1, 2), N, replace = TRUE)
    B <- sample(c(1, 2), N, replace = TRUE)
    X_LR[A == 1 & B == 1] <- -X_LR[A == 1 & B == 1]
    X_QM[A == 1 & B == 1] <- -X_QM[A == 1 & B == 1]
    
    pValuesLR[r] <- pValue(A, B, X_LR, Y_LR)
    pValuesQM[r] <- pValue(A, B, X_QM, Y_QM)
}

library(MASS)

truehist(pValuesLR, prob = FALSE, breaks = (0:20)/20)

truehist(pValuesQM, prob = FALSE, breaks = (0:20)/20)

sum(pValuesLR < 0.05)/Nrep

## [1] 0.0419

sum(pValuesQM < 0.05)/Nrep

## [1] 0.9763