はじめに

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

　この記事では、混合ポアソン分布のグラフ作成ついてR言語を使って確認します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
混合ポアソン分布の作図
- グラフの作成
参考文献
おわりに

混合ポアソン分布の作図

　混合ポアソン分布(mixture Poisson distribution)のグラフを作成します。この記事では、PythonのMatplotlibライブラリを利用して作図します。
　ポアソン分布については「【Python】ポアソン分布の作図 - からっぽのしょこ」、混合ポアソン分布については「混合ポアソン分布の定義式の導出 - からっぽのしょこ」、確率計算については「【Python】混合ポアソン分布の計算 - からっぽのしょこ」、Rを利用する場合は「【R】混合ポアソン分布の作図 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import numpy as np
from scipy.stats import poisson
import matplotlib.pyplot as plt

グラフの作成

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

パラメータの設定

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

# 期待値パラメータを指定
lambda_k = np.array([1.0, 5.5, 10.0, 16.8])
print(lambda_k)

[ 1.   5.5 10.  16.8]

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

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

# 混合比率パラメータを指定
pi_k = np.array([0.2, 0.3, 0.1, 0.4])
print(pi_k)
print(np.sum(pi_k))

[0.2 0.3 0.1 0.4]
1.0

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

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

# クラスタ数を設定
K = len(lambda_k)
print(K)

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

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

# x軸の範囲を設定
u = 5.0
x_max = np.max(lambda_k) # 基準値を指定
x_max *= 1.5 # 倍率を指定
x_max = np.ceil(x_max /u)*u # u単位で切り上げ
#x_max = 30
print(x_max)

# x軸の値を作成
x_vec = np.arange(start=0, stop=x_max+1, step=1)
print(x_vec[:5])

30.0
[0. 1. 2. 3. 4.]

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

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

確率分布の計算

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

# クラスタごとの重み付け確率を計算
weighted_prob_lt = [
    pi_k[k] * poisson.pmf(k=x_vec, mu=lambda_k[k]) for k in range(K)
]
print(len(weighted_prob_lt))
print(weighted_prob_lt[0][:5])

(4, 31)
[0.07357589 0.07357589 0.03678794 0.01226265 0.00306566]

　クラスタ( lambda_k, pi_k の要素)ごとに、重み付け確率を計算して、リストに格納します。リスト内包表記により、パラメータ $\boldsymbol{\lambda}, \boldsymbol{\pi}$ のペアと変数 $x$ の全ての組み合わせに対して計算できます。
　ポアソン分布の確率は、sp.stats.poisson.pmf() で計算できます。確率変数の引数 k に $x$ (のベクトル)、パラメータの引数 mu に $\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}$

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

# 周辺確率を計算
prob_vec = np.sum(weighted_prob_lt, axis=0)
print(prob_vec[:5])

(31,)
[0.07480648 0.0803648  0.05556152 0.04703213 0.0517701 ]

　変数( x_vec の要素)ごとに、重み付け確率の和を計算します。
　列ごとの和は、np.sum(axis=0) で計算できます。

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

$\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}$

確率分布の作図

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

# ラベル用の文字列を作成
lambda_str = ', '.join([f'{lmd}' for lmd in lambda_k])
pi_str     = ', '.join([f'{pi}' for pi in pi_k])
param_lbl = f'$K = {K}, \\lambda = ({lambda_str}), \\pi = ({pi_str})$'
print(param_lbl)

$K = 4, \lambda = (1.0, 5.5, 10.0, 16.8), \pi = (0.2, 0.3, 0.1, 0.4)$

　ギリシャ文字などの記号を使った数式を表示する場合は、LaTeXの記法(LaTeXコマンド)を使います。数式に変換する範囲(コマンドの文字列)を $ で挟んで、'$コマンド$' の形式で文字列を指定します。ただし、メタ文字(改行文字 \n など)として認識されないために、バックスラッシュを追加(2つ \\ に)する必要がある場合があります。
　オブジェクト(プログラム上の変数)の値を使う場合は、f文字列(や format() )を使います。変数名を波括弧 {} で挟んで、f'{変数名}' の形式で文字列を指定します。(丸め込みなど)数値の文字列を調整する場合は、format() の記法で指定します。
　リスト内の文字列を結合する場合は、join() を使います。'区切り文字列'.join(リスト) 形式で指定します。

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

# 混合ポアソン分布を作図:棒グラフ
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
plt.bar(
    x=x_vec, height=prob_vec, 
    color='#00A968'
) # 周辺確率
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('probability') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('mixture Poisson distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.show()

　棒グラフは、plt.bar() で描画できます。縦軸の引数(第2引数) height に確率を指定します。

　幹図を作成します。

# 混合ポアソン分布を作図:幹図
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
plt.vlines(
    x=x_vec, ymin=0.0, ymax=prob_vec, 
    color='black', linewidth=2.0, zorder=0
) # 線分
plt.scatter(
    x=x_vec, y=prob_vec, 
    color='#00A968', s=100, zorder=1
) # 点
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('probability') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('mixture Poisson distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.show()

　幹図は、plt.vlines() と plt.scatter() を使って描画できます。plt.vlines() の x 引数にx座標、ymin 引数に0、ymax 引数に確率値を指定して、線分を描画します。先端に重ねるように、plt.scatter() で確率値の位置に点を描画します。
　描画順(重ねる順番)は zorder 引数の値の大小で指定できます。

　折れ線グラフを作成します。

# 混合ポアソン分布を作図:折れ線グラフ
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
plt.plot(
    x_vec, prob_vec, 
    color='black', linewidth=2.0, zorder=0
) # 線分
plt.scatter(
    x=x_vec, y=prob_vec, 
    color='#00A968', s=100, zorder=1
) # 点
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('probability') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('mixture Poisson distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.show()

　折れ線グラフは、plt.plot() で描画できます。
　混合ポアソン分布は離散型の確率分布なので、変数がとり得る値(離散値)の位置に点を描画します。

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

# 混合ポアソン分布を作図:積み上げ棒グラフ
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
for k in range(K):
    weighted_prob_vec = weighted_prob_lt[k] # 重み付け確率
    bottom_vec = np.sum(weighted_prob_lt[:k], axis=0) # 累積重み付け確率
    plt.bar(
        x=x_vec, bottom=bottom_vec, height=weighted_prob_vec, 
        label=f'$k = {k+1}$'
    ) # 累積重み付け確率
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('probability') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('mixture Poisson distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.legend(title='cluster') # 凡例
plt.show()

　plt.bar() の bottom 引数にクラスタ $k$ までの重み付け確率の和、height 引数にクラスタ $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)$ に対応します。
　混合比率などの影響については「【R】混合ポアソン分布のパラメータの可視化 - からっぽのしょこ」を参照してください。

　以上で、基本的な作図処理を確認しました。

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

参考文献

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

作者:須山敦志
講談社

Amazon

おわりに

　LLMの登場以降、アクセス数の激減が留まるところを知らず、もうブログ更新のモチベが萎えなえです。手元のスクリプトやノートとしてはこれまで通り書けているのですが、はてなブログに投稿するための作業(数式を表示するために細々とした調整が必要なんですや図の貼付、補完する他の記事とのリンク張りなどなど)の手間を乗り越えられず、更新が滞っています。
　随分遅いかと思いますが先月頃から自分でも検索替わりに使ってみているのですが、まぁ便利ですねググることが減りました……。
　解説を書くことに関しては、自分の理解に突っ込みを入れながら書いたり消したりすること自体が私にとっての勉強(自分に合った勉強法)なので、(つまりブログは目的ではなく手段というか副産物なので)今後も手作業で書くのですが、もう時代についていける気がしません、なむ。

　そんなこんなで気が滅入るので、軽い運動を始めたところ(1年ほど軽い食事制限を主に続けてきたので狭義の意味での)ダイエットが捗っています。7年くらい続けてきた勉強習慣が危うい状況ですが、いい方向に進みだした健康的な食事習慣や運動習慣をこのまま続けて身に付けたいです。

　話がとっ散らかっていますが、そんな今日この頃です。

　2025年9月11日は、juice=Juiceのメジャーデビュー12周年の日です。

　12年かぁ、オリメンも全員卒業してグループとして一回りした感じはありますね。新メンバーもデビュー前にソロ曲を担当しているとかまた凄そうな人で、これからも楽しそうでなによりです。
　ところで、えびジュースの1次選考に落選しまして、今日2次選考を申し込みました。ずっと待ち望んでいたコラボライブなんです、どうか当選のほどよろしくお願いしますっ！！

【次の内容】

zenn.dev