どこぞやらにあるプログラム。指数分布とポアソン分布の関係を示そうとするシミュレーション・プログラムだ。
毎回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)
