はじめに

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

　この記事では、混合ポアソン分布の確率計算についてPythonを使って確認します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
混合ポアソン分布の計算
- 確率の計算
参考文献
おわりに

混合ポアソン分布の計算

　混合ポアソン分布(mixture Poisson distribution)に関する計算を実装します。この記事では、Pythonを利用して計算します。
　ポアソン分布については「【Python】ポアソン分布の計算 - からっぽのしょこ」、混合ポアソン分布については「混合ポアソン分布の定義式の導出 - からっぽのしょこ」、Rを利用する場合は「【R】混合ポアソン分布の計算 - からっぽのしょこ」を参照してください。

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

# ライブラリを読込
import numpy as np
import scipy as sp

　また、スクラッチ実装で利用するクラスや関数を読み込みます。

# クラス・関数を読込
from scipy.stats import poisson # ポアソン分布
from scipy.special import gamma # ガンマ関数

　この記事では(コードがごちゃごちゃしないように)、SciPyライブラリからクラスなどを個別に読み込んで利用します。

確率の計算

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

パラメータの設定

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

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

[ 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.35, 0.25, 0.4])
print(pi_k)
print(np.sum(pi_k))

[0.35 0.25 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 = 10.0
print(x)

10.0

　単位時間における発生回数(非負の整数) $x \in \{0, 1, 2, \cdots\}$ を指定します。

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

スクラッチ実装による計算

　定義式を実装して、確率を計算します。

## 定義式により確率を計算

# クラスタごとのポアソン分布の確率を計算
cluster_prob_k = np.exp(-lambda_k) * lambda_k**x / gamma(x + 1)
print(cluster_prob_k.round(3))

# クラスタごとの重み付け確率を計算
weighted_prob_k = pi_k * cluster_prob_k
print(weighted_prob_k.round(3))

# 周辺確率を計算
prob = np.sum(weighted_prob_k)
print(prob)

[0.029 0.125 0.025]
[0.01  0.031 0.01 ]
0.051244258225012516

　混合ポアソン分布(の周辺確率)は、次の式で定義されます。

$\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) \\ &= \sum_{k=1}^K \pi_k \exp(-\lambda_k) \frac{\lambda_k^x}{x!} \end{aligned}$

　階乗 $x!$ の計算は、ガンマ関数 $\Gamma(x + 1) = x!$ に置き換えて行います。ガンマ関数 $\Gamma(x)$ は、sp.special.gamma() で計算できます。
　ネイピア数 $e$ の指数関数 $e^x = \exp(x)$ は、np.exp() で計算できます。

　上の例では、 $K$ 個のクラスタの計算を配列を用いて処理しています。
　次の例では、クラスタ $k$ に関する総和 $\sum_{k=1}^K$ の計算をfor文で処理します。

## 定義式により確率を計算

# 周辺確率を初期化
prob = 0.0

# クラスタごとに処理
for k in range(K):
    
    # クラスタkのパラメータを取得
    lmd = lambda_k[k]
    pi  = pi_k[k]
    
    # クラスタkのポアソン分布の確率を計算
    cluster_prob = np.exp(-lmd) * lmd**x / gamma(x + 1)
    
    # クラスタkの重み付け確率を計算
    weighted_prob = pi * cluster_prob
    
    # 周辺確率を計算
    prob += weighted_prob
print(prob)

0.051244258225012516

　パラメータの配列から順番に値(要素)を取り出して、重み付け確率を計算し、prob に足していきます。始めに prob の値を 0 に初期化(値が0のオブジェクトを作成)しておく必要があります。

クラスによる計算

　scipy のクラスメソッドを使って、確率を計算します。

# クラスにより確率を計算
prob = np.sum(pi_k * poisson.pmf(k=x, mu=lambda_k))
print(prob)

0.05124425822501263

　クラスタとごとの確率計算と、周辺確率(重み付け確率の和)の計算を一度に行います。
　ポアソン分布の確率は、sp.stats.poisson.pmf() で計算できます。確率変数の引数 k に $x$ 、パラメータの引数 mu に $\boldsymbol{\lambda}$ を指定します。
　x 引数にスカラ、mu 引数に配列を指定すると、複数パラメータの計算を一度に処理できます。

　for文を使って、確率を計算します。

## クラスにより確率を計算

# 周辺確率を初期化
prob = 0.0

# クラスタごとに処理
for k in range(K):
    
    # クラスタkのパラメータを取得
    lmd = lambda_k[k]
    pi  = pi_k[k]
    
    # 周辺確率を計算
    prob += pi * poisson.pmf(k=x, mu=lmd)
print(prob)

0.05124425822501263

　「スクラッチ実装による計算」のクラスタごとの確率の計算を、確率メソッド pmf() に置き換えて処理します。

　リスト内包表記を使って、各種の確率を計算します。

# クラスタごとのポアソン分布の確率を計算
cluster_prob_k = np.array(
    [poisson.pmf(k=x, mu=lambda_k[k]) for k in range(K)]
)
print(cluster_prob_k.round(3))

# クラスタごとの重み付け確率を計算
weighted_prob_k = np.array(
    [pi_k[k] * poisson.pmf(k=x, mu=lambda_k[k]) for k in range(K)]
)
print(weighted_prob_k.round(3))

# 周辺確率を計算
prob = np.sum(
    [pi_k[k] * poisson.pmf(k=x, mu=lambda_k[k]) for k in range(K)]
)
print(prob)

[0.029 0.125 0.025]
[0.01  0.031 0.01 ]
0.05124425822501263

　クラスタ( lambda_k, pi_k の要素)ごとに、確率や重み付け確率を計算して、リストに格納します。リスト内包表記により、各パラメータ $\lambda_k$ に対して計算できます。

スポンサードリンク

複数データの場合

　続いて、作図などでの利用を想定して、変数が複数の場合の処理を確認します。

　複数個の確率変数の実現値 $\mathbf{x}$ を設定します。

# 確率変数の値を指定
x_n = np.array([10.0, 6.0, 17.0, 13.0, 15.0])
print(x_n)

# データ数を設定
N = len(x_n)
print(N)

[10.  6. 17. 13. 15.]
5

　 $N$ 個の変数 $\mathbf{x} = \{x_1, \cdots, x_N\}$ を数値ベクトルで指定します。
　計算の処理用に、データ数 $N$ を作成します。

　途中計算用に、変数 $\mathbf{x}$ とパラメータ $\boldsymbol{\lambda}$ の値を複製します。

# 確率変数を複製
x_nk = np.tile(x_n, reps=(K, 1)).T
print(x_nk.shape)
print(x_nk)

# パラメータを複製
lambda_nk = np.tile(lambda_k, reps=(N, 1))
print(lambda_nk.shape)
print(lambda_nk)

(5, 3)
[[10. 10. 10.]
 [ 6.  6.  6.]
 [17. 17. 17.]
 [13. 13. 13.]
 [15. 15. 15.]]
(5, 3)
[[ 5.5 10.  16.8]
 [ 5.5 10.  16.8]
 [ 5.5 10.  16.8]
 [ 5.5 10.  16.8]
 [ 5.5 10.  16.8]]

　 $N$ 個の変数 $n = 1, \dots, N$ が行、 $K$ 個のクラスタ $k = 1, \dots, K$ が列に対応するように、 $\mathbf{x}, \boldsymbol{\lambda}$ の1次元配列をそれぞれ複製して、 $N$ 行 $K$ 列の2次元配列を作成します。
　または、np.meshgrid() で複製します。

# 確率変数・パラメータを複製
lambda_nk, x_nk = np.meshgrid(lambda_k, x_n)
print(x_nk.shape)
print(x_nk)
print(lambda_nk.shape)
print(lambda_nk)

(5, 3)
[[10. 10. 10.]
 [ 6.  6.  6.]
 [17. 17. 17.]
 [13. 13. 13.]
 [15. 15. 15.]]
(5, 3)
[[ 5.5 10.  16.8]
 [ 5.5 10.  16.8]
 [ 5.5 10.  16.8]
 [ 5.5 10.  16.8]
 [ 5.5 10.  16.8]]

　配列の渡し方によって、出力される配列の形状が変わります。

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

# クラスタごとのポアソン分布の確率を計算
cluster_prob_nk = poisson.pmf(k=x_nk, mu=lambda_nk)
print(cluster_prob_nk.round(3))

# クラスタごとの重み付け確率を計算
weighted_prob_nk = pi_k * cluster_prob_nk
print(weighted_prob_nk.round(3))

# 周辺確率を計算
prob_n = np.sum(weighted_prob_nk, axis=1)
print(prob_n.round(3))

[[0.029 0.125 0.025]
 [0.157 0.063 0.002]
 [0.    0.013 0.096]
 [0.003 0.073 0.069]
 [0.    0.035 0.093]]
[[0.01  0.031 0.01 ]
 [0.055 0.016 0.001]
 [0.    0.003 0.038]
 [0.001 0.018 0.028]
 [0.    0.009 0.037]]
[0.051 0.071 0.042 0.047 0.046]

　「関数による計算」では、sp.stats.poisson.pmf() に変数 $x$ のスカラとパラメータ $\boldsymbol{\lambda}$ の配列を指定して、複数パラメータの計算を行いました。
　今回は、それぞれ複製して形状を整えた配列の変数 $\mathbf{x}$ とパラメータ $\boldsymbol{\lambda}$ を指定して、各変数 $x_n$ と各パラメータ $\lambda_k$ の全ての組み合わせに対して計算します。
　pmf() の各引数に同じ形状の配列を指定すると、同じ位置の要素ごとに確率を計算できます。
　または、重み付け和の計算を、np.dot() によりベクトル演算として処理します。

# 周辺確率を計算
prob_n = np.dot(cluster_prob_nk, pi_k)
print(prob_n.round(3))

[0.051 0.071 0.042 0.047 0.046]

　配列の渡し方によって、計算結果が変わることや、計算できないことがあります。この計算は、次のようにも行えます。

# 周辺確率を計算
prob_n = np.dot(pi_k, cluster_prob_nk.T)
print(prob_n.round(3))

[0.051 0.071 0.042 0.047 0.046]

　ここからは、中間オブジェクトを作らずにそれぞれの計算を行います。

　for文のリスト内包表記を使って、各種の確率を計算します。

# クラスタごとのポアソン分布の確率を計算
cluster_prob_nk = np.array(
    [poisson.pmf(k=x, mu=lambda_k) for x in x_n]
)
print(cluster_prob_nk.round(3))

# クラスタごとの重み付け確率を計算
weighted_prob_nk = np.array(
    [pi_k * poisson.pmf(k=x, mu=lambda_k) for x in x_n]
)
print(weighted_prob_nk.round(3))

# 周辺確率を計算
prob_n = np.array(
    [np.sum(pi_k * poisson.pmf(k=x, mu=lambda_k)) for x in x_n]
)
print(prob_n.round(3))

[[0.029 0.125 0.025]
 [0.157 0.063 0.002]
 [0.    0.013 0.096]
 [0.003 0.073 0.069]
 [0.    0.035 0.093]]
[[0.01  0.031 0.01 ]
 [0.055 0.016 0.001]
 [0.    0.003 0.038]
 [0.001 0.018 0.028]
 [0.    0.009 0.037]]
[0.051 0.071 0.042 0.047 0.046]

　クラスタ( lambda_k, pi_k の要素)ごとに、確率や重み付け確率を計算して、リストに格納します。リスト内包表記により、各変数 $x_n$ と各パラメータ $\lambda_k$ の全ての組み合わせに対して計算できます。

　NumPyのブロードキャストを使って、各種の確率を計算します。

# クラスタごとのポアソン分布の確率を計算
cluster_prob_nk = poisson.pmf(k=x_n[:, np.newaxis], mu=lambda_k)
print(cluster_prob_nk.round(3))

# クラスタごとの重み付け確率を計算
weighted_prob_nk = pi_k * poisson.pmf(k=x_n[:, np.newaxis], mu=lambda_k)
print(weighted_prob_nk.round(3))

# 周辺確率を計算
prob_n = np.sum(
    pi_k * poisson.pmf(k=x_n[:, np.newaxis], mu=lambda_k), 
    axis=1
)
print(prob_n.round(3))

[[0.029 0.125 0.025]
 [0.157 0.063 0.002]
 [0.    0.013 0.096]
 [0.003 0.073 0.069]
 [0.    0.035 0.093]]
[[0.01  0.031 0.01 ]
 [0.055 0.016 0.001]
 [0.    0.003 0.038]
 [0.001 0.018 0.028]
 [0.    0.009 0.037]]
[0.051 0.071 0.042 0.047 0.046]

　要素数 $N$ の1次元配列 x_n を x_n[:, np.newaxis] で $N$ 行 $1$ 列の2次元配列に変換して、要素数 $K$ の1次元配列 lambda_k と計算すると、x_n の各要素が各行の全ての要素に、lambda_k の各要素が各列の全ての要素に、それぞれブロードキャストされ、計算結果が $N$ 行 $K$ 列の2次元配列になります。ブロードキャストにより、各変数 $x_n$ と各パラメータ $\lambda_k$ の全ての組み合わせに対して計算できます。