library("dplyr")
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library("MASS")
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
set.seed(123456)
results <- rep(0, 100)
A <- c(rep(1, 50), rep(2, 50))
B <- c(rep(1, 25), rep(2, 25), rep(1, 25), rep(2, 25))
X <- rep(0, 100)
Y <- rep(0, 100)
N <- 0:15
bit1 <- N %%2
rem1 <- N %/% 2
bit2 <- rem1 %% 2
rem2 <- rem1 %/% 2
bit3 <- rem2 %% 2
rem3 <- rem2 %/% 2
bit4 <- rem3 %% 2
bit1 + bit2 * 2 + bit3 * 4 + bit4 * 8
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
X1 <- 2 * bit1 - 1
X2 <- 2 * bit2 - 1
Y1 <- 2 * bit3 - 1
Y2 <- 2 * bit4 - 1
S <- X1 * Y1 + X1 * Y2 + X2 * Y1 - X2 * Y2
lhvModel <- data.frame(X1, X2, Y1, Y2)[S == +2, ]
lhvModel
## X1 X2 Y1 Y2
## 1 -1 -1 -1 -1
## 3 -1 1 -1 -1
## 7 -1 1 1 -1
## 8 1 1 1 -1
## 9 -1 -1 -1 1
## 10 1 -1 -1 1
## 14 1 -1 1 1
## 16 1 1 1 1
for (i in 1:100){
LHVlong <- slice_sample(lhvModel, n = 100, replace = TRUE)
X[A == 1] <- LHVlong$X1[A == 1]
X[A == 2] <- LHVlong$X2[A == 2]
Y[B == 1] <- LHVlong$Y1[B == 1]
Y[B == 2] <- LHVlong$Y2[B == 2]
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.28 0.68 0.52 -0.60 2.08
## [1] 2.00 0.60 0.12 0.20 -0.68 1.60
## [1] 3.00 0.36 0.36 0.28 -0.44 1.44
## [1] 4.00 0.68 0.52 0.36 -0.44 2.00
## [1] 5.00 0.68 0.36 0.60 -0.20 1.84
## [1] 6.00 0.52 0.20 0.60 -0.60 1.92
## [1] 7.00 0.44 0.36 0.36 -0.44 1.60
## [1] 8.00 0.36 0.68 0.60 -0.44 2.08
## [1] 9.00 0.28 0.60 0.20 -0.60 1.68
## [1] 10.00 0.68 0.52 0.44 -0.44 2.08
## [1] 11.00 0.76 0.44 0.52 -0.60 2.32
## [1] 12.00 0.20 0.36 0.60 -0.60 1.76
## [1] 13.00 0.52 0.44 0.20 -0.44 1.60
## [1] 14.00 0.36 0.44 0.52 -0.52 1.84
## [1] 15.00 0.76 0.52 0.52 -0.52 2.32
## [1] 16.00 0.68 0.68 0.60 -0.28 2.24
## [1] 17.00 0.52 0.52 0.04 -0.20 1.28
## [1] 18.00 0.52 0.52 0.04 -0.60 1.68
## [1] 19.00 0.52 0.44 0.76 -0.52 2.24
## [1] 20.00 0.36 0.44 0.68 -0.44 1.92
## [1] 21.00 0.44 0.52 0.60 -0.52 2.08
## [1] 22.00 0.52 0.52 0.52 -0.12 1.68
## [1] 23.00 0.52 0.44 0.52 -0.52 2.00
## [1] 24.00 0.60 0.76 0.52 -0.44 2.32
## [1] 25.00 0.52 0.36 0.92 -0.44 2.24
## [1] 26.00 0.44 0.28 0.44 -0.60 1.76
## [1] 27.00 0.52 0.60 0.36 -0.36 1.84
## [1] 28.00 0.60 0.44 0.52 -0.60 2.16
## [1] 29.00 0.28 0.68 0.68 -0.60 2.24
## [1] 30.00 0.20 0.04 0.36 -0.12 0.72
## [1] 31.00 0.44 0.44 0.68 -0.36 1.92
## [1] 32.00 0.36 0.44 0.52 -0.44 1.76
## [1] 33.00 0.20 0.28 0.44 -0.68 1.60
## [1] 34.00 0.28 0.36 0.60 -0.36 1.60
## [1] 35.00 0.52 0.36 0.44 -0.52 1.84
## [1] 36.00 0.52 0.52 0.84 -0.20 2.08
## [1] 37.00 0.44 0.68 0.84 -0.52 2.48
## [1] 38.00 0.68 0.60 0.52 -0.44 2.24
## [1] 39.00 0.20 0.36 0.36 -0.76 1.68
## [1] 40.00 0.68 0.36 0.36 -0.68 2.08
## [1] 41.00 0.44 0.20 0.60 -0.12 1.36
## [1] 42.00 0.28 0.60 0.44 -0.44 1.76
## [1] 43.00 0.44 0.44 0.60 -0.76 2.24
## [1] 44.00 0.36 0.44 0.60 -0.76 2.16
## [1] 45.00 0.60 0.60 0.68 -0.44 2.32
## [1] 46.00 0.60 0.68 0.52 -0.60 2.40
## [1] 47.00 0.60 0.76 0.36 -0.36 2.08
## [1] 48.00 0.60 0.60 0.52 -0.44 2.16
## [1] 49.00 0.60 0.36 0.60 -0.44 2.00
## [1] 50.00 0.76 0.44 0.36 -0.52 2.08
## [1] 51.00 0.28 0.60 0.52 -0.60 2.00
## [1] 52.00 0.44 0.20 0.60 -0.68 1.92
## [1] 53.00 0.52 0.44 0.36 -0.60 1.92
## [1] 54.00 0.36 0.36 0.44 -0.68 1.84
## [1] 55.00 0.52 0.44 0.60 -0.60 2.16
## [1] 56.00 0.52 0.28 0.44 -0.52 1.76
## [1] 57.00 0.52 0.52 0.44 -0.44 1.92
## [1] 58.00 0.52 0.44 0.28 -0.68 1.92
## [1] 59.00 0.60 0.60 0.52 -0.76 2.48
## [1] 60.00 0.44 0.84 0.84 -0.36 2.48
## [1] 61.00 0.60 0.44 0.44 -0.36 1.84
## [1] 62.00 0.52 0.68 0.52 -0.52 2.24
## [1] 63.00 0.60 0.52 0.76 -0.52 2.40
## [1] 64.00 0.36 0.52 0.36 -0.44 1.68
## [1] 65.00 0.52 0.68 0.60 -0.44 2.24
## [1] 66.00 0.44 0.84 0.20 -0.52 2.00
## [1] 67.00 0.68 0.52 0.36 -0.52 2.08
## [1] 68.00 0.52 0.68 0.52 -0.36 2.08
## [1] 69.00 0.60 0.44 0.76 -0.68 2.48
## [1] 70.00 0.44 0.68 0.44 -0.44 2.00
## [1] 71.00 0.52 0.36 0.44 -0.28 1.60
## [1] 72.00 0.36 0.44 0.20 -0.44 1.44
## [1] 73.00 0.52 0.44 0.68 -0.36 2.00
## [1] 74.00 0.52 0.44 0.28 -0.44 1.68
## [1] 75.00 0.44 0.36 0.28 -0.92 2.00
## [1] 76.00 0.84 0.44 0.44 -0.60 2.32
## [1] 77.00 0.36 0.28 0.52 -0.52 1.68
## [1] 78.00 0.52 0.76 0.36 -0.44 2.08
## [1] 79.00 0.84 0.76 0.60 -0.52 2.72
## [1] 80.00 0.52 0.60 0.76 -0.60 2.48
## [1] 81.00 0.20 0.68 0.60 -0.68 2.16
## [1] 82.00 0.44 0.52 0.04 -0.60 1.60
## [1] 83.00 0.60 0.68 0.52 -0.52 2.32
## [1] 84.00 0.68 0.68 0.52 -0.44 2.32
## [1] 85.00 0.44 0.52 0.60 -0.28 1.84
## [1] 86.00 0.52 0.60 0.84 -0.60 2.56
## [1] 87.00 0.68 0.52 0.60 -0.60 2.40
## [1] 88.00 0.52 0.52 0.52 -0.60 2.16
## [1] 89.00 0.44 0.76 0.28 -0.44 1.92
## [1] 90.00 0.44 0.44 0.60 -0.68 2.16
## [1] 91.00 0.60 0.44 0.60 -0.68 2.32
## [1] 92.00 0.52 0.52 0.60 -0.84 2.48
## [1] 93.00 0.28 0.28 0.68 -0.36 1.60
## [1] 94.00 0.60 0.76 0.44 -0.36 2.16
## [1] 95.00 0.68 0.60 0.36 -0.36 2.00
## [1] 96.00 0.52 0.60 0.60 -0.28 2.00
## [1] 97.00 0.84 0.12 0.60 -0.60 2.16
## [1] 98.00 0.76 0.68 0.36 -0.52 2.32
## [1] 99.00 0.44 0.52 0.68 -0.36 2.00
## [1] 100.00 0.60 0.68 0.52 -0.52 2.32
truehist(results, breaks = seq(from = 0.5, to = 3, by = 0.1), prob = FALSE)
