はじめに

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

　この記事では、混合ポアソン分布のグラフを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 を読み込む必要があります。

グラフの作成

　ggplot2 パッケージを使って、混合ポアソン分布のグラフを作成します。

パラメータの設定

　混合ポアソン分布のパラメータ $\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$ を作成します。
　この例では、パラメータの要素数を使います。

　確率変数がとり得る値 $x$ を設定します。

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

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

[1] 0 1 2 3 4 5

　単位時間における発生回数(0以上の整数) $x \in \{0, 1, 2, \cdots\}$ の数値ベクトルを作成します。
　この例では、指定したパラメータの最大値を使って範囲を設定します。

　続いて、設定した値に従う確率を計算し、グラフを作成していきます。

確率分布の計算

　 $x$ の値ごとに混合ポアソン分布に従う重み付け確率 $p(x, s = k \mid \lambda_k, \boldsymbol{\pi})$ を計算します。

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

# A tibble: 124 × 5
       k     x lambda    pi weighted_prob
   <int> <dbl>  <dbl> <dbl>         <dbl>
 1     1     0      1   0.2   0.0736     
 2     1     1      1   0.2   0.0736     
 3     1     2      1   0.2   0.0368     
 4     1     3      1   0.2   0.0123     
 5     1     4      1   0.2   0.00307    
 6     1     5      1   0.2   0.000613   
 7     1     6      1   0.2   0.000102   
 8     1     7      1   0.2   0.0000146  
 9     1     8      1   0.2   0.00000182 
10     1     9      1   0.2   0.000000203
# ℹ 114 more rows

　パラメータ $\boldsymbol{\lambda}, \boldsymbol{\pi}$ と変数 $x$ の全ての組み合わせを expand_grid() で作成して、組み合わせ(行)ごとに重み付け確率を計算します。
　ポアソン分布の確率は、dpois() で計算できます。確率変数の引数 x に $x$ の列、パラメータの引数 lambda に $\lambda_k$ の列を指定します。

　混合ポアソン分布の重み付け確率(クラスタの混合比率で割り引いた確率)は、次の式で計算できます。

$\displaystyle \begin{aligned} p(x, s = k \mid \lambda_k, \boldsymbol{\pi}) &= p(s = k \mid \boldsymbol{\pi}) p(x \mid s = k, \lambda_k) \\ &= \pi_k \mathrm{Poisson}(x \mid \lambda_k) \end{aligned}$

　クラスタ番号 $k$ を k 列、確率変数 $x$ を x 列、パラメータ $\lambda_k$ を lambda 列、混合比率 $\pi_k$ を pi 列、重み付け確率を weighted_prob 列として、データフレームに格納します。

　 $x$ の値ごとに混合ポアソン分布に従う周辺確率 $p(x \mid \boldsymbol{\lambda}, \boldsymbol{\pi})$ を計算します。

# 周辺確率を計算
prob_df <- cluster_prob_df |> 
  dplyr::summarise(
    marginal_prob = sum(weighted_prob), # 周辺確率
    .by = x
  )
prob_df

# A tibble: 31 × 2
       x marginal_prob
   <dbl>         <dbl>
 1     0        0.0748
 2     1        0.0804
 3     2        0.0556
 4     3        0.0470
 5     4        0.0518
 6     5        0.0560
 7     6        0.0542
 8     7        0.0476
 9     8        0.0399
10     9        0.0340
# ℹ 21 more rows

　変数 $x$ ごとに重み付け確率の和を計算します。

　混合ポアソン分布の周辺確率(全てのクラスタを考慮した生成確率)は、次の式で計算できます。

$\displaystyle \begin{aligned} p(x \mid \boldsymbol{\lambda}, \boldsymbol{\pi}) &= \sum_{k=1}^K p(x, s = k \mid \lambda_k, \boldsymbol{\pi}) \\ &= \sum_{k=1}^K p(s = k \mid \boldsymbol{\pi}) p(x \mid s = k, \lambda_k) \\ &= \sum_{k=1}^K \pi_k \mathrm{Poisson}(x \mid \lambda_k) \end{aligned}$

　確率変数 $x$ を x 列、周辺確率を marginal_prob 列として、データフレームに格納します。

確率分布の作図

　パラメータラベルを作成します。装飾用の処理です。

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

expression(list(K == 4, lambda == (list(1, 5.5, 10, 16.8)), pi == 
    (list(0.2, 0.3, 0.1, 0.4))))

　ギリシャ文字などの記号を使った数式を表示する場合は、expression() の記法を使います。等号は "=="、複数の(数式上の)変数は list(変数1, 変数2) で並べて描画できます。
　オブジェクト(プログラム上の変数)の値を使う場合は、文字列を作成して parse() の text 引数に渡します。

　混合ポアソン分布のグラフを作成します。

# 混合ポアソン分布を作図
ggplot() + 
  geom_bar(
    data = prob_df, 
    mapping = aes(x = x, y = marginal_prob), 
    stat = "identity", position = "identity", 
    width = 1
  ) + # 周辺確率
  labs(
    title = "mixture Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", color = "cluster", 
    x = expression(x), 
    y = "probability"
  )

　棒グラフの描画関数 geom_bar() の縦軸の引数 y に周辺確率を指定します。

クラスタと分布の関係

　混合ポアソン分布のグラフを作成します。

# ラベル用の文字列を作成
cluster_lbl <- paste0("k == ", 1:K) |> 
  parse(text = _)
prob_lbl <- expression(
  sum({}, k == 1, K) * p(x, s == k ~'|'~ lambda[k], pi) == sum({}, k == 1, K) * pi[k] * p(x ~'|'~ s == k, lambda[k])
)

# 混合ポアソン分布を作図
ggplot() + 
  geom_bar(
    data = cluster_prob_df, 
    mapping = aes(x = x, y = weighted_prob, fill = factor(k)), 
    stat = "identity", position = "stack"
  ) + # 周辺確率
  scale_fill_hue(label = cluster_lbl) + # クラスタ番号
  labs(
    title = "mixture Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", 
    x = expression(x), 
    y = prob_lbl
  )

　縦軸の引数 y に重み付け確率、バーの並べ方の引数 position に積み上げ "stack" を指定します。
　周辺確率(の列)を用意しなくても、重み付け確率の積み上げで周辺確率を描画できます。

　クラスタごとに色付けして、各クラスタの重み付け確率を示しています。
　バーの積み上げが、クラスタの周辺化(重み付け和の計算) $\sum_{k=1}^K p(x, s = k \mid \lambda_k, \boldsymbol{\pi}) = \sum_{k=1}^K \pi_k \mathrm{Poi}(x \mid \lambda_k)$ に対応しています。
　混合比率などの影響については「【R】混合ポアソン分布のパラメータの可視化 - からっぽのしょこ」を参照してください。

　以上で、作図処理を確認しました。

　この記事では、混合ポアソン分布のグラフを作成しました。次の記事では、パラメータの影響を可視化します。

参考文献

機械学習スタートアップシリーズ　ベイズ推論による機械学習入門

作者:須山敦志
講談社

Amazon

おわりに

　初稿などこれまでは「確率分布の作図」の中でパラメータの影響も扱っていたのですが、今回の改修からはタイトル通りの作図処理の解説に留め、パラメータの影響については「確率分布のパラメータの可視化」として別記事に独立させることにしました。
　1記事1テーマに絞った方が執筆作業的には扱いやすくなった気がしますが、記事の管理的には面倒臭くなりそうな気がしてきました。書くときに書いた後のことなんか考えないんだよなー。

　最後に、前回の記事でほんの少し触れたJuice=Juiceの新メンバーの林仁愛さんのソロ楽曲をどうぞ♪

　デビュー前でこの歌唱力は凄いなぁ。

【次の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【R】混合ポアソン分布の作図