## First simulation: S = 2.4, N = 250

set.seed(1234)
probs <- c(0.4, 0.1, 0.1, 0.4)
corrs <- numeric(4)
results <- numeric(1000)
for (i in 1:1000){
    outcomes <- sample(4, 250, replace = TRUE, prob = probs)
    settings <- sample(4, 250, replace = TRUE)
    for (set in c(1:4)){
        corrs[set] <- (sum(outcomes %in% c(1, 4) & settings == set) - 
            sum(outcomes %in% c(2, 3) & settings == set)) /
            sum(settings == set)}
    results[i] <- sum(corrs)
    }
library(MASS)
truehist(results, main = "True S = 2.4, N = 250", xlab = "Observed value of S")

sd(results)
## [1] 0.1992223
sd(results)
## [1] 0.1992223
sum(results > 2.4)/1000
## [1] 0.517
## Second simulation: S = 2, N = 250

set.seed(1234)
probs <- c(0.375, 0.125, 0.125, 0.375)
corrs <- numeric(4)
results <- numeric(1000)
for (i in 1:1000){
    outcomes <- sample(4, 250, replace = TRUE, prob = probs)
    settings <- sample(4, 250, replace = TRUE)
    for (set in c(1:4)){
        corrs[set] <- (sum(outcomes %in% c(1, 4) & settings == set) - 
            sum(outcomes %in% c(2, 3) & settings == set)) /
            sum(settings == set)}
    results[i] <- sum(corrs)
    }
library(MASS)
truehist(results, main = "True S = 2, N = 250", xlab = "Observed value of S")

sd(results)
## [1] 0.2209447
sum(results > 2.4)/1000
## [1] 0.025
## Third simulation: S = 2 sqrt 2, N = 250

set.seed(1234)
probs <- c(0.5 - (1 - sqrt(2)/2)/4, (1 - sqrt(2)/2)/4, (1 - sqrt(2)/2)/4,  0.5 - (1 - sqrt(2)/2)/4)
corrs <- numeric(4)
results <- numeric(1000)
for (i in 1:1000){
    outcomes <- sample(4, 250, replace = TRUE, prob = probs)
    settings <- sample(4, 250, replace = TRUE)
    for (set in c(1:4)){
        corrs[set] <- (sum(outcomes %in% c(1, 4) & settings == set) - 
            sum(outcomes %in% c(2, 3) & settings == set)) /
            sum(settings == set)}
    results[i] <- sum(corrs)
    }
library(MASS)
truehist(results, main = "True S = 2 sqrt 2, N = 250", xlab = "Observed value of S")

sd(results)
## [1] 0.1731038
sum(results > 2.4)/1000
## [1] 0.989