Bimodal CLT, R - rextester

Language: Layout:

ppois(3,lambda=1.5,lower=FALSE)
1-ppois(3,1)
1-pbinom(0,12,1-ppois(3,1))

p <- 1 / 2
sigma <- 0.1
n <- 10

cutoff <- 1e-3
dfoo <- function(x) {
    (1 - p) * dnorm(x, 0, sigma) + p * dnorm(x, 1, sigma)
}
pfoo <- function(x) {
    (1 - p) * pnorm(x, 0, sigma) + p * pnorm(x, 1, sigma)
}
qfoo <- function(x) {
    lower <- 0
    upper <- 1
    k <- 1
    repeat {
        if (pfoo(lower) < x) break
        lower <- lower - k
        k <- k + 1
    }
    k <- 1
    repeat {
        if (pfoo(upper) > x) break
        upper <- upper + k
        k <- k + 1
    }
    phred <- function(y) pfoo(y) - x
    out <- uniroot(phred, lower = lower, upper = upper)
    out$root
}

from <- qfoo(cutoff / 2)
to <- qfoo(1 - cutoff / 2)
curve(dfoo(x), from = from, to = to, ylab = "f(x)",
    xlab = "x", main = "PDF of individual terms")

xx <- seq(0, n)
pp <- dbinom(xx, n, p)
dfoo <- function(x) {
    yy <- rep(NA, length(x))
    for (i in seq(along = x))
        yy[i] <- sum(pp * dnorm(x[i], xx, sqrt(n) * sigma))
    return(yy)
}
pfoo <- function(x) sum(pp * pnorm(x, xx, sqrt(n) * sigma))

from <- qfoo(cutoff / 2)
to <- qfoo(1 - cutoff / 2)
curve(dfoo(x), from = from, to = to, ylab = "f(x)",
    xlab = "x", n = 5101,
    main = "PDF of sum (black) and normal approx. (red)")
curve(dnorm(x, mean = n * p,
    sd = sqrt(n * p * (1 - p) + n * sigma^2)),
    col = "red", add = TRUE)

[ + ] Show input
Absolute running time: 1.04 sec, cpu time: 0.86 sec, memory peak: 33 Mb, absolute service time: 1,06 sec