### This script first of all takes the aggregate data from the Christensen et al
### experiment and then estimates the four 2x2 tables of probabilities under
### the assumption of no-signalling. Thus the marginal distributions match
### properly, but each 2x2 table has a log odds ratio which is equal to the
### log odds ratio observed in the experiment. Thus the two measurement outcomes
### have, according to the estimated probabilities,
### the same amount of "statistical dependence" as was observed in the
### experiment, but the estimates of their marginal distributions have been
### improved by assuming no-signalling and pooling the data from measurements
### with the same setting.
### I next compute the optimal test of local realism and the expected
### number of standard errors when each setting pair is measured the
### same number of times, which I take to be 28 000 000 times. I compare
### it with the CHSH test and the CH test (Christensen et al. version).
### Notation: "d" stands for detection, "n" for non-detection
### So "dd" means that both Alice and Bob have a detection ...
### Fitting the probabilities is done with the so-called
### "Iterative Proportional Fitting" method, which is
### asymptotically equivalent to maximum likelihood.
### See Little and Wu (1991) http://www.jstor.org/stable/2289718
### JASA 86, 87-95.
load("results.RData")
result <- apply(counts, c(1, 2), sum)
tables <- rbind(
"dd" = result["C", ],
"dn" = result["SA", ] - result["C", ],
"nd" = result["SB", ] - result["C", ],
"nn" = result["Z", ])
tables
## 11 12 21 22
## dd 29173 34146 34473 1869
## dn 16897 13931 116364 148844
## nd 17029 112097 12975 142502
## nn 27089919 28192171 27663499 27633773
probs <- tables / (outer(c(1, 1, 1, 1), result["N", ], ))
probs
## 11 12 21 22
## dd 0.0010743925 0.0012043448 0.0012388189 6.692451e-05
## dn 0.0006222881 0.0004913527 0.0041816473 5.329755e-03
## nd 0.0006271494 0.0039537118 0.0004662686 5.102663e-03
## nn 0.9976761699 0.9943505908 0.9941132652 9.895007e-01
AliceProbs1 <- (result["SA", "11"] + result["SA", "12"])/(result["N", "11"] + result["N", "12"])
AliceProbs2 <- (result["SA", "21"] + result["SA", "22"])/(result["N", "21"] + result["N", "22"])
BobProbs1 <- (result["SB", "11"] + result["SB", "21"])/(result["N", "11"] + result["N", "21"])
BobProbs2 <- (result["SB", "12"] + result["SB", "22"])/(result["N", "12"] + result["N", "22"])
AliceProbs1 <- c("d" = AliceProbs1, "n" = 1 - AliceProbs1)
AliceProbs2 <- c("d" = AliceProbs2, "n" = 1 - AliceProbs2)
BobProbs1 <- c("d" = BobProbs1, "n" = 1 - BobProbs1)
BobProbs2 <- c("d" = BobProbs2, "n" = 1 - BobProbs2)
P11mat <- t(matrix(tables[ , "11"], 2, 2)) / result["N", "11"]
dimnames(P11mat) <- list(alice = c("d", "n"), bob = c("d", "n"))
P11mat
## bob
## alice d n
## d 0.0010743925 0.0006222881
## n 0.0006271494 0.9976761699
P12mat <- t(matrix(tables[ , "12"], 2, 2)) / result["N", "12"]
dimnames(P12mat) <- list(alice = c("d", "n"), bob = c("d", "n"))
P12mat
## bob
## alice d n
## d 0.001204345 0.0004913527
## n 0.003953712 0.9943505908
P21mat <- t(matrix(tables[ , "21"], 2, 2)) / result["N", "21"]
dimnames(P21mat) <- list(alice = c("d", "n"), bob = c("d", "n"))
P21mat
## bob
## alice d n
## d 0.0012388189 0.004181647
## n 0.0004662686 0.994113265
P22mat <- t(matrix(tables[ , "22"], 2, 2)) / result["N", "22"]
dimnames(P22mat) <- list(alice = c("d", "n"), bob = c("d", "n"))
P22mat
## bob
## alice d n
## d 6.692451e-05 0.005329755
## n 5.102663e-03 0.989500658
fixProbs <- function(probs, alice, bob) {
for(i in 1: 10000) {
probs <- alice * probs / apply(probs, 1, sum)
pt <- t(probs)
probs <- t(bob * pt / apply(pt, 1, sum))
}
probs
}
P11 <- as.vector(t(fixProbs(P11mat, AliceProbs1, BobProbs1)))
P12 <- as.vector(t(fixProbs(P12mat, AliceProbs1, BobProbs2)))
P21 <- as.vector(t(fixProbs(P21mat, AliceProbs2, BobProbs1)))
P22 <- as.vector(t(fixProbs(P22mat, AliceProbs2, BobProbs2)))
Ns <- result["N", ]
FixedProbs <- c(P11, P12, P21, P22)
# dimnames(FixedProbs) <- dimnames(tables)
FixedProbs
## [1] 0.0010748891 0.0006212894 0.0006284474 0.9976753742 0.0012051243
## [6] 0.0004910541 0.0039586540 0.9943451675 0.0012367450 0.0041718065
## [11] 0.0004665914 0.9941248570 0.0000669958 0.0053415557 0.0050967825
## [16] 0.9894946659
## Here are the expected counts under the new estimated model.
## i.e. no-signalling has been assumed and incorporated in the estimation.
## Compare with the original
FixedProbsMat <- matrix(FixedProbs, 4, 4,
dimnames = list(outcomes = c("dd", "dn", "nd", "nn"),
settings = c(11, 12, 21, 22)))
FixedCountsMat <- FixedProbsMat * outer(rep(1, 4), Ns)
dimnames(FixedCountsMat) <- dimnames(tables)
FixedCountsMat
## 11 12 21 22
## dd 29186.48 34168.10 34415.29 1870.991
## dn 16869.88 13922.54 116090.16 149173.563
## nd 17064.24 112237.12 12983.99 142337.784
## nn 27089897.39 28192017.24 27663821.57 27633605.662
tables
## 11 12 21 22
## dd 29173 34146 34473 1869
## dn 16897 13931 116364 148844
## nd 17029 112097 12975 142502
## nn 27089919 28192171 27663499 27633773
### From now on we do a theoretical exercise in which the "true probabilities"
### are given by fixed probabilities which we just computed. For simplicity we
### will also take the number of trials at each setting combination to be the same.
### Because the true probabilities satisfy no-signalling, CHSH and CH give the same
### answer (when normalised properly). However when we "estimate" them by doing the
### experiment, the estimates have different standard errors. Hence the expected
### number of standard deviations departure from local realism can be different.
### Moreover, there are uncountably many alternative tests. I will compute the best
### one in the sense of having the smallest statistical error.
VecNames <- as.vector(t(outer(colnames(FixedProbsMat), rownames(FixedProbsMat), paste, sep = "")))
names(FixedProbs) <- VecNames
CHSH <- c(c(1, -1, -1, 1), c(1, -1, -1, 1), c(1, -1, -1, 1), c(-1, 1, 1, -1))
names(CHSH) <- VecNames
## value of CHSH
sum(CHSH * FixedProbs)
## [1] 2.000201
Aplus <- c(1, 1, 0, 0)
Aminus <- - Aplus
Bplus <- c(1, 0, 1, 0)
Bminus <- - Bplus
zero <- c(0, 0, 0, 0)
NSa1 <- c(Aplus, Aminus, zero, zero)
NSa2 <- c(zero, zero, Aplus, Aminus)
NSb1 <- c(Bplus, zero, Bminus, zero)
NSb2 <- c(zero, Bplus, zero, Bminus)
NS <- cbind(NSa1 = NSa1, NSa2 = NSa2, NSb1 = NSb1, NSb2 = NSb2)
rownames(NS) <- VecNames
NS
## NSa1 NSa2 NSb1 NSb2
## 11dd 1 0 1 0
## 11dn 1 0 0 0
## 11nd 0 0 1 0
## 11nn 0 0 0 0
## 12dd -1 0 0 1
## 12dn -1 0 0 0
## 12nd 0 0 0 1
## 12nn 0 0 0 0
## 21dd 0 1 -1 0
## 21dn 0 1 0 0
## 21nd 0 0 -1 0
## 21nn 0 0 0 0
## 22dd 0 -1 0 -1
## 22dn 0 -1 0 0
## 22nd 0 0 0 -1
## 22nn 0 0 0 0
cov11 <- diag(FixedProbsMat[ , "11"]) - outer(FixedProbsMat[ , "11"], FixedProbsMat[ , "11"])
cov12 <- diag(FixedProbsMat[ , "12"]) - outer(FixedProbsMat[ , "12"], FixedProbsMat[ , "12"])
cov21 <- diag(FixedProbsMat[ , "21"]) - outer(FixedProbsMat[ , "21"], FixedProbsMat[ , "21"])
cov22 <- diag(FixedProbsMat[ , "22"]) - outer(FixedProbsMat[ , "22"], FixedProbsMat[ , "22"])
Cov <- matrix(0, 16, 16)
rownames(Cov) <- VecNames
colnames(Cov) <- VecNames
Cov[1:4, 1:4] <- cov11
Cov[5:8, 5:8] <- cov12
Cov[9:12, 9:12] <- cov21
Cov[13:16, 13:16] <- cov22
varS <- t(CHSH) %*% Cov %*% CHSH
covSN <- t(CHSH) %*% Cov %*% NS
covNS <- t(covSN)
covNN <- t(NS) %*% Cov %*% NS
varS
## [,1]
## [1,] 0.08249738
N <- 28000000
seS <- sqrt(varS / N)
S <- sum(CHSH * FixedProbs)
## Number of standard deviations when we use plain CHSH
(S - 2) / seS
## [,1]
## [1,] 3.70284
covSN %*% solve(covNN) %*% covNS
## [,1]
## [1,] 0.06692506
seSopt <- sqrt((varS - covSN %*% solve(covNN) %*% covNS)/N)
S <- sum(CHSH * FixedProbs)
## Number of standard deviations when we use optimal test
(S - 2) / seSopt
## [,1]
## [1,] 8.522722
covSN %*% solve(covNN)
## NSa1 NSa2 NSb1 NSb2
## [1,] -0.5316002 -1.813628 -0.5856347 -1.802788
B <- c(c(1, 0, 0, 0), c(1, 0, 0, 0), c(1, 0, 0, 0), c(-1, 0, 0, 0))
B <- B - (c(Aplus, Aplus, zero, zero) + c(Bplus, zero, Bplus, zero)) / 2
sum(B * FixedProbs)
## [1] 5.024768e-05
varB <- t(B) %*% Cov %*% B / N
seB <- sqrt(varB)
seB
## [,1]
## [1,] 6.624193e-06
## Number of standard deviations when we use CH (Christensen et al version)
sum(B * FixedProbs) / seB
## [,1]
## [1,] 7.58548
### Check constraints
sum(c(Aplus, -Aplus, zero, zero) * FixedProbs)
## [1] 0
sum(c(zero, zero, Aplus, -Aplus) * FixedProbs)
## [1] -1.355253e-20
sum(c(Bplus, zero, -Bplus, zero) * FixedProbs)
## [1] -1.084202e-19
sum(c(zero, Bplus, zero, -Bplus) * FixedProbs)
## [1] 2.032879e-19