## First I write out the data for the "instruction sets".
## In a standard Bell experiment, Alice and Bob each choose one of two setting values.
## Bell's own simple local hidden variables model allows us to calculate,
## for any value of *his* hidden variable lambds, the values of the outcomes
## which Alice and Bob would observe, under each of those settings.
## It turns out that the resulting vector of four potential outcomes +/-1
## always takes one of 8 possible values, out of all the 16 = 2^4 theoretically
## possible. Moreover it turns out that each of the 8 has the same probability, 1/8.
## Knowing this, allows one to simulate an actual CHSH experiment, with "reality"
## following Bell's example model, using as hidden variable lambda = "which of the 8",
## chosen uniformly at random from the set {1, 2, ..., 8}
set.seed(12345678) ## Set random seed
library(MASS) ## Install a useful library
## To start with I write out the 8 vectors and do some simple checks.
x1 <- c(-1, -1, -1, -1, 1, 1, 1, 1)
x2 <- c(-1, -1, 1, 1, -1, -1, 1, 1)
y1 <- c(-1, -1, -1, 1, -1, 1, 1, 1)
y2 <- c(-1, 1, -1, -1, 1, 1, -1, 1)
x12 <- cbind(x1, x2)
y12 <- cbind(y1, y2)
s <- x1 * y1 + x1 * y2 + x2 * y1 - x2 * y2
s
## [1] 2 2 2 2 2 2 2 2
cbind(x1, x2, y1, y2)
## x1 x2 y1 y2
## [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
## [7,] 1 1 1 -1
## [8,] 1 1 1 1
mean(x1); mean(x2); mean(y1); mean(y2)
## [1] 0
## [1] 0
## [1] 0
## [1] 0
mean(x1 *y1); mean(x1 * y2); mean(x2 * y1); mean(x2 *y2)
## [1] 0.5
## [1] 0.5
## [1] 0.5
## [1] -0.5
round(table(x1, y1) / 8, 3)
## y1
## x1 -1 1
## -1 0.375 0.125
## 1 0.125 0.375
round(table(x1, y2) / 8, 3)
## y2
## x1 -1 1
## -1 0.375 0.125
## 1 0.125 0.375
round(table(x2, y1) / 8, 3)
## y1
## x2 -1 1
## -1 0.375 0.125
## 1 0.125 0.375
round(table(x2, y2) / 8, 3)
## y2
## x2 -1 1
## -1 0.125 0.375
## 1 0.375 0.125
## Next I set up the experiment I am going to do.
## Fix a sample size N & generate Alice and Bob's random settings, elements of {1, 2}
## The vectors X and Y which will later be filled with observed measurement
## outcomes are filled with zero's for the time being
N <- 10^4
X <- rep(0, N)
Y <- rep(0, N)
A <- sample(c(1, 2), N, replace = TRUE)
B <- sample(c(1, 2), N, replace = TRUE)
## This is where the physics starts:
## the local hidden variable takes values in the set {1, 2, ..., 8}
lambda <- sample(1:8, N, replace = TRUE)
table(A); table(B); table(lambda)
## A
## 1 2
## 4987 5013
## B
## 1 2
## 4931 5069
## lambda
## 1 2 3 4 5 6 7 8
## 1226 1242 1275 1272 1250 1228 1277 1230
## More physics: Alice's outcomes functions of her setting and of the LHV;
## Alice's outcomes functions of his setting and of the LHV;
## x12 and y12 are 8 x 2 matrices:
## rows indexed by values of lambda, columns by setting label
for (i in 1:N) X[i] <- x12[lambda[i], A[i]]
for (i in 1:N) Y[i] <- y12[lambda[i], B[i]]
## Now we bring the data together and computer four correlations
## and a value of S
corr11 <- mean(((X*Y)[A == 1 & B == 1]))
corr12 <- mean(((X*Y)[A == 1 & B == 2]))
corr21 <- mean(((X*Y)[A == 2 & B == 1]))
corr22 <- mean(((X*Y)[A == 2 & B == 2]))
corr11; corr12; corr21; corr22
## [1] 0.508576
## [1] 0.5024194
## [1] 0.4975248
## [1] -0.4978756
S <- corr11 + corr12 + corr21 - corr22
S
## [1] 2.006396
## I will now repeat the experiment 10 thousand times and draw a histogram of the
## values of S
Nrep <- 10^4
results <- rep(0, Nrep)
for (j in 1: Nrep) {
A <- sample(c(1, 2), N, replace = TRUE)
B <- sample(c(1, 2), N, replace = TRUE)
lambda <- sample(1:8, N, replace = TRUE)
for (i in 1:N) X[i] <- x12[lambda[i], A[i]]
for (i in 1:N) Y[i] <- y12[lambda[i], B[i]]
corr11 <- mean(((X*Y)[A == 1 & B == 1]))
corr12 <- mean(((X*Y)[A == 1 & B == 2]))
corr21 <- mean(((X*Y)[A == 2 & B == 1]))
corr22 <- mean(((X*Y)[A == 2 & B == 2]))
results[j] <- corr11 + corr12 + corr21 - corr22
}
truehist(results)
