mamas6.R

thetas <- sample(c(0, pi/2), 100000, replace = TRUE)
xunifs <- runif(100000)
yunifs <- runif(100000)
betas <- seq(from = 0, to = pi, length = 1000)
corrs <- rep(0, 1000)
for (alpha in seq(from = 0, to  = pi, length = 21)) {
  for(i in 1:1000){
    beta <- betas[i]
    x <- ifelse(xunifs < cos(thetas - alpha)^2, 1, -1)
    y <- ifelse(yunifs < cos(thetas + pi/2 - beta)^2, 1, -1)
    corrs[i] <- mean(x * y)
  }
  plot(betas * 180 / pi, corrs, type = "l", ylim = c(-1, +1), main = paste("Alpha = ", alpha *180 / pi, " degrees"))
  abline(h = c(-1, 0, 1))
}