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

『ポアソン分布・ポアソン回帰・ポアソン過程』（島谷健一郎著, 近代科学社, 2017）の第4章「AICの根拠をシミュレーションで納得する」から、カルバック・ライブラー情報量や平均対数尤度のシミュレーションをRで再現してみました。

準備

作図のためggplot2パッケージを読み込みます。

また、擬似乱数を固定し、フォントも指定します。

library(ggplot2)

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

カルバック・ライブラー情報量

数列の違いを表す尺度として、差の平方和と、差の絶対値の和を計算する関数を定義します（式4.1, 4.2）。

dist_square <- function(p, q) {
  sum((p - q)^2)
}

dist_abs <- function(p, q) {
  sum(abs(p - q))
}

つづいて、カルバック・ライブラー情報量を計算する関数を定義します（式4.3）。

KL <- function(p, q) {
  sum(p * log(p / q))
}

強度4のポアソン分布と、強度3および5のポアソン分布とを比較します。まずは確率分布を図示します（図4.1(a)）。

X <- 0:20

p3 <- dpois(X, 3)
p4 <- dpois(X, 4)
p5 <- dpois(X, 5)

df <- data.frame(X = rep(X, 3),
           Y = c(p4, p3, p5),
           Intensity = factor(rep(c(4, 3, 5), each = length(X)),
                              levels = c(4, 3, 5)))
ggplot(df, aes(x = X, y = Y, colour = Intensity)) +
  geom_line(linewidth = 1) +
  geom_line(data = dplyr::filter(df, Intensity == 4), linewidth = 2) +
  geom_point(aes(shape = Intensity), size = 3) +
  scale_colour_discrete(name = "強度") +
  scale_shape(name = "強度") +
  theme_classic(base_family = jp_font)

カルバック・ライブラー情報量、平方和、差の絶対値の和をそれぞれ比較します。

KL(p4, p3)
## [1] 0.1507283
KL(p4, p5)
## [1] 0.1074258
dist_square(p4, p3)
## [1] 0.02237314
dist_square(p4, p5)
## [1] 0.01515742
dist_abs(p4, p3)
## [1] 0.4275235
dist_abs(p4, p5)
## [1] 0.3766872

カルバック・ライブラー情報量の値は図4.1(a)と一致するのですが、残りの2つは一致しません。どうも本のほうがなんだか間違っている気がするのですが（abs(dpois(1, 4) - dpois(1, 3)) = 0.0760986なのですが、本では0.13になっています）。

つづいて、強度4のポアソン分布と、強度3.3および4.8のポアソン分布の比較です。まず分布を図示します（図4.1(b)）。

p33 <- dpois(X, 3.3)
p48 <- dpois(X, 4.8)

df2 <- data.frame(X = rep(X, 3),
                  Y = c(p4, p33, p48),
                  Intensity = factor(rep(c("4", "3.3", "4.8"),
                                         each = length(X)),
                              levels = c("4", "3.3", "4.8")))
ggplot(df2, aes(x = X, y = Y, colour = Intensity)) +
  geom_line(linewidth = 1) +
  geom_line(data = dplyr::filter(df, Intensity == "4"),
            linewidth = 2) +
  geom_point(aes(shape = Intensity), size = 3) +
  scale_colour_discrete(name = "強度") +
  scale_shape(name = "強度") +
  theme_classic()

カルバック・ライブラー情報量、平方和、差の絶対値の和をそれぞれ比較します。

p33 <- dpois(X, 3.3)
p48 <- dpois(X, 4.8)

KL(p4, p33)
## [1] 0.06948756
KL(p4, p48)
## [1] 0.07071378
dist_square(p4, p33)
## [1] 0.01038763
dist_square(p4, p48)
## [1] 0.01011685
dist_abs(p4, p33)
## [1] 0.2937362
dist_abs(p4, p48)
## [1] 0.3051563

こちらも平方和の値が本と一致しません。

平均対数尤度

ポアソン分布モデルの平均対数尤度と最大対数尤度の比較をおこないます。まず、最大対数尤度と平均対数尤度を計算する関数を定義します。

mll_fun <- function(x, div = 1) {
  n <- length(x)
  xx <- split(x, ceiling(seq_along(x) / (n / div)))
  l <- sapply(1:div,
              function(i) {
                r <- mean(xx[[i]])
                sum(-r + xx[[i]] * log(r) - log(factorial(xx[[i]])))
              })
  return(sum(l))
}

all_fun <- function(x, lambda, div = 1, max_k = 50) {
  n <- length(x)
  xx <- split(x, ceiling(seq_along(x) / (n / div)))
  l <- sapply(1:div,
              function(i) {
                r <- mean(xx[[i]])
                n <- length(xx[[i]])
                k <- 0:max_k
                n * sum(dpois(k, lambda) * (-r + k * log(r) - log(factorial(k))))
              })
  return(sum(l))
}

つづいて、シミュレーションを実行する関数を定義します。サンプルサイズ(n)と強度(lambda)、パラメータ数を1にするかどうか、すなわち強度がすべて同じというモデルにするかどうか(n_param)を引数にとり（max_kは平均対数尤度の計算の上限値）、最大対数尤度と平均対数尤度の差を返します。

sim_fun <- function(n = 100, lambda = 4, one_param = FALSE, max_k = 50) {
  x <- sapply(lambda, function(r) rpois(n %/% length(lambda), r))
  
  if (!one_param) {
    # max log likelihood
    mll <- sapply(seq_along(lambda),
                  function(i) {
                    r <- mean(x[, i])
                    sum(-r + x[, i] * log(r) - log(factorial(x[, i])))
                  }) |>
      sum()

    # average log likelihood
    k <- 0:max_k
    all <- sapply(seq_along(lambda),
                  function(i) {
                    r <- mean(x[, i])
                    n <- length(x[, i])
                    n * sum(dpois(k, lambda[i]) *
                              (-r + k * log(r) - log(factorial(k))))
                }) |>
      sum()
  } else { # number of parameters = 1
    # max log likelihood
    r0 <- mean(x)
    mll <- sum(-r0 + x * log(r0) - log(factorial(x)))

    # average log likelihood
    k <- 0:max_k
    all <- sapply(seq_along(lambda),
                  function(i) {
                    n <- length(x[, i])
                    n * sum(dpois(k, lambda[i]) *
                              (-r0 + k * log(r0) - log(factorial(k))))
                }) |>
      sum()
  }
  return(mll - all)
}

シミュレーションごとの値を図示します（図4.2）。

R <- 100
sim <- replicate(R, sim_fun(n = 1000, lambda = 4))
data.frame(n = seq_len(R), y = sim) |>
  ggplot(aes(x = n, y = y)) +
  geom_point() +
  geom_hline(aes(yintercept = 0)) +
  labs(x = "シミュレーション番号", y = "最大対数尤度-平均対数尤度") +
  theme_bw(base_family = jp_font)

強度とサンプルサイズを変化させて、シミュレーションで差の平均を計算します（表4.2）。

result <- list()
R <- 10000
set <- data.frame(num_params = rep(c(2, 1), each = 6),
                  lambda_A = rep(4, 12),
                  lambda_B = c(rep(6, 3), rep(3, 3)),
                  size = rep(c(2000, 200, 20), 2))

result <- purrr::map(1:nrow(set),
            \(i) replicate(R,
                           sim_fun(n = set$size[i],
                                   lambda = c(set$lambda_A[i],
                                              set$lambda_B[i]),
                                   one_param = (set$num_params[i] == 1))) |>
              mean())

set |>
  dplyr::mutate(result = unlist(result) |>
                           round(2)) |>
  tidyr::pivot_wider(names_from = c("lambda_A", "lambda_B", "size"),
                     values_from = "result")
## # A tibble: 2 × 7
##   num_params `4_6_2000` `4_6_200` `4_6_20` `4_3_2000` `4_3_200` `4_3_20`
##        <dbl>      <dbl>     <dbl>    <dbl>      <dbl>     <dbl>    <dbl>
## 1          2       2.29      1.88     1.98       1.62      1.92     2.02
## 2          1       1.26      1.09     0.98       0.63      1.03     1.02

パラメータ数が2のときはおおよそ差は2、1のときの差はおおよそ1になりました。サンプルサイズが大きいほうが、結果を安定させるために必要なシミュレーションの繰り返し回数が大きくなるようです。

つづいて強度とサンプルサイズを変えて、パラメータ数3と1で比較します（表4.3(a)）。

result <- list()
R <- 10000
set <- data.frame(num_params = rep(c(1, 3), each = 2),
                  size = rep(c(30, 300), 2))

result <- purrr::map(1:nrow(set),
            \(i) replicate(R,
                           sim_fun(n = set$size[i],
                                   lambda = c(4, 5, 6),
                                   one_param = (set$num_params[i] == 1))) |>
              mean())


set |>
  dplyr::mutate(result = unlist(result) |>
                           round(2)) |>
  tidyr::pivot_wider(names_from = "size",
                     values_from = "result")
## # A tibble: 2 × 3
##   num_params  `30` `300`
##        <dbl> <dbl> <dbl>
## 1          1  0.98  1.02
## 2          3  3.07  2.91

パラメータ数が3のときの差はおよそ3になりました。

パラメータ数を1、2、3、4と変化させたときの差の値です（表4.3(b)）。

R <- 10000
result <- list()
set <- data.frame(one_param = rep(c(TRUE, FALSE), each = 3),
                  num_params = rep(c(2, 3, 4), 2))

result <- purrr::map(1:nrow(set),
            \(i) replicate(R,
                           sim_fun(n = 100,
                                   lambda = rep(4, set$num_params[i]),
                                   one_param = set$one_param[i])) |>
              mean())

set |>
  dplyr::mutate(result = unlist(result) |>
                           round(2) |>
                           as.character(),
                num_params = c(1, 1, 1, 2, 3, 4),
                size = rep(c(200, 300, 400), 2),
                .keep = "used") |>
  tidyr::pivot_wider(names_from = "size",
                     values_from = "result",
                     values_fill = "")
## # A tibble: 4 × 4
##   num_params `200`  `300`  `400` 
##        <dbl> <chr>  <chr>  <chr> 
## 1          1 "0.96" "1.01" "0.92"
## 2          2 "2.08" ""     ""    
## 3          3 ""     "3.11" ""    
## 4          4 ""     ""     "4.07"

おおよそパラメータ数に応じた値になっています。

結論

最大対数尤度 - 平均対数尤度 ≒ パラメータ数

となることが確認できました。