## 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

############# NIST #############

table11 <- matrix(c(6378, 3282, 3189, 43897356), 
    2, 2, byrow = TRUE,
    dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))
table12 <- matrix(c(6794, 2821, 23243, 43276943), 
    2, 2, byrow = TRUE,
    dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))
table21 <- matrix(c(6486, 21334, 2843, 43338281), 
    2, 2, byrow = TRUE,
    dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))    
table22 <- matrix(c(106, 27539, 30040, 42502788), 
    2, 2, byrow = TRUE,
    dimnames = list(Alice = c("d", "n"), Bob = c("d", "n")))

table11
##      Bob
## Alice    d        n
##     d 6378     3282
##     n 3189 43897356
table12
##      Bob
## Alice     d        n
##     d  6794     2821
##     n 23243 43276943
table21
##      Bob
## Alice    d        n
##     d 6486    21334
##     n 2843 43338281
table22
##      Bob
## Alice     d        n
##     d   106    27539
##     n 30040 42502788
## Check of the total number of trials

# "The number of valid trials is N = 173 149 423"
sum(table11) + sum(table12) + sum(table21) + sum(table22)
## [1] 173149423
## 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
##       ++     6378     6794     6486      106
##       +-     3282     2821    21334    27539
##       -+     3189    23243     2843    30040
##       -- 43897356 43276943 43338281 42502788
## The total number of trials for each setting pair

Ns <- apply(tables, 2, sum)
Ns
##       11       12       21       22 
## 43910205 43309801 43368944 42560473
## observed relative frequencies, one 4x4 matrix

rawProbsMat <- tables / outer(rep(1,4), Ns)
rawProbsMat
##         settings
## outcomes           11           12           21           22
##       ++ 1.452510e-04 1.568698e-04 1.495540e-04 2.490574e-06
##       +- 7.474345e-05 6.513537e-05 4.919188e-04 6.470558e-04
##       -+ 7.262549e-05 5.366684e-04 6.555382e-05 7.058192e-04
##       -- 9.997074e-01 9.992413e-01 9.992930e-01 9.986446e-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.452510e-04 7.474345e-05 7.262549e-05 9.997074e-01 1.568698e-04 
##         12+-         12-+         12--         21++         21+- 
## 6.513537e-05 5.366684e-04 9.992413e-01 1.495540e-04 4.919188e-04 
##         21-+         21--         22++         22+-         22-+ 
## 6.555382e-05 9.992930e-01 2.490574e-06 6.470558e-04 7.058192e-04 
##         22-- 
## 9.986446e-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.000092
## 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,] 2.47335e-10
sqrt(varCHSH / varCHSHopt)
##          [,1]
## [1,] 2.355303
## 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.017925 -0.4640094 -0.5180658    1 2.294644 -1.535991 0.8306344    1
##          21++      21+-      21-+ 21--      22++       22+-       22-+
## [1,] 2.336275 0.8182087 -1.481934    1 -4.648843 -0.8182087 -0.8306344
##      22--
## [1,]   -1
CHSHopt <- sum(Sopt * rawProbsVec)
CHSHopt
## [1] 2.000051
## p-values assuming approximate normality for testing CHSH inequality

pnorm((CHSH - 2)/ sqrt(varCHSH), lower.tail = FALSE)
##              [,1]
## [1,] 2.062969e-09
pnorm((CHSHopt - 2)/ sqrt(varCHSHopt), lower.tail = FALSE)
##              [,1]
## [1,] 1.110193e-14
## 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] 1.207121e-05
## Next, its estimated variance and resulting p-value

varJ <- t(J) %*% Cov %*% J
sum(J * rawProbsVec) / sqrt(varJ)
##          [,1]
## [1,] 4.778576
pnorm(sum(J * rawProbsVec) / sqrt(varJ), lower.tail = FALSE)
##              [,1]
## [1,] 8.827054e-07
## 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,] 6.381229e-12
sqrt(varJ / varJopt)
##          [,1]
## [1,] 1.513269
## 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.2544812 -0.3660023 -0.3795165    0 0.3236609 -0.6339977 -0.04234141
##      12--      21++        21+-       21-+ 21--       22++       22+-
## [1,]    0 0.3340686 -0.04544781 -0.6204835    0 -0.9122108 0.04544781
##            22-+ 22--
## [1,] 0.04234141    0
## Observed estimate of J, and improved estimate of J

sum(J * rawProbsVec)
## [1] 1.207121e-05
sum(Jopt * rawProbsVec)
## [1] 1.274879e-05
## 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,] 8.827054e-07
pnorm(sum(Jopt * rawProbsVec) / sqrt(varJopt), lower.tail = FALSE)
##              [,1]
## [1,] 1.110193e-14
## The p-value of the optimized J got a lot better, and the estimate got a bit bigger
## The optimization procedure for CHSH made an enormous difference
## The deviation from no-signalling is small; it is responsible for these small changes