R: 『ポアソン分布・ポアソン回帰・ポアソン過程』のグラフをRで再現してみる(2)

つづきです。第1章ののこり。

長さ10000のなかに40000個の点をランダムに配置し、0と1の間に はいっている個数をかぞえるという設定です。

library(ggplot2)
library(dplyr)
set.seed(1)

# 図1.5
sim3 <- function(R = 10000, N = 40000, L = 10000, M = 1) {
  x <- replicate(R, sum(runif(N, 0, L) < M))

  types <- factor(c("simulated", "expected"),
                  levels = c("simulated", "expected"))
  # simulated
  df <- data.frame(table(factor(x, levels = 0:max(x)))) %>%
    transmute(x = as.numeric(Var1) - 1,
           dens = Freq / R,
           type = types[1])
  # expected
  df2 <- data.frame(x = df$x) %>%
    mutate(dens = dpois(x, N * M / L),
           type = types[2]) %>%
    bind_rows(df)

  px <- 0:max(x)
  py <- pretty(df2$dens)
  df2 %>%
    ggplot(mapping = aes(x = x, y = dens, fill = type)) +
    geom_bar(position = "dodge", stat = "identity", colour = "black") +
    scale_fill_manual(values = c("black", "white")) +
    labs(x = paste("長さ", M, "の部分に入った個数"), y = "割合") +
    scale_x_continuous(breaks = px, labels = px) +
    scale_y_continuous(limits = range(py), expand = c(0, 0)) +
    theme_classic(base_family = "IPAexGothic")
}

sim3(10000, 40000, 10000, 1)

図1.5
白い部分は、ポアソン分布による確率。

(a) 長さ10000のなかに20000回のイベント、長さ1のなかに はいった個数という設定。(b) 長さ10000のなかに40000回のイベントで、長さ2の間に はいった個数という設定。

# 図1.6
sim3(10000, 20000, 10000, 1)
sim3(10000, 40000, 10000, 2)

図1.6(a)

図1.6(b)

頻度4でイベントが発生するとした場合に、全体を10000等分したときのイベント発生回数の合計。

# 図1.7
sim4 <- function(R = 200, N = 10000, r = 4) {
  rh <- r / N
  x <- replicate(R, sum(rbinom(N, 1, rh)))

  max.x <- max(x, qpois(0.999, r))
  types <- factor(c("simulated", "expected"),
                  levels = c("simulated", "expected"))
  # simulated
  df <- data.frame(table(factor(x, levels = 0:max.x))) %>%
    transmute(x = as.numeric(Var1) - 1,
           dens = Freq / R,
           type = types[1])
  # expected
  df2 <- data.frame(x = df$x) %>%
    mutate(dens = dpois(x, r),
           type = types[2]) %>%
    bind_rows(df)

  px <- 0:max.x
  py <- pretty(df2$dens)
  df2 %>%
    ggplot(mapping = aes(x = x, y = dens, fill = type)) +
    geom_bar(position = "dodge", stat = "identity", colour = "black") +
    scale_fill_manual(values = c("black", "white")) +
    labs(x = "全体での総数", y = paste(R, "回中の割合")) +
    scale_x_continuous(breaks = px, labels = px) +
    scale_y_continuous(limits = range(py), expand = c(0, 0)) +
    theme_classic(base_family = "IPAexGothic")
}

sim4(200, 10000, 4)

図1.7

ポアソン分布にしたがう乱数を発生させる。

# 図1.8
sim5 <- function(R = 200, r = 4, max = 20) {
  p <- ppois(0:max, r)
  x <- runif(R, 0, 1)
  n <- sapply(x, function(i) sum(i > p))

  max.n <- max(n, qpois(0.9999, r))
  types <- factor(c("simulated", "expected"),
                  levels = c("simulated", "expected"))
  # simulated
  df <- data.frame(table(factor(n, levels = 0:max.n))) %>%
    transmute(x = as.numeric(Var1) - 1,
              freq = Freq,
              type = types[1])
  # expected
  df2 <- data.frame(x = df$x) %>%
    mutate(freq = R * dpois(x, r),
           type = types[2]) %>%
    bind_rows(df)

  px <- 0:max.n
  py <- pretty(df2$freq)
  df2 %>%
    ggplot(mapping = aes(x = x, y = freq, fill = type)) +
    geom_bar(position = "dodge", stat = "identity", colour = "black") +
    scale_fill_manual(values = c("black", "white")) +
    labs(x = "生成された乱数", y = "生成された回数") +
    scale_x_continuous(breaks = px, labels = px) +
    scale_y_continuous(limits = range(py), expand = c(0, 0)) +
    theme_classic(base_family = "IPAexGothic")
}

sim5(10000)

図1.8

タグ：R