chantal.R

richard — Feb 15, 2014, 8:18 PM

set.seed(1234)

## This is a warm-up for Chantal2.R: http://rpubs.com/gill1109/13329
## Please see that script for explanation of what is going on here.

angles <- seq(from = 0, to = 360, by = 7.5) * 2 * pi / 360
K <- length(angles)
corrs <- numeric(K)
Ns <- numeric(K)

beta <- 0 * 2 * pi / 360

M <- 10^6

theta <- runif(M, 0, 2*pi)

for (i in 1:(K-1)) {
    alpha <- angles[i]
    sA <- sin(theta + alpha)
    sB <- sin(theta + beta)
    good <- (runif(M) < 2.3 * abs(sA)) & (runif(M) < 2.3 * abs(sB))
    N <- sum(good) 

    corrs[i] <- sum(sign(sA[good])*sign(sB[good]))/N
    Ns[i] <- N
}

corrs[K] <- corrs[1]
Ns[K] <- Ns[1]


plot(angles*180/pi, corrs, type="l", col="blue", pch = ".", cex = 2,
    main = "Two correlation functions",
    xlab = "Angle (degrees)",
    ylab = "Correlation")
points(angles*180/pi, corrs, col="blue", pch = ".", cex = 2)
lines(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2)
points(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2)

legend(x=0, y=-0.5, legend = c("Chantal", "cosine"), text.col=c("blue", "black"), lty=1, col = c("blue", "black"))

angles[16]*180/pi

[1] 112.5

round(corrs[16], 4)

[1] -0.3476

round(cos(angles[16]),  4)

[1] -0.3827

round(1/sqrt(Ns[16]), 4)

[1] 0.0012

Ns

 [1] 812092 801844 777572 751702 731718 722220 719060 719206 718623 719225
[11] 718489 718958 718730 718460 718871 718428 718417 718418 718825 722280
[21] 731941 750908 777464 801234 812209 801668 777694 750911 731320 722176
[31] 719164 719259 719245 718839 718597 718316 718614 718876 718470 718390
[41] 717863 719172 718800 722036 731656 751290 777766 801570 812092

Ns/M

 [1] 0.8121 0.8018 0.7776 0.7517 0.7317 0.7222 0.7191 0.7192 0.7186 0.7192
[11] 0.7185 0.7190 0.7187 0.7185 0.7189 0.7184 0.7184 0.7184 0.7188 0.7223
[21] 0.7319 0.7509 0.7775 0.8012 0.8122 0.8017 0.7777 0.7509 0.7313 0.7222
[31] 0.7192 0.7193 0.7192 0.7188 0.7186 0.7183 0.7186 0.7189 0.7185 0.7184
[41] 0.7179 0.7192 0.7188 0.7220 0.7317 0.7513 0.7778 0.8016 0.8121

angles*180/pi

 [1]   0.0   7.5  15.0  22.5  30.0  37.5  45.0  52.5  60.0  67.5  75.0
[12]  82.5  90.0  97.5 105.0 112.5 120.0 127.5 135.0 142.5 150.0 157.5
[23] 165.0 172.5 180.0 187.5 195.0 202.5 210.0 217.5 225.0 232.5 240.0
[34] 247.5 255.0 262.5 270.0 277.5 285.0 292.5 300.0 307.5 315.0 322.5
[45] 330.0 337.5 345.0 352.5 360.0

length(angles)

[1] 49