### 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 then compute the expected values both the Clauser-Horne statistic B
### and of the CHSH quantity S.
### I check that B = (S - 2) / 4
### Finally, using multinomial variance and covariance formulas I calculate the
### standard errors of the actually observed B and the actually observed S,
### and hence the expected amount of violation of local realism measured as
### number of standard deviations.
### It turns out that B is more than twice as good as S.
### Notation: "d" stands for detection, "n" for non-detection
### So "dd" means that both Alice and Bob have a detection ...
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
oddsRatios <- tables["dd", ] * tables["nn", ] / (tables["dn", ] * tables["nd", ])
logOddsRatios <- log(oddsRatios)
logOddsRatios
## 11 12 21 22
## 7.9181067 6.4239639 6.4482979 0.8899427
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)
Alpha11 <- AliceProbs1
Beta11 <- BobProbs1
gamma11 <- logOddsRatios["11"] / 4
Gamma11 <- matrix(0, 2, 2)
Gamma11[1, 1] <- Gamma11[2, 2] <- exp(gamma11)
Gamma11[1, 2] <- Gamma11[2, 1] <- exp(- gamma11)
for (i in 1:1000000){
Alpha11 <- AliceProbs1 / (Beta11 %*% Gamma11); Beta11 <- BobProbs1 / (Alpha11 %*% Gamma11)
}
Alpha12 <- AliceProbs1
Beta12 <- BobProbs2
gamma12 <- logOddsRatios["12"] / 4
Gamma12 <- matrix(0, 2, 2)
Gamma12[1, 1] <- Gamma12[2, 2] <- exp(gamma12)
Gamma12[1, 2] <- Gamma12[2, 1] <- exp(- gamma12)
for (i in 1:1000000){
Alpha12 <- AliceProbs1 / (Beta12 %*% Gamma12); Beta12 <- BobProbs2 / (Alpha12 %*% Gamma12)
}
Alpha21 <- AliceProbs2
Beta21 <- BobProbs1
gamma21 <- logOddsRatios["21"] / 4
Gamma21 <- matrix(0, 2, 2)
Gamma21[1, 1] <- Gamma21[2, 2] <- exp(gamma21)
Gamma21[1, 2] <- Gamma21[2, 1] <- exp(- gamma21)
for (i in 1:1000000){
Alpha21 <- AliceProbs2 / (Beta21 %*% Gamma21); Beta21 <- BobProbs1 / (Alpha21 %*% Gamma21)
}
Alpha22 <- AliceProbs2
Beta22 <- BobProbs2
gamma22 <- logOddsRatios["22"] / 4
Gamma22 <- matrix(0, 2, 2)
Gamma22[1, 1] <- Gamma22[2, 2] <- exp(gamma22)
Gamma22[1, 2] <- Gamma22[2, 1] <- exp(- gamma22)
for (i in 1:1000000){
Alpha22 <- AliceProbs2 / (Beta22 %*% Gamma22); Beta22 <- BobProbs2 / (Alpha22 %*% Gamma22)
}
P11 <- rep(Alpha11, 2) * as.vector(Gamma11) * rep(Beta11, c(2, 2))
P12 <- rep(Alpha12, 2) * as.vector(Gamma12) * rep(Beta12, c(2, 2))
P21 <- rep(Alpha21, 2) * as.vector(Gamma21) * rep(Beta21, c(2, 2))
P22 <- rep(Alpha22, 2) * as.vector(Gamma22) * rep(Beta22, c(2, 2))
Ns <- result["N", ]
FixedProbs <- cbind(P11, P12, P21, P22)
dimnames(FixedProbs) <- dimnames(tables)
FixedProbs
## 11 12 21 22
## dd 0.0010748891 0.0012051243 0.0012367450 0.0000669958
## dn 0.0006284474 0.0039586540 0.0004665914 0.0050967825
## nd 0.0006212894 0.0004910541 0.0041718065 0.0053415557
## nn 0.9976753742 0.9943451675 0.9941248570 0.9894946659
## Here are the expected counts under the new estimated model.
## i.e. no-signalling has been assumed and incorporated in the estimation.
FixedProbs * outer(rep(1, 4), Ns)
## 11 12 21 22
## dd 29186.48 34168.10 34415.29 1870.991
## dn 17064.24 112237.12 12983.99 142337.784
## nd 16869.88 13922.54 116090.16 149173.563
## nn 27089897.39 28192017.24 27663821.57 27633605.662
B <- FixedProbs["dd", "11"] + FixedProbs["dd", "12"] + FixedProbs["dd", "21"] -
(FixedProbs["dd", "22"] + AliceProbs1[1] + BobProbs1[1])
names(B) <- NULL
B
## [1] 5.024768e-05
corr11 <- 2 * (FixedProbs["dd", "11"] + FixedProbs["nn", "11"]) - 1
corr12 <- 2 * (FixedProbs["dd", "12"] + FixedProbs["nn", "12"]) - 1
corr21 <- 2 * (FixedProbs["dd", "21"] + FixedProbs["nn", "21"]) - 1
corr22 <- 2 * (FixedProbs["dd", "22"] + FixedProbs["nn", "22"]) - 1
corrs <- c(corr11, corr12, corr21, corr22)
signs <- c(1, 1, 1, -1)
S <- sum(corrs * signs)
names(S) <- NULL
S
## [1] 2.000201
(S - 2) / 4
## [1] 5.024768e-05
Ns <- result["N", ]
seS <- sqrt(sum((1 - corrs^2)/Ns))
names(seS) <- NULL
(S - 2)/seS
## [1] 3.699296
Ns <- result["N", ]
N11 <- Ns["11"]
N12 <- Ns["12"]
N21 <- Ns["21"]
N22 <- Ns["22"]
C11 <- FixedProbs["dd", "11"]
S11A <- AliceProbs1[1]
S11B <- BobProbs1[1]
C12 <- FixedProbs["dd", "12"]
S12A <- AliceProbs1[1]
S12B <- BobProbs2[1]
C21 <- FixedProbs["dd", "21"]
S21A <- AliceProbs2[1]
S21B <- BobProbs1[1]
C22 <- FixedProbs["dd", "22"]
S22A <- AliceProbs2[1]
S22B <- BobProbs2[1]
S1dotA <- AliceProbs1[1]
Sdot1B <- BobProbs1[1]
N1dot <- N11 + N12
Ndot1 <- N11 + N21
B <- (C11 + C12 + C21 ) -
(C22 + S1dotA +Sdot1B )
names(B) <- NULL
B
## [1] 5.024768e-05
Bprime <- (C11 + C12 + C21 - C22 ) /
(S1dotA + Sdot1B )
names(Bprime) <- NULL
Bprime
## [1] 1.014781
probs11 <- c(C11, S11A - C11, S11B - C11)
coefs11 <- c(1 - N11 / N1dot - N11 / Ndot1, - N11 / N1dot, - N11 / Ndot1)
covmat11 <- diag(probs11) - outer(probs11, probs11)
var11 <- sum(covmat11 * outer(coefs11, coefs11)) / N11
probs12 <- c(C12, S12A - C12)
coefs12 <- c(1 - N12 / N1dot, - N12 / N1dot)
covmat12 <- diag(probs12) - outer(probs12, probs12)
var12 <- sum(covmat12 * outer(coefs12, coefs12)) / N12
probs21 <- c(C21, S21B - C21)
coefs21 <- c(1 - N21 / Ndot1, - N21 / Ndot1)
covmat21 <- diag(probs21) - outer(probs21, probs21)
var21 <- sum(covmat21 * outer(coefs21, coefs21)) / N21
probs22 <- C22
var22 <- probs22 * (1 - probs22) / N22
seB <- sqrt(var11 + var12 + var21 + var22)
names(seB) <- NULL
B2 <- sum(probs11 * coefs11) + sum(probs12 * coefs12) + sum(probs21 * coefs21) - probs22
B; seB; B/seB; Bprime
## [1] 5.024768e-05
## [1] 6.582954e-06
## [1] 7.632999
## [1] 1.014781
S; (S - 2) / 4; (S - 2) / seS
## [1] 2.000201
## [1] 5.024768e-05
## [1] 3.699296