## Comparison of CHSH and J for recent Bell experiments,
## together with optimally noise-reduced versions of both.
## Theory: https://pub.math.leidenuniv.nl/~gillrd/Peking/Peking_4.pdf
## In short: assume four multinomial samples,
## estimate covariance matrix of estimated relative frequencies,
## use sample deviations from no-signalling to optimally reduce
## the noise in the estimate of Bell's S or Eberhard's J
## AKA: generalized least squares
############# VIENNA #############
## The basic data, four 2x2 tables
table11 <- matrix(c(141439, 73391, 76224, 875392736),
2, 2, byrow = TRUE,
dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))
table12 <- matrix(c(146831, 67941, 326768, 874976534),
2, 2, byrow = TRUE,
dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))
table21 <- matrix(c(158338, 425067, 58742, 875239860),
2, 2, byrow = TRUE,
dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))
table22 <- matrix(c(8392, 576445, 463985, 874651457),
2, 2, byrow = TRUE,
dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))
table11
## Bob
## Alice d n
## d 141439 73391
## n 76224 875392736
table12
## Bob
## Alice d n
## d 146831 67941
## n 326768 874976534
table21
## Bob
## Alice d n
## d 158338 425067
## n 58742 875239860
table22
## Bob
## Alice d n
## d 8392 576445
## n 463985 874651457
## Check of the total number of trials
# "The number of valid trials is N = 3 502 784 150"
sum(table11) + sum(table12) + sum(table21) + sum(table22)
## [1] 3502784150
## The same data now in one 4x4 table
tables <- cbind(as.vector(t(table11)), as.vector(t(table12)), as.vector(t(table21)), as.vector(t(table22)))
dimnames(tables) = list(outcomes = c("++", "+-", "-+", "--"),
settings = c(11, 12, 21, 22))
tables
## settings
## outcomes 11 12 21 22
## ++ 141439 146831 158338 8392
## +- 73391 67941 425067 576445
## -+ 76224 326768 58742 463985
## -- 875392736 874976534 875239860 874651457
## The total number of trials for each setting pair
Ns <- apply(tables, 2, sum)
Ns
## 11 12 21 22
## 875683790 875518074 875882007 875700279
## observed relative frequencies, one 4x4 matrix
rawProbsMat <- tables / outer(rep(1,4), Ns)
rawProbsMat
## settings
## outcomes 11 12 21 22
## ++ 1.615183e-04 1.677076e-04 1.807755e-04 9.583188e-06
## +- 8.380993e-05 7.760091e-05 4.853017e-04 6.582675e-04
## -+ 8.704512e-05 3.732282e-04 6.706611e-05 5.298445e-04
## -- 9.996676e-01 9.993815e-01 9.992669e-01 9.988023e-01
## Convert the relative frequencies to one vector of length 16
VecNames <- as.vector(t(outer(colnames(rawProbsMat), rownames(rawProbsMat), paste, sep = "")))
rawProbsVec <- as.vector(rawProbsMat)
names(rawProbsVec) <- VecNames
VecNames
## [1] "11++" "11+-" "11-+" "11--" "12++" "12+-" "12-+" "12--" "21++" "21+-"
## [11] "21-+" "21--" "22++" "22+-" "22-+" "22--"
rawProbsVec
## 11++ 11+- 11-+ 11-- 12++
## 1.615183e-04 8.380993e-05 8.704512e-05 9.996676e-01 1.677076e-04
## 12+- 12-+ 12-- 21++ 21+-
## 7.760091e-05 3.732282e-04 9.993815e-01 1.807755e-04 4.853017e-04
## 21-+ 21-- 22++ 22+- 22-+
## 6.706611e-05 9.992669e-01 9.583188e-06 6.582675e-04 5.298445e-04
## 22--
## 9.988023e-01
## Building up the 4 no-signalling constraints, combined in one 16 x 4 matrix "NS"
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
## 11++ 1 0 1 0
## 11+- 1 0 0 0
## 11-+ 0 0 1 0
## 11-- 0 0 0 0
## 12++ -1 0 0 1
## 12+- -1 0 0 0
## 12-+ 0 0 0 1
## 12-- 0 0 0 0
## 21++ 0 1 -1 0
## 21+- 0 1 0 0
## 21-+ 0 0 -1 0
## 21-- 0 0 0 0
## 22++ 0 -1 0 -1
## 22+- 0 -1 0 0
## 22-+ 0 0 0 -1
## 22-- 0 0 0 0
## Build the 16x16 estimated covariance matrix of the 16 observed relative frequencies
cov11 <- diag(rawProbsMat[ , "11"]) - outer(rawProbsMat[ , "11"], rawProbsMat[ , "11"])
cov12 <- diag(rawProbsMat[ , "12"]) - outer(rawProbsMat[ , "12"], rawProbsMat[ , "12"])
cov21 <- diag(rawProbsMat[ , "21"]) - outer(rawProbsMat[ , "21"], rawProbsMat[ , "21"])
cov22 <- diag(rawProbsMat[ , "22"]) - outer(rawProbsMat[ , "22"], rawProbsMat[ , "22"])
Cov <- matrix(0, 16, 16)
rownames(Cov) <- VecNames
colnames(Cov) <- VecNames
Cov[1:4, 1:4] <- cov11/Ns["11"]
Cov[5:8, 5:8] <- cov12/Ns["12"]
Cov[9:12, 9:12] <- cov21/Ns["21"]
Cov[13:16, 13:16] <- cov22/Ns["22"]
## The vector "S" is used to compute the CHSH statistic "CHSH"
## The sum of the first three sample correlations minus the fourth
## Note: the experiment is designed to favour use of Eberhard's J !
S <- c(c(1, -1, -1 ,1), c(1, -1, -1 , 1), c(1, -1, -1, 1), - c(1, -1, -1, 1))
names(S) <- VecNames
CHSH <- sum(S * rawProbsVec)
CHSH
## [1] 2.000028
## Compute the estimated variance of the CHSH statistic,
## its estimated covariances with the observed deviations from no-signalling,
## and the 4x4 estimated covariance matrix of those deviations.
## We'll later also need the inverse of the latter.
varS <- t(S) %*% Cov %*% S
covNN <- t(NS) %*% Cov %*% NS
covSN <- t(S) %*% Cov %*% NS
covNS <- t(covSN)
InvCovNN <- solve(covNN)
## Estimated variance of the CHSH statistic,
## and estimated variance of the optimally "noise reduced" CHSH statistic.
varCHSH <- varS
varCHSHopt <- varS - covSN %*% InvCovNN %*% covNS
## The variance, and the improvement as ratio of standard deviations
varS
## [,1]
## [1,] 1.078084e-11
sqrt(varCHSH / varCHSHopt)
## [,1]
## [1,] 2.055586
## The coefficients of the noise reduced CHSH statistic and the resulting improved estimate
Sopt <- S - covSN %*% InvCovNN %*% t(NS)
Sopt
## 11++ 11+- 11-+ 11-- 12++ 12+- 12-+ 12--
## [1,] 2.011642 -0.5819319 -0.406426 1 2.302625 -1.418068 0.720693 1
## 21++ 21+- 21-+ 21-- 22++ 22+- 22-+ 22--
## [1,] 2.188951 0.7825251 -1.593574 1 -4.503218 -0.7825251 -0.720693 -1
CHSHopt <- sum(Sopt * rawProbsVec)
CHSHopt
## [1] 2.000028
## p-values assuming approximate normality for testing CHSH inequality
pnorm((CHSH - 2)/ sqrt(varCHSH), lower.tail = FALSE)
## [,1]
## [1,] 5.43703e-18
pnorm((CHSHopt - 2)/ sqrt(varCHSHopt), lower.tail = FALSE)
## [,1]
## [1,] 4.745794e-69
## Now we repeat for the Eberhard J statistic
## First, the coefficients in the vector "J"
## and the observed value of the statistic
J <- c(c(1, 0, 0 ,0), c(0, -1, 0 ,0), c(0, 0, -1, 0), c(-1, 0, 0, 0))
names(J) <- VecNames
sum(J * rawProbsVec)
## [1] 7.26814e-06
## Next, its estimated variance and resulting p-value
varJ <- t(J) %*% Cov %*% J
sum(J * rawProbsVec) / sqrt(varJ)
## [,1]
## [1,] 12.10426
pnorm(sum(J * rawProbsVec) / sqrt(varJ), lower.tail = FALSE)
## [,1]
## [1,] 5.013606e-34
## The covariances between J and the observed deviations from no-signaling
## The variance of the usual estimate of J and of the improved estimate of J
## The improvement as a ration of standard deviations
covJN <- t(J) %*% Cov %*% NS
covNJ <- t(covJN)
varJopt <- varJ - covJN %*% InvCovNN %*% covNJ
varJ
## [,1]
## [1,] 3.605539e-13
sqrt(varJ / varJopt)
## [,1]
## [1,] 1.503676
## The coefficients of an improved estimataor of Eberhard's J
Jopt <- J - covJN %*% InvCovNN %*% t(NS)
Jopt
## 11++ 11+- 11-+ 11-- 12++ 12+- 12-+
## [1,] 0.2529105 -0.395483 -0.3516065 0 0.3256562 -0.604517 -0.06982674
## 12-- 21++ 21+- 21-+ 21-- 22++ 22+-
## [1,] 0 0.2972378 -0.05436871 -0.6483935 0 -0.8758045 0.05436871
## 22-+ 22--
## [1,] 0.06982674 0
## Observed estimate of J, and improved estimate of J
sum(J * rawProbsVec)
## [1] 7.26814e-06
sum(Jopt * rawProbsVec)
## [1] 6.997615e-06
## p-values based on J and on improved J
## Note that the p-value based on improved J is the same as that of improved CHSH
pnorm(sum(J * rawProbsVec) / sqrt(varJ), lower.tail = FALSE)
## [,1]
## [1,] 5.013606e-34
pnorm(sum(Jopt * rawProbsVec) / sqrt(varJopt), lower.tail = FALSE)
## [,1]
## [1,] 4.745794e-69
## The p-value of the optimized J got a lot better, even though the estimate got a bit smaller
## The optimization procedure for CHSH made an enormous difference
## The deviation from no-signalling is small; it is responsible for these small changes