## This simulation programme generates an N x 4 spreadsheet
## containing the numbers -1 and +1, and two lists of binary
## settings for Alice and Bob (setting 1 and 2 corresponding to
## outcomes of 2 N fair coin tosses) such that the
## quantum mechanical prediction of the value of S for
## the CHSH measurements on the singlet state is approximately 
## reproduced, when the settings are used to select outcomes
## and compute the experimentally observed CHSH quantity.

## The spreadsheet contains a set of N counterfactual outcomes
## for each of the N iterations of the experiment, ie the outcomes 
## which would be observed if either of Alice's two settings and 
## either of Bob's two settings had been chosen.

## Combined with the two sets of actually chosen settings we
## obtain the actual experimental record: for each iteration,
## a setting and an outcome for Alice and for Bob each.

## After this, we compute four observed correlations
## and the CHSH quantity S.

## The N x 4 spreadsheet of sets of counterfactual outcomes,
## and the N x 4 spreadsheet with the actually chosen settings
## and resulting experimental record, are both output as csv files.

## A little simulation experiment shows that the observed result was
## actually extremely lucky. Repeating 1000 times with new random
## settings, we never even see one violation of CHSH. However, with
## a different random seed, it could occasionally happen ...
## which is the point of the whole thing here.


## Reference:

## Gill, Richard D. "Statistics, Causality and Bell’s Theorem". 
## Statist. Sci. 29 (2014), no. 4, 512--528. doi:10.1214/14-STS490. 
## http://projecteuclid.org/euclid.ss/1421330545

## http://arxiv.org/abs/1207.5103
## https://pubpeer.com/publications/D985B475C637F666CC1D3E3A314522

## Files (one R script; two output files produced by the script):

## http://www.math.leidenuniv.nl/~gill/Spreadsheet/spreadsheet.R
## http://www.math.leidenuniv.nl/~gill/Spreadsheet/spreadsheet.csv
## http://www.math.leidenuniv.nl/~gill/Spreadsheet/experimentalRecord.csv

set.seed(1234)

N <- 800
S <- 2 * sqrt(2)
rho <- S / 4
quantumProbs <- c((1 + rho)/2, (1 - rho)/2, (1 - rho)/2, (1 + rho)/2)

settingsAlice <- sample(c(1, 2), N, replace = TRUE)
settingsBob <- sample(c(1, 2), N, replace = TRUE)

spreadsheet <- matrix(sample(c(-1, +1), N * 4, replace = TRUE), N, 4)

quantumOutcomes <- sample(4, N, quantumProbs, replace = TRUE)
quantumAlice <- 2 * (quantumOutcomes %in% c(1, 2)) - 1
quantumAlice[settingsAlice == 2 & settingsBob == 2] <- 
  - quantumAlice[settingsAlice == 2 & settingsBob == 2]
quantumBob <- 2 * (quantumOutcomes %in% c(1, 3)) - 1

spreadsheet[settingsAlice == 1, 1] <- quantumAlice[settingsAlice == 1]
spreadsheet[settingsAlice == 2, 2] <- quantumAlice[settingsAlice == 2]
spreadsheet[settingsBob == 1, 3] <- quantumBob[settingsBob == 1]
spreadsheet[settingsBob == 2, 4] <- quantumBob[settingsBob == 2]

colnames(spreadsheet) <- c("A", "A'", "B", "B'")
head(spreadsheet)
##       A A'  B B'
## [1,]  1 -1  1 -1
## [2,] -1 -1  1  1
## [3,] -1  1 -1 -1
## [4,]  1  1  1  1
## [5,] -1 -1 -1 -1
## [6,] -1  1 -1 -1
tail(spreadsheet)
##         A A'  B B'
## [795,] -1 -1 -1  1
## [796,]  1  1 -1  1
## [797,] -1  1  1  1
## [798,]  1  1  1 -1
## [799,]  1 -1  1  1
## [800,] -1  1 -1 -1
outcomesAlice <- integer(N)
outcomesBob <- integer(N)

for (i in 1:N) {
  outcomesAlice[i] <- spreadsheet[i, settingsAlice[i]]
  outcomesBob[i] <- spreadsheet[i, settingsBob[i] + 2]
}

experimentalRecord <- cbind(settingsAlice, outcomesAlice, settingsBob, outcomesBob)

head(experimentalRecord)
##      settingsAlice outcomesAlice settingsBob outcomesBob
## [1,]             1             1           1           1
## [2,]             2            -1           2           1
## [3,]             2             1           2          -1
## [4,]             2             1           1           1
## [5,]             2            -1           1          -1
## [6,]             2             1           2          -1
tail(experimentalRecord)
##        settingsAlice outcomesAlice settingsBob outcomesBob
## [795,]             2            -1           1          -1
## [796,]             1             1           2           1
## [797,]             2             1           1           1
## [798,]             2             1           2          -1
## [799,]             1             1           1           1
## [800,]             2             1           2          -1
rho11 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==1 & settingsBob ==1])
rho12 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==1 & settingsBob ==2])
rho21 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==2 & settingsBob ==1])
rho22 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==2 & settingsBob ==2])

rho11
## [1] 0.7628866
rho12
## [1] 0.6470588
rho21
## [1] 0.7232143
rho22
## [1] -0.7641026
S <- rho11 + rho12 + rho21 - rho22

S
## [1] 2.897262
2 * sqrt(2)
## [1] 2.828427
write.csv(spreadsheet, "spreadsheet.csv")
write.csv(experimentalRecord, "experimentalRecord.csv")


## Now we show that the results we just obtained were rather exceptional
## Keep the same spreadsheet but generate new settings, again and again

library(MASS)
Nrep <- 1000
results <- numeric(1000)
for (rep in 1:Nrep){
  settingsAlice <- sample(c(1, 2), N, replace = TRUE)
  settingsBob <- sample(c(1, 2), N, replace = TRUE)
  for (i in 1:N) {
    outcomesAlice[i] <- spreadsheet[i, settingsAlice[i]]
    outcomesBob[i] <- spreadsheet[i, settingsBob[i] + 2]
  }
  rho11 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==1 & settingsBob ==1])
  rho12 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==1 & settingsBob ==2])
  rho21 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==2 & settingsBob ==1])
  rho22 <- mean((outcomesAlice*outcomesBob)[settingsAlice ==2 & settingsBob ==2])
  results[rep] <- rho11 + rho12 + rho21 - rho22
}
truehist(results, xlim = c(0, 2))