はじめに

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

　この記事では、混合ポアソン分布のパラメータの影響について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 を読み込む必要があります。

周辺化の計算

　まずは、クラスタの周辺化(重み付け確率の和・クラスタごとの分布と重みの線形結合)の計算を図で確認します。

　パラメータを指定して、混合ポアソン分布に従う確率を計算します。

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

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

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

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

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

# 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],     # 混合比率
    cluster_prob  = dpois(x = x, lambda = lambda), # 確率
    weighted_prob = pi * cluster_prob # 重み付け確率
  )

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

　クラスタごとのポアソン分布 $\mathrm{Poi}(x \mid \lambda_k)$ をグラフで確認します。

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

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

# 確率軸の範囲を設定
u <- 0.005
cluster_prob_max <- cluster_prob_df |> 
  dplyr::pull(cluster_prob) |> 
  max() |> 
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ
marginal_prob_max <- cluster_prob_df |> 
  dplyr::summarise(
    marginal_prob = sum(weighted_prob), # 周辺確率
    .by =　x
  ) |> 
  dplyr::pull(marginal_prob) |> 
  max() |> 
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ
prob_max <- max(cluster_prob_max, marginal_prob_max)

# ポアソン分布を作図
ggplot() + 
  geom_bar(
    data = cluster_prob_df, 
    mapping = aes(x = x, y = cluster_prob, fill = factor(k)), 
    stat = "identity", position = "identity", 
    alpha = 0.5
  ) + # 周辺確率
  scale_fill_hue(label = cluster_lbl) + # クラスタ番号
  coord_cartesian(ylim = c(0, prob_max)) + # 確率軸の比較用
  labs(
    title = "Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", 
    x = expression(x), 
    y = prob_lbl
  )

　クラスタごとに色分けしたグラフを並べて各クラスタを示しています。
　それぞれの山が1つのポアソン分布であり、それぞれの総和が $\sum_{x=0}^{\infty} \mathrm{Poi}(x \mid \lambda_k) = 1$ です。

　クラスタごとに重み付けしたポアソン分布 $\pi_k \mathrm{Poi}(x \mid \lambda_k)$ をグラフで確認します。

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

# ラベル用の文字列を作成
prob_lbl <- expression(
  p(x, s == k ~'|'~ lambda[k], pi) == 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 = "identity", 
    alpha = 0.5
  ) + # 周辺確率
  scale_fill_hue(label = cluster_lbl) + # クラスタ番号
  coord_cartesian(ylim = c(0, prob_max)) + # 確率軸の比較用
  labs(
    title = "weighted Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", 
    x = expression(x), 
    y = prob_lbl
  )

　それぞれの山に $0 \leq \pi_k \leq 1$ 、 $\sum_{k=1}^K \pi_k = 1$ である重み(混合比率)を掛けることで、山全体(各バーの高さ)が $\pi_k$ 倍され、それぞれの総和が $\sum_{x=0}^{\infty} \pi_k \mathrm{Poi}(x \mid \lambda_k) = \pi_k$ となります。
　よって、 $K$ 個の山の総和が $\sum_{k=1}^K \sum_{x=0}^{\infty} \pi_k \mathrm{Poi}(x \mid \lambda_k) = 1$ になります。

　クラスタごとに重み付けしたポアソン分布の和 $\sum_{k=1}^K \pi_k \mathrm{Poi}(x \mid \lambda_k)$ をグラフで確認します。

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

# ラベル用の文字列を作成
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) + # クラスタ番号
  coord_cartesian(ylim = c(0, prob_max)) + # 確率軸の比較用
  labs(
    title = "mixture Poisson distribution", 
    subtitle = param_lbl, 
    fill = "cluster", 
    x = expression(x), 
    y = prob_lbl
  )

　ここまではクラスタごとに並べていたバーを、 $x$ の値ごとに積み上げて描画しています。
　 $K$ 個のバーの積み上げが、全てのクラスタの和(周辺化) $\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)$ に対応します。

　以上の操作(計算)により、 $x$ がとり得る値それぞれの確率が0から1の値で総和が1になることから、(離散型の)確率分布の条件を満たしました。

スポンサードリンク

パラメータの影響

　次は、パラメータの値を少しずつ変化させてグラフの変化をアニメーションで確認します。
　作図コードについては「Probability-Distribution/code/mixture_poisson/parameter.R at main · anemptyarchive/Probability-Distribution · GitHub」を参照してください。

パラメータと形状の関係

　パラメータ $\boldsymbol{\lambda}$ を変化させたときの混合ポアソン分布の形状の変化をアニメーションにします。

　クラスタごとの期待値を破線、期待値を中心に標準偏差1つ分の範囲を線分で示します。

　各クラスタのパラメータ(発生回数の期待値) $\lambda_k$ が大きくなるに従って、確率変数の値(発生回数) $x$ が大きいほど確率が高くなり(山が右に移動し)分散が大きくなる(山の裾が広がる)のが分かります。

混合比率と形状の関係

　混合比率 $\boldsymbol{\pi}$ を変化させたときの混合ポアソン分布の形状の変化をアニメーションにします。

　各クラスタの混合比率(重み) $\pi_k$ が大きくなるに従って、対応するクラスタのポアソン分布の影響が大きくなる(山が高くなる)のが分かります。ただし、総和が1になる必要があるので、他のクラスタの影響が小さく(山が低く)なります。

　以上で、パラメータの影響を確認しました。

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

参考文献

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

作者:須山敦志
講談社

Amazon

おわりに

　これまでは「確率分布の作図」の中でパラメータの影響を可視化するアニメーションを掲載して作図処理を解説していたのですが、よくよく考えるとアニメーションを見たい人はいても作りたい人は少なかろうと思い至り、「確率分布の作図」(前の記事)ではグラフ作成の解説に留めて、グラフ自体に注目するのは「パラメータの可視化」(この記事)ですることにしました。本質的ではない内容なのに書くのが面倒な上に需要もなさそうなアニメーション作成の解説は止めました。面倒な処理の解説の一番の読者は将来の自分なので、備忘録的な意味合いで書いていた解説を止めてしまって一番困るのは将来の自分なわけですが、まぁごめん恨まないで。
　これまでに書いた分布についても現在進めている修正時に分割していくつもりですが、パラメータが1つの分布なんかだと動画を1つ載せただけの記事になりそうで今から悩ましいです。混合分布だと混合パラメータがあるので複数個の図を確保できる上に周辺化の図も載せたので文量的に問題なかったのですが、同じシリーズの全ての記事で構成を統一するにはどのレベルを基準にするかも考えどころです。大抵の場合は簡単なものから書いていくので全部載せ記事を書いて、レベルが上がるとともに崩壊します。
　今回は摘まみ食い的に修正なり追記なり新規で書いたりしているのでどうなるか分かりませんが、前回の改修時にはPython版まで手が回らなかったっぽいので全部やり切りたいと思いつつ現時点ではPythonコードは一文字も書けていません。

　最後に、先日公開された宮本佳林さんのカバー曲をどうぞ♪

　普段のパフォーマンスの雰囲気的にも好きなことをパフォーマンスにするという意味でもボカロ曲をカバーをしてくれるのは嬉しいのですが、微妙に世代が違うので知らない曲なんだよなぁ。宮本さんの歌唱で知れてよかったとも言えるのですが、老いを感じてしまいます。
　どうも、ボカロ曲というとメルトが思い浮かぶ世代です。このブログの読者層の世代構成はどんな感じなのでしょうか。

【次の内容】

からっぽのしょこ

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

【R】混合ポアソン分布のパラメータの可視化

はじめに

混合ポアソン分布のパラメータの可視化

周辺化の計算

パラメータの影響

パラメータと形状の関係

混合比率と形状の関係

参考文献

おわりに