『ポアソン分布・ポアソン回帰・ポアソン過程』第2章の図をRで再現する

R
作者

伊東宏樹

公開

2025年4月3日

『ポアソン分布・ポアソン回帰・ポアソン過程』（島谷健一郎著, 近代科学社, 2017）の第2章「ポアソン分布モデルと最尤法」から、尤度関数の図と、モデルのあてはまり、最尤推定値の一致性・漸近正規性の図をRで再現してみました。

準備

作図のためggplot2とpatchworkパッケージを読み込みます。また、擬似乱数を固定します。フォントも設定します。

library(ggplot2)
library(patchwork)

set.seed(123)
jp_font <- "YuGothic"

データ

表2.2から、18個の調査区（5m×10m）におけるナナカマドとホオノキの本数のデータをそれぞれベクトルに格納します。調査区の面積(a)はareaに格納します。

nanakamado <- c(3, 3, 2, 3, 7, 10, 0, 0, 0,
                1, 1, 2, 0, 1, 1, 1, 4, 0)
hoonoki <- c(1, 2, 3, 3, 1, 0, 0, 2, 0,
             0, 0, 2, 0, 0, 3, 2, 0, 1)
area <- 0.5  # a (100 m^2)

まずはデータを確認します。このグラフは本にはありません。

data.frame(Site = factor(rep(1:18, 2)),
           N = c(nanakamado, hoonoki),
           Species = rep(c("ナナカマド", "ホオノキ"), each = 18)) |>
  ggplot(aes(x = Site, y = N, fill = Species)) +
  geom_col(position = "dodge") +
  scale_y_continuous(breaks = seq(0, 10, 2)) +
  theme_bw(base_family = jp_font)

図 2.3 尤度関数と対数尤度関数

ホオノキの本数データにポアソン分布をあてはめたときの尤度関数を定義します。引数は以下のとおりです。

  • r: 密度 (本/a)

  • X: 調査区内の本数

  • A: 調査区の面積 (a)

likelihood_fun <- function(r, X, A) {
  sapply(r,
         function(r_i)
           (exp(-r_i * A) * (r_i * A)^X / factorial(X)) |>
           prod()
  )
}

横軸を密度rとして、尤度関数を表示します。

data.frame(x = c(0, 6)) |>
  ggplot(aes(x = x)) +
  geom_function(fun = likelihood_fun,
                args = list(X = hoonoki, A = area),
                n = 201) +
  labs(title = "尤度関数", x = "r", y = "尤度") +
  theme_classic(base_family = jp_font)

同じく、対数尤度関数を定義します。

loglik_fun <- function(r, X, A) {
  sapply(r,
         function(r_i)
           (exp(-r_i * A) * (r_i * A)^X / factorial(X)) |>
           log() |>
           sum()
  )
}

横軸を密度rとして、対数尤度関数を表示します。

data.frame(x = c(0, 6)) |>
  ggplot(aes(x = x)) +
  geom_function(fun = loglik_fun,
                args = list(X = hoonoki, A = area),
                n = 201) +
  labs(title = "対数尤度関数", x = "r", y = "対数尤度") +
  theme_classic(base_family = jp_font)

最尤推定値を求めます。

mle <- optim(par = 1, fn = loglik_fun, method = "Brent",
      lower = 0, upper = 6,
      control = list(fnscale = -1),
      X = hoonoki, A = area)
print(mle$par)
## [1] 2.222222

密度の最尤推定値は2.222で、これはホオノキの密度mean(hoonoki) / area = 2.222と一致しました。

図 2.4 モデルのあてはまり

最尤推定値からポアソン分布にしたがう乱数を発生させ、それが観測データと矛盾しないかをみます。

ホオノキ

まず、ホオノキです。

最尤推定値は平均と一致することがわかっているので、平均から計算します。

ランダムに発生させた各調査区中の本数kの頻度分布を計算し、1000回のシミュレーションで中央値と95％区間を求めます。図ではそれぞれ点線と灰色の領域です。これに、観測値を実線と点で重ねます。

# simulation function
sim_fun <- function(n, lambda,
                    max_k = qpois(0.999, lambda)) {
  rpois(n, lambda) |>
    factor(levels = 0:max_k) |>
    table()
}

n <- 18
lambda <- mean(hoonoki)
R <- 1000
max_k <- 6
hoonoki_tbl <- hoonoki |>
  factor(levels = 0:max_k) |>
  table()

sim_hoo <- replicate(R, sim_fun(n, lambda, max_k)) |>
  apply(1, quantile, probs = c(0.025, 0.5, 0.975), type = 1)

data.frame(k = 0:max_k,
           obs = c(hoonoki_tbl),
           lower = sim_hoo["2.5%", ],
           median = sim_hoo["50%", ],
           upper = sim_hoo["97.5%", ]) |>
  ggplot() +
  geom_ribbon(aes(x = k, ymin = lower, ymax = upper),
              color = "gray", alpha = 0.3) +
  geom_line(aes(x = k, y = median), linetype = 2) +
  geom_line(aes(x = k, y = obs)) +
  geom_point(aes(x = k, y = obs), size = 2) +
  scale_x_continuous(name = "本数 k",
                     breaks = seq(0, max_k, 2)) +
  scale_y_continuous(name = "調査区の数",
                     breaks = seq(0, 12, 2)) +
  labs(title = "ホオノキ") +
  theme_classic(base_family = jp_font)

ホオノキでは95％の区間に収まりました。

ナナカマド

同様にナナカマドです。

lambda <- mean(nanakamado)
max_k <- 10
nanakamado_tbl <- nanakamado |>
  factor(levels = 0:max_k) |>
  table()

sim_nana <- replicate(R, sim_fun(n, lambda, max_k)) |>
  apply(1, quantile, probs = c(0.025, 0.5, 0.975), type = 1)

data.frame(k = 0:max_k,
                  obs = c(nanakamado_tbl),
                  lower = sim_nana["2.5%", ],
                  median = sim_nana["50%", ],
                  upper = sim_nana["97.5%", ]) |>
  ggplot() +
  geom_ribbon(aes(x = k, ymin = lower, ymax = upper),
              color = "gray", alpha = 0.3) +
  geom_line(aes(x = k, y = median), linetype = 2) +
  geom_line(aes(x = k, y = obs)) +
  geom_point(aes(x = k, y = obs), size = 2) +
  scale_x_continuous(name = "本数 k",
                     breaks = seq(0, max_k, 2)) +
  scale_y_continuous(name = "調査区の数",
                     breaks = 0:9, minor_breaks = NULL) +
  labs(title = "ナナカマド") +
  theme_classic(base_family = jp_font)

ナナカマドでは、k=10の調査区がシミュレーションの95％の区間からはずれています。

図 2.5 最尤推定値の一致性

強度4のポアソン分布にしたがう乱数を発生させて、サンプルサイズが大きくなるにつれて最尤推定値が真値に近づくことを確認します。

まずはn=4, 16, 64, 256のときです。

sim_fun <- function(R, n, lambda) {
  replicate(R, rpois(n, lambda)) |>
    apply(2, mean)
}

sizes <- 4^(1:7)
lambda <- 4
R <- 250

sim_result <- sapply(sizes,
                     function(n)
                       sim_fun(R, n, lambda)) |>
  data.frame()
level <- paste("n", sizes, sep = "=")
colnames(sim_result) <- level
cols <- palette.colors(n = 7, palette = "Okabe-Ito")

sim_long <- sim_result |>
  tidyr::pivot_longer(cols = 1:7) |>
  dplyr::mutate(name = factor(name, levels = level))

sim_long |>
  dplyr::filter(name %in% level[1:4]) |>
  ggplot(aes(x = value, colour = name)) +
  geom_line(stat = "bin", breaks = seq(0.2, 8.0, 0.4)) +
  scale_x_continuous(name = "最尤推定値",
                     breaks = seq(0.2, 8.0, 0.4),
                     minor_breaks = NULL) +
  scale_y_continuous(name = "回数") +
  scale_colour_manual(name = "サンプルサイズ",
                      values = cols[1:4]) +
  theme_classic(base_family = jp_font)

つづいてn=256, 1024, 4096, 16384の場合です。横軸のスケールが違うことに注意です。

sim_long |>
  dplyr::filter(name %in% level[4:8]) |>
  ggplot(aes(x = value, colour = name)) +
  geom_line(stat = "bin", breaks = seq(3.6, 4.4, 0.05)) +
  scale_x_continuous(name = "最尤推定値",
                     breaks = seq(3.6, 4.4, 0.05),
                     minor_breaks = NULL) +
  scale_y_continuous(name = "回数") +
  scale_colour_manual(name = "サンプルサイズ",
                      values = cols[4:7]) +
  theme_classic(base_family = jp_font)

図 2.6 最尤推定値の漸近正規性

サンプルサイズを大きくすると、最尤推定量が正規分布に近づくこと確認します。

前の図の作成に用いたデータのうち、n=4, 64, 1024, 16384について、データのヒストグラムと、対応する正規分布の確率密度関数を重ねて表示します。

plot_fun <- function(data, name, breaks) {
  n <- substring(name, 3) |>
    as.numeric()
  ggplot(data, aes(x = value)) +
  geom_histogram(breaks = breaks, fill = "gray") +
  geom_function(fun = ~ 250 * (breaks[2] - breaks[1]) *
                        dnorm(.x, 4, sqrt(4 / n))) +
  scale_x_continuous(breaks = breaks,
                     minor_breaks = NULL) +
  labs(title = name, x = "", y = "") +
  guides(x = guide_axis(angle = 45)) +
  theme_classic(base_family = jp_font)
}

breaks <- list("n=4" = seq(1.4, 7.0, 0.4),
               "n=64" = seq(2.6, 5.4, 0.2),
               "n=1024" = seq(3.65, 4.35, 0.05),
               "n=16384" = seq(3.93, 4.07, 0.01))
plots <- lapply(level[c(1, 3, 5, 7)],
                function(lv)
                  sim_long |>
                  dplyr::filter(name == lv) |>
                  plot_fun(lv, breaks = breaks[[lv]]))
(plots[[1]] + plots[[2]]) / (plots[[3]] + plots[[4]])