はじめに

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

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

【前の内容】

zenn.dev

【他の内容】

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$ を作成します。
　この例では、パラメータの要素数を使います。

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

乱数の生成

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

# サンプルサイズを指定
N = 300

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

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

# クラスタを生成
s_nk = np.random.multinomial(n=1, pvals=pi_k, size=N)
print(s_nk.shape)
print(s_nk[:5])

(300, 4)
[[0 0 0 1]
 [1 0 0 0]
 [0 0 0 1]
 [0 1 0 0]
 [0 1 0 0]]

　 $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$ を表します。
　カテゴリ分布の乱数は、多項分布の乱数生成関数 np.random.multinomial() の試行回数の引数 size に 1 を指定することで生成できます。割り当て確率の引数 pvals に $\boldsymbol{\pi}$ 、サンプルサイズの引数 size に $N$ を指定します。
　行がサンプル $n = 1, \dots, N$ 、列がクラスタ $k = 1, \dots, K$ に対応する $N \times K$ の2次元配列を $\mathbf{S}$ として扱います。

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

# クラスタ番号を抽出
_, s_n = np.where(s_nk == 1)
print(s_n[:5])

[3 0 3 1 1]

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

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

# サンプルを生成
x_n = np.random.poisson(lam=lambda_k[s_n], size=N)
print(x_n.shape)
print(x_n[:10])

(300,)
[17  2 16  6  5 18 23 16  1 18]

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

# 割り当てクラスタのパラメータを抽出
lambda_n = lambda_k[s_n]
print(lambda_n[:5])

[16.8  1.  16.8  5.5  5.5]

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

# 割り当てクラスタのパラメータを抽出
lambda_n = np.prod(lambda_k**s_nk, axis=1)
print(lambda_n[:5])

[16.8  1.  16.8  5.5  5.5]

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

$\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乗は $x^0 = 1$ なので、 $s_{n,k} = 1$ 以外の項が1になり、 $K$ 個の項を掛け合わせると $k$ 番目のパラメータ $\lambda_k$ のみが残ります。

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

乱数の集計

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

# x軸の範囲を設定
u = 5.0
x_max = np.max(x_n) # 基準値を指定
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\}$ の数値ベクトルを作成します。
　この例では、生成したサンプル x_n の最大値を使って、範囲を設定します。

　サンプル集合 $\mathbf{x}$ を集計します。

# 度数を集計
freq_vec = np.array([np.sum(x_n == x) for x in x_vec])
print(freq_vec[:10])

# 相対度数を計算
relfreq_vec = freq_vec / N
print(relfreq_vec[:10])

[20 23 26 12 17 18 14 16  7 10]
[0.06666667 0.07666667 0.08666667 0.04       0.05666667 0.06
 0.04666667 0.05333333 0.02333333 0.03333333]

　 $x$ がとり得る値( x_vec の要素)ごとに、サンプル x_n に含まれる要素数をカウントして(度数を求めて)、配列に格納します。True は 1、False は 0 として扱える(計算できる)ので、値が x の要素数は np.sum(x_n == x) で得られます。リスト内包表記により、全ての $x$ の値に対してカウントできます。
　度数を $N_x$ ( 値が $1$ の度数を $N_1$ で表す)として、サンプルサイズ $N$ で割ると相対度数 $\frac{N_x}{N}$ が求まります。

　または、np.unique() を使って度数をカウントします。

# 度数を集計
#x_vec, freq_vec = np.unique(ar=x_n, return_counts=True) # (未観測値を補完しない場合)
values, counts = np.unique(ar=x_n, return_counts=True)
print(values[:10])
print(counts[:10])

[0 1 2 3 4 5 6 7 8 9]
[20 23 26 12 17 18 14 16  7 10]

　return_counts=True を指定すると、度数を出力します。ただし、未観測値については配列に含まれません。
　観測されなかった変数の値も配列に含める場合は、次の処理を追加します。

# 未観測値を補完
freq_vec = np.zeros_like(x_vec, dtype='int') # 受け皿を作成
print(freq_vec[:10])
freq_vec[values] = counts.copy() # 観測値を格納
print(freq_vec[:10])

[0 0 0 0 0 0 0 0 0 0]
[20 23 26 12 17 18 14 16  7 10]

　確率変数 $x$ が0以上の離散値をとる場合、全ての要素が 0 の配列を作成して、階級値 values をインデックスとして使い、度数 counts の値を代入します。

　クラスタごとに集計する場合は、次のように処理します。

# クラスタごとの度数を集計
cluster_freq_lt = [
    np.array([np.sum(x_n[s_n == k] == x) for x in x_vec]) for k in range(K)
]
print(len(cluster_freq_lt))
print(cluster_freq_lt[0][:10])

# クラスタごとの相対度数を計算
cluster_relfreq_lt = [
    cluster_freq_lt[k] / N for k in range(K)
]
print(len(cluster_relfreq_lt))
print(cluster_relfreq_lt[0][:10])

4
[20 21 17  3  2  0  0  0  0  0]
4
[0.06666667 0.07       0.05666667 0.01       0.00666667 0.
 0.         0.         0.         0.        ]

　サンプル x_n からクラスタ番号が k のサンプルを取り出して、サンプル値が x の要素数をカウントして、度数をリストに格納します。
　クラスタごとに、度数をサンプルサイズで割って、相対度数をリストに格納します。

乱数の作図

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

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

# 階級幅を設定
bin_size = 1.0

# サンプルの度数を作図
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
plt.bar(
    x=x_vec, height=freq_vec, 
    width=bin_size, align='center', 
    color='#00A968'
) # 度数
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('frequency') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('mixture Poisson distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.show()

　棒グラフの描画関数 plt.bar() でヒストグラムを作成します。横軸の引数(第1引数) x に階級値、縦軸の引数(第2引数) height に度数を指定します。

　または、集計前のサンプルデータを用いて、ヒストグラムを作成します。

# x軸の範囲を設定
x_min = 0

# 階級数を設定
bin_num = int(x_max - x_min + 1)
#bin_num = len(x_vec)

# サンプルの度数を作図
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
plt.hist(
    x=x_n, 
    bins=bin_num, range=(x_min-0.5*bin_size, x_max+0.5*bin_size), 
    color='#00A968'
) # 度数
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('frequency') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('Negative Binomial Distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.show()

　ヒストグラムは、plt.hist() で描画できます。データの引数 x にサンプルデータ、バーの数の引数 bins に階級数、バーの範囲の引き数 range に階級の最小値と最大値を(階級値が離散値となるように調整して)指定します。

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

# サンプルの度数を作図
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
for k in range(K):
    cluster_freq_vec = cluster_freq_lt[k] # 度数
    bottom_vec = np.sum(cluster_freq_lt[:k], axis=0) # 累積度数
    plt.bar(
        x=x_vec, bottom=bottom_vec, height=cluster_freq_vec, 
        label=f'$k = {k+1}$'
    ) # 累積度数
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('frequency') # 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$ の度数を指定します。

　混合ポアソン分布に従う確率を計算します。装飾用の処理です。

## 混合ポアソン分布の確率を計算

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

# 周辺確率を計算
prob_vec = np.sum(weighted_prob_lt, axis=0)

　詳しくは「分布の作図」を参照してください。

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

# サンプルの相対度数を作図
plt.figure(figsize=(8, 6), dpi=100, facecolor='white')
plt.bar(
    x=x_vec, height=relfreq_vec, 
    color='#00A968', alpha=0.5, 
    label='random number'
) # 相対度数
plt.bar(
    x=x_vec, height=prob_vec, 
    facecolor='none', edgecolor='green', linewidth=1.0, linestyle='--', 
    label='generator'
) # 生成確率
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('relative frequency') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('mixture Poisson distribution', fontsize=20) # 図タイトル
plt.legend(title='distribution') # 凡例
plt.grid() # グリッド線
plt.show()

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

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

　第2軸(右側)に度数目盛を表示するグラフを作成します。

# 相対度数軸の範囲を設定
u = 0.05
relfreq_max = np.max(freq_vec/N)
relfreq_max = np.ceil(relfreq_max /u)*u # u単位で切り上げ
#relfreq_max = 0.1

# サンプルの相対度数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white')
fig.suptitle('mixture Poisson distribution', fontsize=20)

# 第1軸の図
ax.bar(
    x=x_vec, height=relfreq_vec, 
    color='#00A968', alpha=0.5, 
    label='random number'
) # 相対度数
ax.bar(
    x=x_vec, height=prob_vec, 
    facecolor='none', edgecolor='green', linewidth=1.0, linestyle='--', 
    label='generator'
) # 確率
ax.set_xlabel('$x$')
ax.set_ylabel('relative frequency, probability') # 1軸ラベル
ax.set_title(param_lbl, loc='left')
ax.legend(title='distribution')
ax.grid()
ax.set_ylim(ymin=0.0, ymax=relfreq_max) # (目盛の共通化用)

# 2軸用のインスタンスを作成
ax2 = ax.twinx()

# 相対度数軸(1軸)目盛を取得
relfreq_vals = ax.get_yticks()

# 度数軸(2軸)目盛に変換
freq_vals = relfreq_vals * N

# 第2軸の図
ax2.set_yticks(ticks=freq_vals, labels=[f'{y:.1f}' for y in freq_vals]) # 2軸目盛
ax2.set_ylabel('frequency') # 2軸ラベル
ax2.set_ylim(ymin=0.0, ymax=relfreq_max*N) # (目盛の共通化用)

plt.show()

　(2軸の図は、pyplot.function() の記法では作成できないようなので、axes.function() の記法で作成します。)

　第2軸は、ax.twinx() を使って描画できます。x軸を共有するの描画用のインスタンスが出力して、第1軸と同様に軸の設定を追加します。
　2つの軸の目盛の位置を共通化するために、第1軸の目盛の値を ax.get_yticks() で取り出して、第2軸の目盛の値として ax.set_yticks() で設定します。この例では、第1軸は相対度数なので、データ数(サンプルサイズ)を掛けて度数に変換して第2軸に使います。
　軸の範囲を ax.set_ylim() で対応させます。

　サンプルクラスタ(カテゴリ分布の乱数)のヒストグラムを作成します。

# k軸の値を作成
k_vec = np.arange(start=1, stop=K+1, step=1)

# カラーマップを作成:(配色の共通化用)
cmap = plt.get_cmap('tab10') # カラーマップを指定
color_num = 10               # カラーマップの色数を設定
color_lt = [cmap(k%color_num) for k in range(K)] # クラスタごとの色

# ラベル用の文字列を作成
N_str  = ', '.join([f'{np.sum(s_n == k)}' for k in range(K)])
pi_str = ', '.join([f'{pi}' for pi in pi_k])
param_lbl = f'$N = ({N_str}), K = {K}, \\pi = ({pi_str})$'

# サンプルクラスタの相対度数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white')
fig.suptitle('Categorical distribution', fontsize=20)
ax.bar(
    x=k_vec, height=np.sum(s_nk, axis=0)/N, 
    color=color_lt, alpha=0.5, 
    zorder=1
) # 相対度数
ax.bar(
    x=k_vec, height=pi_k, 
    facecolor='none', edgecolor=color_lt, linewidth=1.0, linestyle='--', 
    zorder=2
) # 生成確率
ax.set_xticks(ticks=k_vec)
ax.set_xlabel('$k$')
ax.set_ylabel('$\\frac{N_k}{N}, p(s_n = k \\mid \\pi)$')
ax.set_title(param_lbl, loc='left')
ax.grid(zorder=0)
plt.show()

　カテゴリ分布に従う確率を破線の棒グラフで示します。

　塗りつぶしの各バーの高さは、各クラスタ $k$ の相対度数 $\frac{N_k}{N}$ を表します。ここで、 $N_k$ はクラスタ $k$ の度数(サンプル数)です。
　破線の各バーの高さは、各クラスタ $k$ が割り当てられる確率( 各データ $x_n$ のクラスタが $k$ となる確率) $p(s_n = k \mid \boldsymbol{\pi})$ を表します。

　クラスタごとにサンプルデータ(ポアソン分布の乱数)のヒストグラムを作成します。

# カラーマップを作成:(配色の共通化用)
cmap = plt.get_cmap('tab10') # カラーマップを指定
color_num = 10               # カラーマップの色数を設定

# ラベル用の文字列を作成
N_str = ', '.join([f'{np.sum(s_n == k)}' for k in range(K)])
param_lbl = f'$N = ({N_str}), K = {K}$'

# サンプルデータの相対度数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white')
fig.suptitle('Poisson distribution', fontsize=20)

for k in range(K):
    
    # クラスタごとの値を作成
    lmd = lambda_k[k] # パラメータ
    cluster_N        = np.sum(s_n == k) # サンプル数
    cluster_freq_vec = np.array([np.sum(x_n[s_n == k] == x) for x in x_vec]) # 度数
    cluster_prob_vec = poisson.pmf(k=x_vec, mu=lmd) #生成確率
    color_tp = cmap(k%color_num) # 色
    
    ax.bar(
        x=x_vec, height=cluster_freq_vec/cluster_N, 
        color=color_tp, alpha=0.5, 
        label=f'$k = {k+1}, \\lambda_{{k}} = {lmd:.2f}$', 
        zorder=1
    ) # 相対度数
    ax.bar(
        x=x_vec, height=cluster_prob_vec, 
        facecolor='none', edgecolor=color_tp, linewidth=1.0, linestyle='--', 
        zorder=2
    ) # 生成確率
    plt.plot(
        x_vec, cluster_prob_vec, 
        color=color_tp, linewidth=1.0, 
        zorder=3
    ) # 生成確率
    plt.scatter(
        x=x_vec, y=cluster_prob_vec, 
        color=color_tp, s=20, 
        zorder=4
    ) # 生成確率

ax.set_xlabel('$x$')
ax.set_ylabel('$\\frac{N_{kx}}{N_k}, p(x_n = x \\mid s_n = k, \\lambda_k)$')
ax.set_title(param_lbl, loc='left')
ax.legend(title='parameter')
ax.grid(zorder=0)
plt.show()

　クラスタごとに(積み上げずに並べた)ポアソン分布に従う(重み付けしていない)確率を破線の棒グラフや折れ線グラフで示します。

　塗りつぶしの各バーの高さは、各クラスタ $k$ における各値 $x$ の相対度数 $\frac{N_{kx}}{N_k}$ を表します。ここで、 $N_{kx}$ はクラスタ $k$ における値 $x$ の度数です。
　破線の各バーの高さは、割り当てられたクラスタ $s_n = k$ に対応する値が $x$ となる確率 $p(x_n = x \mid s_n = k, \lambda_k)$ を表します。