決めうちしないプログラム その2

どこぞやらにあるプログラム。指数分布とポアソン分布の関係を示そうとするシミュレーション・プログラムだ。

毎回1個ずつ乱数を発生させるのはもったいない。かといって，毎回何個発生させるのがベストなのか。取りあえずほどほどの個数を発生させておいて，足りなければ追加する。まとめて発生させるメリットと，いろいろな余分な操作のデメリットのせめぎ合いで，良さそうな所を見付ける。

オリジナルはこんな感じのプログラム

VisitorCounter <- function(lambda)
{
    counter <- 0
    time.arrival <- rexp(1, rate = lambda)
    while (time.arrival < 1) {
        counter <- counter + 1
        time.arrival <- time.arrival + rexp(1, rate = lambda)
    }
    return(counter)
}
lambda <- 5
N <- 10^5
x <- sapply(1:N, function(i) {VisitorCounter(lambda)})
barplot(rbind(table(x)/N, dpois(0:max(x), lambda)),
    col = c("violetred1", "slateblue4"),
    legend.text = c("Simulation", "Theoretical"),
    args.legend = list(x = "topright"),
    beside = TRUE)

書き換えると 5 倍弱の速さのプログラムになる

VisitorCounter2 <- function(N, lambda)
{
    result <- integer(N)
    repo <- matrix(rexp(15*N, rate = lambda), ncol=N)
    colsum <- colSums(repo)
    for (i in 1:N) {
        x <- repo[, i]
        if (colsum[i] < 1) {
            repeat {
                x <- c(x, rexp(10, rate = lambda))
                if (sum(x) > 1) break
            }
        }
        result[i] <- sum(cumsum(x) <= 1)
    }
    return(result)
}
lambda <- 5
N <- 10^5
x <- VisitorCounter2(N, lambda)
barplot(rbind(table(factor(x, levels=0:max(x)))/N, dpois(0:max(x), lambda)),
    col = c("violetred1", "slateblue4"),
    legend.text = c("Simulation", "Theoretical"),
    args.legend = list(x = "topright"),
    beside = TRUE)