results <- rep(0, 100)
set.seed(1234)
for (i in 1:100){
data <- sample(1:4, 400, replace = TRUE, prob = c(3, 1, 1, 3))
X <- 1 * (data == 1 | data == 2) + (-1) * (data == 3 | data == 4)
Y <- 1 * (data == 1 | data == 3) + (-1) * (data == 2 | data == 4)
X[301:400] <- -X[301:400]
A <- c(rep(1, 200), rep(2, 200))
B <- c(rep(1, 100), rep(2, 100), rep(1, 100), rep(2, 100))
rho11 <- mean((X*Y)[A == 1 & B == 1])
rho12 <- mean((X*Y)[A == 1 & B == 2])
rho21 <- mean((X*Y)[A == 2 & B == 1])
rho22 <- mean((X*Y)[A == 2 & B == 2])
S <- rho11 + rho12 + rho21 - rho22
print(c(i, rho11, rho12, rho21, rho22, S))
results[i] <- S
}
## [1] 1.00 0.62 0.46 0.56 -0.48 2.12
## [1] 2.00 0.44 0.34 0.44 -0.50 1.72
## [1] 3.00 0.54 0.46 0.54 -0.66 2.20
## [1] 4.00 0.48 0.54 0.48 -0.66 2.16
## [1] 5.00 0.60 0.56 0.56 -0.44 2.16
## [1] 6.00 0.60 0.60 0.42 -0.46 2.08
## [1] 7.00 0.42 0.50 0.52 -0.44 1.88
## [1] 8.00 0.46 0.54 0.44 -0.52 1.96
## [1] 9.00 0.52 0.60 0.50 -0.36 1.98
## [1] 10.00 0.60 0.50 0.44 -0.58 2.12
## [1] 11.00 0.48 0.40 0.54 -0.40 1.82
## [1] 12.00 0.48 0.50 0.52 -0.58 2.08
## [1] 13.00 0.50 0.44 0.58 -0.46 1.98
## [1] 14.00 0.50 0.44 0.54 -0.50 1.98
## [1] 15.00 0.48 0.54 0.56 -0.56 2.14
## [1] 16.00 0.48 0.50 0.52 -0.44 1.94
## [1] 17.00 0.60 0.44 0.46 -0.48 1.98
## [1] 18.00 0.58 0.46 0.48 -0.64 2.16
## [1] 19.00 0.46 0.62 0.56 -0.54 2.18
## [1] 20.00 0.52 0.46 0.60 -0.40 1.98
## [1] 21.00 0.46 0.52 0.50 -0.32 1.80
## [1] 22.00 0.66 0.56 0.58 -0.44 2.24
## [1] 23.00 0.50 0.48 0.54 -0.48 2.00
## [1] 24.00 0.42 0.54 0.56 -0.36 1.88
## [1] 25.00 0.48 0.52 0.56 -0.60 2.16
## [1] 26.00 0.52 0.40 0.48 -0.64 2.04
## [1] 27.00 0.48 0.34 0.40 -0.48 1.70
## [1] 28.00 0.62 0.52 0.54 -0.40 2.08
## [1] 29.00 0.28 0.48 0.66 -0.60 2.02
## [1] 30.00 0.32 0.44 0.52 -0.70 1.98
## [1] 31.00 0.46 0.54 0.42 -0.50 1.92
## [1] 32.00 0.52 0.44 0.50 -0.66 2.12
## [1] 33.00 0.54 0.56 0.42 -0.40 1.92
## [1] 34.00 0.44 0.60 0.60 -0.46 2.10
## [1] 35.00 0.50 0.42 0.46 -0.52 1.90
## [1] 36.00 0.44 0.46 0.46 -0.48 1.84
## [1] 37.00 0.54 0.52 0.50 -0.38 1.94
## [1] 38.00 0.56 0.44 0.62 -0.48 2.10
## [1] 39.00 0.56 0.62 0.60 -0.56 2.34
## [1] 40.00 0.44 0.50 0.62 -0.46 2.02
## [1] 41.00 0.50 0.42 0.46 -0.48 1.86
## [1] 42.00 0.36 0.44 0.48 -0.52 1.80
## [1] 43.00 0.60 0.60 0.38 -0.34 1.92
## [1] 44.00 0.66 0.50 0.44 -0.52 2.12
## [1] 45.00 0.46 0.38 0.52 -0.48 1.84
## [1] 46.00 0.52 0.46 0.46 -0.60 2.04
## [1] 47.00 0.62 0.40 0.54 -0.46 2.02
## [1] 48.00 0.58 0.46 0.68 -0.44 2.16
## [1] 49.00 0.52 0.54 0.68 -0.56 2.30
## [1] 50.00 0.52 0.60 0.62 -0.46 2.20
## [1] 51.00 0.46 0.58 0.58 -0.34 1.96
## [1] 52.00 0.58 0.34 0.54 -0.52 1.98
## [1] 53.00 0.62 0.46 0.48 -0.52 2.08
## [1] 54.00 0.48 0.42 0.52 -0.28 1.70
## [1] 55.00 0.54 0.44 0.40 -0.40 1.78
## [1] 56.00 0.48 0.42 0.56 -0.46 1.92
## [1] 57.00 0.46 0.52 0.52 -0.30 1.80
## [1] 58.00 0.50 0.52 0.44 -0.52 1.98
## [1] 59.00 0.38 0.62 0.54 -0.54 2.08
## [1] 60.00 0.68 0.44 0.46 -0.42 2.00
## [1] 61.00 0.54 0.50 0.42 -0.52 1.98
## [1] 62.00 0.50 0.54 0.52 -0.36 1.92
## [1] 63.00 0.48 0.38 0.66 -0.62 2.14
## [1] 64.00 0.48 0.42 0.46 -0.54 1.90
## [1] 65.00 0.46 0.44 0.42 -0.60 1.92
## [1] 66.00 0.48 0.48 0.56 -0.64 2.16
## [1] 67.00 0.56 0.52 0.40 -0.50 1.98
## [1] 68.00 0.42 0.48 0.40 -0.50 1.80
## [1] 69.00 0.50 0.62 0.54 -0.48 2.14
## [1] 70.00 0.52 0.50 0.52 -0.38 1.92
## [1] 71.00 0.42 0.44 0.56 -0.56 1.98
## [1] 72.00 0.48 0.44 0.52 -0.50 1.94
## [1] 73.00 0.46 0.28 0.54 -0.60 1.88
## [1] 74.00 0.64 0.58 0.62 -0.64 2.48
## [1] 75.00 0.32 0.60 0.36 -0.46 1.74
## [1] 76.00 0.56 0.56 0.52 -0.54 2.18
## [1] 77.00 0.64 0.46 0.66 -0.46 2.22
## [1] 78.00 0.48 0.44 0.52 -0.48 1.92
## [1] 79.00 0.46 0.42 0.62 -0.50 2.00
## [1] 80.00 0.46 0.46 0.60 -0.60 2.12
## [1] 81.00 0.50 0.56 0.48 -0.44 1.98
## [1] 82.00 0.60 0.56 0.58 -0.44 2.18
## [1] 83.00 0.46 0.54 0.70 -0.44 2.14
## [1] 84.00 0.48 0.56 0.38 -0.46 1.88
## [1] 85.00 0.60 0.54 0.54 -0.54 2.22
## [1] 86.00 0.60 0.42 0.52 -0.54 2.08
## [1] 87.0 0.5 0.4 0.4 -0.3 1.6
## [1] 88.00 0.56 0.52 0.70 -0.42 2.20
## [1] 89.00 0.60 0.44 0.52 -0.48 2.04
## [1] 90.00 0.44 0.52 0.44 -0.50 1.90
## [1] 91.00 0.48 0.56 0.62 -0.42 2.08
## [1] 92.00 0.62 0.44 0.56 -0.50 2.12
## [1] 93.00 0.48 0.52 0.40 -0.40 1.80
## [1] 94.00 0.60 0.42 0.46 -0.52 2.00
## [1] 95.00 0.28 0.36 0.64 -0.58 1.86
## [1] 96.00 0.44 0.68 0.54 -0.42 2.08
## [1] 97.00 0.62 0.52 0.42 -0.60 2.16
## [1] 98.00 0.48 0.60 0.48 -0.36 1.92
## [1] 99.00 0.52 0.48 0.44 -0.44 1.88
## [1] 100.00 0.56 0.56 0.54 -0.56 2.22
library(MASS)
truehist(results, breaks = seq(from = 1.4, to = 2.8, by = 0.05), prob = FALSE)
for(j in 1 :3){
for (i in 1:100){
data <- sample(1:4, 400, replace = TRUE, prob = c(3, 1, 1, 3))
X <- 1 * (data == 1 | data == 2) + (-1) * (data == 3 | data == 4)
Y <- 1 * (data == 1 | data == 3) + (-1) * (data == 2 | data == 4)
X[301:400] <- -X[301:400]
A <- c(rep(1, 200), rep(2, 200))
B <- c(rep(1, 100), rep(2, 100), rep(1, 100), rep(2, 100))
rho11 <- mean((X*Y)[A == 1 & B == 1])
rho12 <- mean((X*Y)[A == 1 & B == 2])
rho21 <- mean((X*Y)[A == 2 & B == 1])
rho22 <- mean((X*Y)[A == 2 & B == 2])
S <- rho11 + rho12 + rho21 - rho22
results[i] <- S
}
truehist(results, breaks = seq(from = 1.4, to = 2.8, by = 0.05), prob = FALSE)
}
