### After a brief look at the Christensen et al data as reproduced in the
### original paper and recomputed from Graft's data sets, I turn to
### comparison of CHSH and CH in the ideal situation of a perfect Bell-CHSH
### experiment.

### Regarding the Christensen et al. data, 
### there are some small differences which need to be explained!

### For artifical (perfect) data from the ideal experiment,
### I 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 S is 50% better than B. ie 3/2 times as many standard
### deviations away from the local realist bound.

### Notation: before, "d" stood for detection, "n" for non-detection.
### So "dd" meant that both Alice and Bob have a detection ...
### But in the CHSH experiment "d" stands for "+1" and "n" for "-1".




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
result
##      setting
## count       11       12       21       22
##    Z  27089919 28192171 27663499 27633773
##    N  27153018 28352345 27827311 27926988
##    C     29173    34146    34473     1869
##    SA    46070    48077   150837   150713
##    SB    46202   146243    47448   144371
read.table("table1.dat")
##         V1       V2       V3       V4
## 1    46068    48076   150840   150505
## 2    29173    34145    34473     1862
## 3    46039   146205    47447   144070
## 4 27153020 28352350 27827318 27926994
table1 <- t(read.table("table1.dat"))

table1
##      [,1]  [,2]   [,3]     [,4]
## V1  46068 29173  46039 27153020
## V2  48076 34145 146205 28352350
## V3 150840 34473  47447 27827318
## V4 150505  1862 144070 27926994
rownames(table1) <- rownames(t(result)[ ,c("SA", "C", "SB", "N")])
colnames(table1) <- colnames(t(result)[ ,c("SA", "C", "SB", "N")])
table1
##        SA     C     SB        N
## 11  46068 29173  46039 27153020
## 12  48076 34145 146205 28352350
## 21 150840 34473  47447 27827318
## 22 150505  1862 144070 27926994
myTable1 <- t(result)[ ,c("SA", "C", "SB", "N")]
myTable1
##        count
## setting     SA     C     SB        N
##      11  46070 29173  46202 27153018
##      12  48077 34146 146243 28352345
##      21 150837 34473  47448 27827311
##      22 150713  1869 144371 27926988
table1
##        SA     C     SB        N
## 11  46068 29173  46039 27153020
## 12  48076 34145 146205 28352350
## 21 150840 34473  47447 27827318
## 22 150505  1862 144070 27926994
apply(table1, 2, sum)
##        SA         C        SB         N 
##    395489     99653    383761 111259682
apply(myTable1, 2, sum)
##        SA         C        SB         N 
##    395697     99661    384264 111259662
p <- c((1 + cos(pi/4))/4, (1 - cos(pi/4))/4, (1 - cos(pi/4))/4, (1 + cos(pi/4))/4)

p4 <- c((1 - cos(pi/4))/4, (1 + cos(pi/4))/4, (1 + cos(pi/4))/4, (1 - cos(pi/4))/4)
FixedProbs <- cbind(p, p, p, p4)

colnames(FixedProbs) <- c("11", "12", "21", "22")
rownames(FixedProbs) <- c("dd", "dn", "nd", "nn")

FixedProbs
##           11        12        21        22
## dd 0.4267767 0.4267767 0.4267767 0.0732233
## dn 0.0732233 0.0732233 0.0732233 0.4267767
## nd 0.0732233 0.0732233 0.0732233 0.4267767
## nn 0.4267767 0.4267767 0.4267767 0.0732233
Ns <- rep(25000000, 4)


AliceProbs1 <- AliceProbs2 <- BobProbs1 <- BobProbs2 <- c(0.5, 0.5)



B <- FixedProbs["dd", "11"] + FixedProbs["dd", "12"] + FixedProbs["dd", "21"] -
    (FixedProbs["dd", "22"] + AliceProbs1[1] + BobProbs1[1])
names(B) <- NULL
B
## [1] 0.2071068
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.828427
(S - 2) / 4
## [1] 0.2071068
Ns <- result["N", ]
seS <- sqrt(sum((1 - corrs^2)/Ns))
names(seS) <- NULL

(S - 2)/seS
## [1] 3089.057
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] 0.2071068
Bprime <- (C11  + C12  + C21  - C22 ) / 
                (S1dotA  + Sdot1B )
names(Bprime) <- NULL
        
Bprime      
## [1] 1.207107
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] 0.2071068
## [1] 0.000100766
## [1] 2055.324
## [1] 1.207107
S; (S - 2) / 4; (S - 2) / seS          
## [1] 2.828427
## [1] 0.2071068
## [1] 3089.057