toner_krivine.R

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

f0 <- function(x) 3 * (1 - 2 * x / pi) - 4 * (1 - 2 * x / pi)^3
f <- function(y){
    x <- abs(y)
    case_when(
        x < pi/4      ~ +1,
        x < pi/2      ~ f0(x),
        x < 3 *pi / 4 ~ -f0(pi - x),
        TRUE ~ -1
    )
}
n <- 1024
x <- seq(from = -pi, to = +pi, length = n)
plot(x, f(x), type = "l")

times <- 2 * pi * (0:(n-1))/n
fx <- f(x)
# b <- c(0.9284, -0.0597, -0.0041, 0.0017)


ks <- 0:10
cks <- (768 * cos((2 * ks + 1) * (pi / 4)) / (pi^4 * (2 * ks + 1)^3)) * ((1 / (2 * ks + 1)) - (-1)^ks * pi / 4 )
aks <- (-1)^ks * (2 / (pi * (2 * ks + 1))) * cks
bks <- rep(0, 11)
bks[1] <- 1 / (sqrt(2) * aks[1])
bks[2] <- - (1 / aks[1]) * aks[2] * bks[1]
bks[3] <- - (1 / aks[1]) * (aks[3] * bks[1])
bks[4] <- - (1 / aks[1]) * (aks[4] * bks[1])
bks[5] <- - (1 / aks[1]) * (aks[5] * bks[1] + aks[2] * bks[2])
bks[6] <- - (1 / aks[1]) * (aks[6] * bks[1])
bks[7] <- - (1 / aks[1]) * (aks[7] * bks[1])
bks[8] <- - (1 / aks[1]) * (aks[8] * bks[1] + aks[3] * bks[2] + aks[2] * bks[3])
bks[9] <- - (1 / aks[1]) * (aks[9] * bks[1])
bks[10] <- - (1 / aks[1]) * (aks[10] * bks[1])
bks[11] <- - (1 / aks[1]) * (aks[11] * bks[1] + aks[4] * bks[2] + aks[2] * bks[4])


b <- bks
absb <- abs(b)
absb<- c(absb, 1 - sum(absb))
b <- c(b, 1 - sum(abs(b)))
alpha <- 0
betas <- seq(-pi, pi, length = 1000)
thetas <- runif(1000000, -pi, pi)
yunifs <- runif(1000000)
corrs <- rep(0, 1000)
meanYs <- rep(0, 1000)
nn <- sample(1:12, 1000000, replace = TRUE, prob = absb)
X <- sign(b[nn]) * sign(cos((2 * nn - 1) * (thetas - alpha)))
mean(X)

## [1] -0.000264

for (i in 1:1000) {
beta <- betas[i]
index <- (2 * nn - 1) * (thetas - beta)
index <- 1024 * (index + 2 * pi) / (2 * pi)
index <- (round(index) %% 1024) + 1
meanBeta <- fx[index]
Y <- 2 * (yunifs < (1 + meanBeta)/2) - 1
corrs[i] <- mean(X*Y)
meanYs[i] <- mean(Y)
}
plot(betas, meanYs, type = "l")
abline(h = 0)

# jpeg()
plot(betas, corrs, type = "l", ylim = c(-1, +1))
abline(h = c(1 / sqrt(2), - 1 / sqrt(2)))
abline(h = c(-1, 0, 1))

# graphics.off()