はじめに

　機械学習や統計学で登場する各種の確率分布について、「計算式の導出・計算のスクラッチ実装・計算過程や結果の可視化」などの「数式・プログラム・図」を用いた解説により、様々な角度から理解を目指すシリーズです。

　この記事では、混合ポアソン分布の乱数についてR言語を使って確認します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
混合ポアソン分布の乱数生成
- サンプリング
- サンプルサイズの影響
  - サンプルサイズと形状の関係
参考文献
おわりに

混合ポアソン分布の乱数生成

　混合ポアソン分布(mixture Poisson distribution)に従う乱数を生成して、グラフやアニメーションで確認します。この記事では、Rのggplot2パッケージを利用して作図します。
　ポアソン分布については「【R】ポアソン分布の乱数生成 - からっぽのしょこ」、混合ポアソン分布については「混合ポアソン分布の定義式の導出 - からっぽのしょこ」、Rによるグラフ作成については「【R】混合ポアソン分布の作図 - からっぽのしょこ」、Pythonを利用する場合は「Python版」を参照してください。

　利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)

　この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードについては(ごちゃごちゃしないように)パッケージ名を省略するため、ggplot2 を読み込む必要があります。
　また、ネイティブパイプ(baseパイプ)演算子 |> を使います。(基本的には)magrittrパッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

サンプリング

　まずは、混合ポアソン分布に従う乱数を生成します。

パラメータの設定

　混合ポアソン分布のパラメータ $\boldsymbol{\lambda}$ を設定します。

# クラスタごとのパラメータを指定
lambda_k <- c(1, 5.5, 10, 16.8)
lambda_k

[1]  1.0  5.5 10.0 16.8

　各クラスタ(クラス)に対応する $K$ 個のポアソン分布のパラメータ $\boldsymbol{\lambda} = \{\lambda_1, \cdots, \lambda_K\}$ を数値ベクトルで指定します。各値は、単位時間における発生回数の期待値であり、正の値 $\lambda_k \gt 0$ を満たす必要があります。

　クラスタの混合比率 $\boldsymbol{\pi}$ を設定する。

# 混合比率を指定
pi_k <- c(0.2, 0.3, 0.1, 0.4)
pi_k; sum(pi_k)

[1] 0.2 0.3 0.1 0.4
[1] 1

　カテゴリ分布のパラメータ $\boldsymbol{\pi} = (\pi_1, \cdots, \pi_K)$ を数値ベクトルで指定します。各値は、クラスタの割り当て確率であり、0から1の値 $0 \leq \pi_k \leq 1$ で、総和が1になる値 $\sum_{k=1}^K \pi_k = 1$ を満たす必要があります。

　クラスタ数 $K$ を設定します。

# クラスタ数を取得
K <- length(lambda_k)
K

[1] 4

　クラスタ数 $K$ を作成します。
　この例では、パラメータの要素数を使います。

　続いて、設定した値に従う乱数を生成し、グラフを作成していきます。

乱数の生成

　サンプルサイズ $N$ を指定します。

# サンプルサイズを指定
N <- 100

　サンプルサイズ(データ数) $N$ を指定します。

　カテゴリ分布の乱数 $\mathbf{s}_n$ を生成します。

# クラスタを生成
s_nk <- rmultinom(n =  N, size = 1, prob = pi_k) |> 
  t()
head(s_nk)

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    1
[2,]    1    0    0    0
[3,]    0    0    1    0
[4,]    0    1    0    0
[5,]    0    0    1    0
[6,]    0    0    0    1

　 $N$ 個のデータに対応するクラスタのサンプル(乱数)の集合 $\mathbf{S} = \{\mathbf{s}_1, \cdots, \mathbf{s}_N\}$ を作成します。各要素は、one-hotベクトル $\mathbf{s}_n = (s_{n,1}, \cdots, s_{n,K})$ であり、 $s_{n,k} = 1$ のときクラスタ $k$ を表します。
　カテゴリ分布の乱数は、多項分布の乱数生成関数 rmultinom() の試行回数の引数 size に 1 を指定することで生成できます。サンプルサイズの引数 n に $N$ 、割り当て確率(パラメータ)の引数 prob に $\boldsymbol{\pi}$ を指定します。
　出力されるマトリクスを t() で転置して、行をサンプル $n = 1, \dots, N$ 、列をクラスタ $k = 1, \dots, K$ に対応させ、 $\mathbf{S}$ を $N \times K$ のマトリクスとして扱います。

　各データのクラスタ番号 $s_n$ を抽出します。

# クラスタ番号を抽出
s_n <- which(t(s_nk) == 1, arr.ind = TRUE) |> 
  _[, "row"]
head(s_n)

[1] 4 1 3 2 3 4

　 $N$ 個のデータに対応するクラスタ番号の集合 $\mathbf{s} = \{s_1, \cdots, s_N\}$ を作成します。各要素は、データに割り当てられたクラスタ番号 $s_n \in \{1, \dots, K\}$ を表します。
　各データのクラスタ(one-hotベクトル) $\mathbf{s}_n$ について、値が1 $s_{n,k} = 1$ であるクラスタ番号(列番号) $s_n = k$ を、which() を使って取り出します。

　混合ポアソン分布の乱数 $x_n$ を生成します。

# サンプルを生成
x_n <- rpois(n = N, lambda = lambda_k[s_n])
head(x_n)

[1] 4 1 3 2 3 4

　 $N$ 個のデータのサンプル(乱数)の集合 $\mathbf{x} = \{x_1, \cdots, x_N\}$ を作成します。各要素は、単位時間における発生回数(0以上の整数) $x_n \in \{0, 1, 2, \cdots\}$ です。
　ポアソン分布の乱数は、rpois() で生成できます。サンプルサイズの引数 n に $N$ 、パラメータの引数 lambda にクラスタ $s_{n,k} = 1$ に対応するパラメータ $\lambda_{s_n} = \lambda_k$ を指定します。
　 $K$ 個のパラメータ $\boldsymbol{\lambda}$ に対して $N$ 個のクラスタ番号 $\mathbf{s}$ を添字として用いて、 $N$ 個のパラメータ $\lambda_{s_1}, \cdots, \lambda_{s_N}$ を取り出して、引数に指定できます。

# 割り当てクラスタのパラメータを抽出
lambda_n <- lambda_k[s_n]
head(lambda_n)

[1] 4 1 3 2 3 4

　 $N$ 個のクラスタ(one-hotベクトル) $\mathbf{s}_n$ を $\boldsymbol{\lambda}$ の指数として用いても、対応するパラメータを取り出せます。

# 割り当てクラスタのパラメータを抽出
lambda_n <- lambda_k^t(s_nk) |> 
  t() |> 
  apply(MARGIN = 1, FUN = prod)
head(lambda_n)

[1] 16.8  1.0 10.0  5.5 10.0 16.8

　この処理は、次の計算に対応します。

$\displaystyle \begin{aligned} \prod_{k'=1}^K \lambda_{k'}^{s_{n,k'}} &= \lambda_1^{s_{n,1}} * \cdots * \lambda_k^{s_{n,k'}} * \cdots * \lambda_K^{s_{n,K}} \\ &= \lambda_1^0 * \cdots * \lambda_k^1 * \cdots * \lambda_K^0 \\ &= 1 * \cdots * \lambda_k^1 * \cdots * 1 \\ &= \lambda_k \end{aligned}$

　0乗は1 $x^0 = 1$ なので、 $s_{n,k} = 1$ 以外の項が1になり $k$ 番目のパラメータ $\lambda_k$ のみが残ります。

　以上で、乱数を生成できました。続いて、乱数を可視化していきます。

乱数の集計

　サンプル集合 $\mathbf{x}$ をクラスタごとに集計します。

# x軸の範囲を指定
u <- 5
x_max <- x_n |> 
  max() |> 
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ

# クラスタごとにサンプルを集計
cluster_freq_df <- tibble::tibble(
  s = s_n, # クラスタ番号
  x = x_n  # サンプル値
) |> 
  dplyr::count(
    s, x, name = "freq" # 度数
  ) |> 
  dplyr::mutate(
    rel_freq = freq / N # 相対度数
  ) |> 
  tidyr::complete(
    s = 1:K, x = 0:x_max, 
    fill = list(freq = 0, rel_freq = 0)
  ) # 未観測値を補完
cluster_freq_df

# A tibble: 124 × 4
       s     x  freq rel_freq
   <int> <int> <int>    <dbl>
 1     1     0     7     0.07
 2     1     1     3     0.03
 3     1     2     3     0.03
 4     1     3     1     0.01
 5     1     4     1     0.01
 6     1     5     0     0   
 7     1     6     0     0   
 8     1     7     0     0   
 9     1     8     0     0   
10     1     9     0     0   
# ℹ 114 more rows

　サンプルのクラスタと値をデータフレームに格納して、count() で度数を求めます(重複する値をカウントします)。度数を $N_x$ ( 値が $1$ の度数を $N_1$ で表す)として、サンプルサイズ $N$ で割ると相対度数 $\frac{N_x}{N}$ が求まります。
　観測されなかったクラスタや値もデータフレームに含める場合は、complete() を使って補完します。欠損値となる他の列の値を fill 引数に指定できます。

　サンプル集合 $\mathbf{x}$ をクラスタをまとめて集計します。

# クラスタをまとめてサンプルを集計
freq_df <- cluster_freq_df |> 
  dplyr::summarise(
    freq     = sum(freq),     # 度数
    rel_freq = sum(rel_freq), # 相対度数
    .by = x
  ) # 全てのクラスタを合算
freq_df

# A tibble: 31 × 3
       x  freq rel_freq
   <int> <int>    <dbl>
 1     0     8     0.08
 2     1     3     0.03
 3     2     4     0.04
 4     3     8     0.08
 5     4     6     0.06
 6     5     8     0.08
 7     6     9     0.09
 8     7     1     0.01
 9     8     3     0.03
10     9     2     0.02
# ℹ 21 more rows

　各クラスタの度数と相対度数をそれぞれサンプルの値ごとに合計します。
　クラスタごとの集計を用いない場合は、次のように処理します。

# クラスタをまとめてサンプルを集計
freq_df <- tibble::tibble(
  x = x_n  # サンプル値
) |> 
  dplyr::count(
    x, name = "freq" # 度数
  ) |> 
  dplyr::mutate(
    rel_freq = freq / N # 相対度数
  ) |> 
  tidyr::complete(
    x = 0:x_max, 
    fill = list(freq = 0, rel_freq = 0)
  ) # 未観測値を補完
freq_df

# A tibble: 31 × 3
       x  freq rel_freq
   <int> <int>    <dbl>
 1     0     8     0.08
 2     1     3     0.03
 3     2     4     0.04
 4     3     8     0.08
 5     4     6     0.06
 6     5     8     0.08
 7     6     9     0.09
 8     7     1     0.01
 9     8     3     0.03
10     9     2     0.02
# ℹ 21 more rows

　クラスタ番号を含めずに同様に集計します。

　サンプルのクラスタ $s_n$ を s 列、サンプルの値 $x_n$ を x 列、度数 $N_x$ を freq 列、相対度数 $\frac{N_x}{N}$ を rel_freq 列として、データフレームに格納します。

乱数の作図

　混合ポアソン分布の乱数のヒストグラムを作成します。

# ラベル用の文字列を作成
param_lbl <- paste0(
  "list(", 
  "N == ", N, ", ", 
  "K == ", K, ", ", 
  "lambda == (list(", paste0(round(lambda_k, digits = 2), collapse = ", "), ")), ", 
  "pi == (list(", paste0(round(pi_k, digits = 2), collapse = ", "), "))", 
  ")"
) |> 
  parse(text = _)

# サンプルの度数を作図
ggplot() + 
  geom_bar(
    data = freq_df, 
    mapping = aes(x = x, y = freq), 
    stat = "identity", position = "identity", width = 1
  ) + # 度数
  labs(
    title = "mixture Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", 
    x = expression(x), 
    y = "frequency"
  )

　棒グラフの描画関数 geom_bar() の縦軸の引数 y に度数を指定してヒストグラムを作成します。

　クラスタごとに色分けした積み上げ棒グラフとしてヒストグラムを作成します。

# ラベル用の文字列を作成
cluster_lbl <- paste0("k == ", 1:K) |> 
  parse(text = _)
freq_lbl <- paste0(
  "N == (list(", paste0(colSums(s_nk), collapse = ", "), "))"
) |> 
  parse(text = _)

# サンプルの度数を作図
ggplot() + 
  geom_bar(
    data = cluster_freq_df, 
    mapping = aes(x = x, y = freq, fill = factor(s)), 
    stat = "identity", position = "stack", width = 1
  ) + # 度数
  geom_label(
    data = tibble::tibble(freq_lbl = freq_lbl), 
    mapping = aes(x = -Inf, y = Inf, label = freq_lbl),
    parse = TRUE, hjust = 0, vjust = 1, alpha = 0.5
  ) + # 割り当て数のラベル
  scale_fill_hue(label = cluster_lbl) + # クラスタ番号
  labs(
    title = "mixture Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", 
    x = expression(x), 
    y = "frequency"
  )

　y 引数にクラスタごとの度数を指定して、position 引数に "stack" を指定(デフォルトの設定)します。クラスタごとに色付けする場合は、fill (や color )引数に因子型に変換したクラスタ番号を指定します。

　混合ポアソン分布に従う確率を計算します。

作図コード(クリックで展開)

# x軸の値を作成
x_vec <- seq(from = 0, to = x_max, by = 1)

# クラスタごとの重み付けポアソン分布を計算
cluster_prob_df <- tidyr::expand_grid(
  k = 1:K,  # クラスタ番号
  x = x_vec # 確率変数
) |> # クラスタごとに変数を複製
  dplyr::mutate(
    lambda = lambda_k[k], # パラメータ
    pi     = pi_k[k],     # 混合比率
    prob   = pi * dpois(x = x, lambda = lambda) # 重み付け確率
  )

# ポアソン分布を計算
prob_df <- cluster_prob_df |> 
  dplyr::summarise(
    prob = sum(prob), # 周辺確率
    .by = x
  )

　詳しくは「混合ポアソン分布の作図」を参照してください。

　相対度数のグラフを作成します。

# サンプルの相対度数を作図
ggplot() + 
  geom_bar(
    data = prob_df, 
    mapping = aes(x = x, y = prob), 
    stat = "identity", position = "identity", 
    fill = NA, color = "darkgreen", linetype = "dashed"
  ) + # 生成確率
  geom_bar(
    data = freq_df, 
    mapping = aes(x = x, y = rel_freq), 
    stat = "identity", position = "identity", 
    fill = "#00A968", alpha = 0.5
  ) + # 相対度数
  scale_y_continuous(
    sec.axis = sec_axis(transform = ~ .*N, name = "freqency")
  ) + # 度数軸目盛
  labs(
    title = "mixture Poisson distribution", 
    subtitle = param_lbl, 
    x = expression(x), 
    y = "relative frequency"
  )

　形状の比較用に、混合ポアソン分布に従う確率(サンプルの生成分布)を破線で示します。

　y 引数に相対度数(サンプル全体に対する割合)を指定します。

　クラスタごとの相対度数をグラフで確認します。

# サンプルの相対度数を作図
ggplot() + 
  geom_bar(
    data = cluster_freq_df, 
    mapping = aes(x = x, y = rel_freq, fill = factor(s)), 
    stat = "identity", position = "identity", 
    alpha = 0.5
  ) + # 相対度数
  geom_bar(
    data = cluster_prob_df, 
    mapping = aes(x = x, y = prob, color = factor(k)), 
    stat = "identity", position = "identity", 
    fill = NA, linetype = "dashed"
  ) + # 生成分布
  geom_line(
    data = cluster_prob_df, 
    mapping = aes(x = x, y = prob, color = factor(k)), 
    linewidth = 1, linetype = "dashed"
  ) + # 生成分布
  geom_point(
    data = cluster_prob_df, 
    mapping = aes(x = x, y = prob, color = factor(k)), 
    size = 3
  ) + # 生成分布
  geom_label(
    data = tibble::tibble(freq_lbl = freq_lbl), 
    mapping = aes(x = -Inf, y = Inf, label = freq_lbl),
    parse = TRUE, hjust = 0, vjust = 1, alpha = 0.5
  ) + # 割り当て数のラベル
  scale_y_continuous(
    sec.axis = sec_axis(transform = ~ .*N, name = "freqency")
    ) + # 度数軸目盛
  scale_fill_hue(label = cluster_lbl) + # クラスタ番号
  scale_color_hue(label = cluster_lbl) + # クラスタ番号
  labs(
    title = "mixture Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", color = "cluster", 
    x = expression(x), 
    y = "relative frequency"
  )