## Comparison of CHSH and J for recent Bell experiments,
## together with optimally noise-reduced versions of both.
## Theory: https://arxiv.org/abs/2209.00702
## "Optimal statistical analyses of Bell experiments"
## 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
############# DIQKD #############
## Zhang, W., van Leent, T., Redeker, K. et al.
## A device-independent quantum key distribution system for distant users.
## Nature 607, 687–691 (2022). https://doi.org/10.1038/s41586-022-04891-y
## Data supplied by Tim van Leent
##########################
## The basic data, four 2x2 tables
setwd("/Users/richardgill/Desktop/Optimal Bell/Experiments/DIQKD")
rm(list = ls())
data <- read.csv("DIQKD_06-2020_2_psiP.csv")
attach(data)
summary(data)
## atom1_rnd_bit_1 atom1_rnd_bit_2 atom1_readout atom2_rnd_bit
## Min. :0.0000 Min. :0.0000 Min. : 6.00 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 23.00 1st Qu.:0.0000
## Median :1.0000 Median :0.0000 Median : 45.00 Median :0.0000
## Mean :0.5069 Mean :0.4934 Mean : 75.18 Mean :0.4925
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:129.00 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :186.00 Max. :1.0000
## atom2_readout coincidence photon_1 photon_2
## Min. : 0.00 Min. :0.000 Min. :750.0 Min. :750.2
## 1st Qu.: 3.00 1st Qu.:1.000 1st Qu.:753.8 1st Qu.:767.2
## Median : 7.00 Median :4.000 Median :758.7 Median :780.4
## Mean :17.76 Mean :2.683 Mean :762.3 Mean :784.8
## 3rd Qu.:33.00 3rd Qu.:4.000 3rd Qu.:767.1 3rd Qu.:798.7
## Max. :61.00 Max. :4.000 Max. :831.4 Max. :844.9
head(data)
## atom1_rnd_bit_1 atom1_rnd_bit_2 atom1_readout atom2_rnd_bit atom2_readout
## 1 1 0 87 1 3
## 2 0 0 10 1 1
## 3 1 1 92 1 2
## 4 0 1 95 1 22
## 5 0 0 96 1 37
## 6 0 0 116 1 40
## coincidence photon_1 photon_2
## 1 1 763.21 801.87
## 2 4 754.93 771.67
## 3 1 764.87 775.60
## 4 4 761.73 764.43
## 5 1 762.88 782.63
## 6 1 766.81 778.15
settingA <- 1 + atom1_rnd_bit_1 + 2 * atom1_rnd_bit_2
settingB <- 1 + atom2_rnd_bit
outcomeA <- ifelse(atom1_readout >= 62, +1, -1)
outcomeB <- ifelse(atom2_readout >= 12, +1, -1)
corrs <- matrix(0, 4, 2)
for (a in c(1, 2, 3, 4)) for (b in c(1, 2)) corrs[a, b] <- mean((outcomeA*outcomeB)[settingA == a & settingB == b])
ntrials <- matrix(0, 4, 2)
for (a in c(1, 2, 3, 4)) for (b in c(1, 2)) ntrials[a, b] <- sum(settingA == a & settingB == b)
table11 <- table(outcomeB[settingA == 4 & settingB == 1], -outcomeA[settingA == 4 & settingB == 1])
table12 <- table(outcomeB[settingA == 4 & settingB == 2], -outcomeA[settingA == 4 & settingB == 2])
table21 <- table(outcomeB[settingA == 3 & settingB == 1], -outcomeA[settingA == 3 & settingB == 1])
table22 <- table(outcomeB[settingA == 3 & settingB == 2], -outcomeA[settingA == 3 & settingB == 2])
table11
##
## -1 1
## -1 178 44
## 1 29 183
table12
##
## -1 1
## -1 199 36
## 1 28 160
table21
##
## -1 1
## -1 160 47
## 1 31 151
table22
##
## -1 1
## -1 38 160
## 1 166 39
## Check of the total number of trials
# "The number of valid trials is N = ... "
sum(table11) + sum(table12) + sum(table21) + sum(table22)
## [1] 1649
## 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
## ++ 178 199 160 38
## +- 44 36 47 160
## -+ 29 28 31 166
## -- 183 160 151 39
## The total number of trials for each setting pair
Ns <- apply(tables, 2, sum)
Ns
## 11 12 21 22
## 434 423 389 403
## observed relative frequencies, one 4x4 matrix
rawProbsMat <- tables / outer(rep(1,4), Ns)
rawProbsMat
## settings
## outcomes 11 12 21 22
## ++ 0.41013825 0.47044917 0.41131105 0.09429280
## +- 0.10138249 0.08510638 0.12082262 0.39702233
## -+ 0.06682028 0.06619385 0.07969152 0.41191067
## -- 0.42165899 0.37825059 0.38817481 0.09677419
## 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++ 12+- 12-+
## 0.41013825 0.10138249 0.06682028 0.42165899 0.47044917 0.08510638 0.06619385
## 12-- 21++ 21+- 21-+ 21-- 22++ 22+-
## 0.37825059 0.41131105 0.12082262 0.07969152 0.38817481 0.09429280 0.39702233
## 22-+ 22--
## 0.41191067 0.09677419
## 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
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.577832
## 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,] 0.005686275
sqrt(varCHSH / varCHSHopt)
## [,1]
## [1,] 1.000805
## 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,] 1.07514 -0.9606679 -0.9641923 1 0.9448948 -1.039332 -1.015773 1
## 21++ 21+- 21-+ 21-- 22++ 22+- 22-+ 22--
## [1,] 1.026324 -0.9378682 -1.035808 1 -1.046359 0.9378682 1.015773 -1
CHSHopt <- sum(Sopt * rawProbsVec)
CHSHopt
## [1] 2.577653
## p-values assuming approximate normality for testing CHSH inequality
pnorm((CHSH - 2)/ sqrt(varCHSH), lower.tail = FALSE)
## [,1]
## [1,] 9.096171e-15
pnorm((CHSHopt - 2)/ sqrt(varCHSHopt), lower.tail = FALSE)
## [,1]
## [1,] 8.831428e-15
## 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] 0.1510475
## Next, its estimated variance and resulting p-value
varJ <- t(J) %*% Cov %*% J
sum(J * rawProbsVec) / sqrt(varJ)
## [,1]
## [1,] 4.469809
pnorm(sum(J * rawProbsVec) / sqrt(varJ), lower.tail = FALSE)
## [,1]
## [1,] 3.914472e-06
## 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,] 0.001141956
sqrt(varJ / varJopt)
## [,1]
## [1,] 1.79399
## The coefficients of an improved estimataor of Eberhard's J
Jopt <- J - covJN %*% InvCovNN %*% t(NS)
Jopt
## 11++ 11+- 11-+ 11-- 12++ 12+- 12-+ 12--
## [1,] 0.01878495 -0.490167 -0.4910481 0 -0.0137763 -0.509833 -0.5039433 0
## 21++ 21+- 21-+ 21-- 22++ 22+- 22-+ 22--
## [1,] 0.00658103 -0.484467 -0.5089519 0 -0.01158968 0.484467 0.5039433 0
## Observed estimate of J, and improved estimate of J
sum(J * rawProbsVec)
## [1] 0.1510475
sum(Jopt * rawProbsVec)
## [1] 0.1444132
## 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,] 3.914472e-06
pnorm(sum(Jopt * rawProbsVec) / sqrt(varJopt), lower.tail = FALSE)
## [,1]
## [1,] 8.831428e-15
## Bell game
wins <- table11[1, 1] + table11[2, 2] +
table12[1, 1] + table12[2, 2] +
table21[1, 1] + table21[2, 2] +
table22[1, 2] + table22[2, 1]
Ntrials <- sum(table11) + sum(table12) + sum(table21) + sum(table22)
wins; Ntrials
## [1] 1357
## [1] 1649
1 - pbinom(wins, Ntrials, 3/4)
## [1] 5.151435e-13
## Our procedure for CHSH led to small improvement in p-value.
## CHSH is very close to optimised CHSH
## Because of the near perfect experimental symmetry,
## the statistical deviations from no-signalling are hardly correlated with S, at all.