はじめに

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

　この記事では、ポアソン分布のグラフをPythonを使って作成します。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

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

ポアソン分布の作図

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

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

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

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

# クラスを読込
from scipy.stats import poisson

グラフの作成

　まずは、ポアソン分布のグラフを作成します。

パラメータの設定

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

# パラメータを指定
lmd = 4.0
print(lmd)

4.0

　単位時間における発生回数の期待値(正の値) $\lambda \gt 0$ を指定します。
　lambda は予約語のため変数名として使えないので lmd とします。

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

# x軸の範囲を設定
#x_max = 10
u = 5.0
x_max = lmd
x_max *= 1.5 # 倍率を指定
x_max = np.ceil(x_max /u)*u # u単位で切り上げ

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

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

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

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

確率分布の計算

　 $x$ の値ごとにポアソン分布に従う確率 $\mathrm{Poi}(x \mid \lambda)$ を計算します。

# 確率を計算
prob_vec = poisson.pmf(k=x_vec, mu=lmd)
print(prob_vec[:5])

[0.01831564 0.07326256 0.14652511 0.19536681 0.19536681]

　ポアソン分布の確率は、sp.stats.poisson.pmf() で計算できます。確率変数の引数 k に $x$ (のベクトル)、パラメータの引数 mu に $\lambda$ を指定します。

確率分布の作図

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

# ポアソン分布を作図:棒グラフ
plt.figure(figsize=(8, 6), facecolor='white')
plt.bar(
    x=x_vec, height=prob_vec, 
    color='#00A968'
) # バー
plt.xticks(ticks=x_vec) # x軸目盛
plt.grid() # グリッド線
plt.xlabel('x') # x軸ラベル
plt.ylabel('probability') # y軸ラベル
plt.title(f'$\\lambda = {lmd}$', loc='left') # タイトル
plt.suptitle('Poisson distribution', fontsize=20) # 図タイトル
plt.show()

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

　幹図を作成します。

# ポアソン分布を作図:幹図
plt.figure(figsize=(8, 6), 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.xticks(ticks=x_vec) # x軸目盛
plt.grid() # グリッド線
plt.xlabel('x') # x軸ラベル
plt.ylabel('probability') # y軸ラベル
plt.title(f'$\\lambda = {lmd}$', loc='left') # タイトル
plt.suptitle('Poisson distribution', fontsize=20) # 図タイトル
plt.show()

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

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

# ポアソン分布を作図:折れ線グラフ
plt.figure(figsize=(8, 6), 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.xticks(ticks=x_vec) # x軸目盛
plt.grid() # グリッド線
plt.xlabel('x') # x軸ラベル
plt.ylabel('probability') # y軸ラベル
plt.title(f'$\\lambda = {lmd}$', loc='left') # タイトル
plt.suptitle('Poisson distribution', fontsize=20) # 図タイトル
plt.show()

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

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

スポンサードリンク

複数グラフの作成

　次は、複数のパラメータのポアソン分布を並べて描画します。

パラメータの設定

　複数個のパラメータ $\lambda$ を指定します。

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

# パラメータ数を設定
param_num = len(lambda_vals)
print(param_num)

[ 1.   5.5 10.  16.8]
4

　パラメータ $\lambda$ を数値ベクトルで指定します。
　計算や作図の処理用に、指定したパラメータ数を作成しておきます。

　確率変数の実現値 $x$ を設定します。

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

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

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

　「グラフの作成」のときと同様に処理します。
　この例では、指定したパラメータ lambda_vals の最大値を使って範囲を設定します。

確率分布の計算

　 $\lambda$ と $x$ の組み合わせごとに確率 $\mathrm{Poi}(x \mid \lambda)$ を計算します。

# 確率計算
prob_lt = [
    poisson.pmf(k=x_vec, mu=lambda_vals[i]) for i in range(param_num)
]
print(len(prob_lt))
print(prob_lt[0][:5])

4
[0.36787944 0.36787944 0.18393972 0.06131324 0.01532831]

　パラメータごとに確率を計算してリストに格納します。リスト内包表記により、パラメータ $\lambda$ と変数 $x$ の全ての組み合わせに対して計算できます。

確率分布の作図

　パラメータごとのポアソン分布を重ねて描画します。

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

# パラメータごとにポアソン分布を作図
plt.figure(figsize=(9, 6), facecolor='white')
for i in range(param_num):
    # 値を取得
    lmd = lambda_vals[i]  # パラメータ
    prob_vec = prob_lt[i] # 確率
    color_tp = cmap(i%color_num) # 色
    
    plt.bar(
        x=x_vec, height=prob_vec, 
        color=color_tp, alpha=0.5, 
        label=f'$\\lambda = {lmd}$'
    ) # バー
    plt.plot(
        x_vec, prob_vec, 
        color=color_tp, linewidth=2.0
    ) # 折れ線
    plt.scatter(
        x=x_vec, y=prob_vec, 
        color=color_tp, s=50
    ) # 点
#plt.xticks(ticks=x_vec) # x軸目盛
plt.grid() # グリッド線
plt.xlabel('x')
plt.ylabel('probability')
plt.suptitle('Poisson distribution', fontsize=20)
plt.legend(title='parameter') # 凡例
plt.show()

　「グラフの作成」のときと同様に各種グラフを描画します。
　グラフ(描画関数)ごとに配色(色付けルール)が異なる場合があるので、配色を統一する処理を行います。plt.get_cmap() にカラーマップ名を指定して、カラーマップインスタンス(cmapオブジェクト)を作成します。cmap(整数) の形式で色番号を指定すると、(0から1に正規化された)RBGA値を出力します。
　カラーマップごとに色数が異なり(決まっており)、色数を超える値を指定すると意図しない配色になるので、色付けがサイクルする(例えば10色のカラーマップの場合は1番目と11番目が同じ色になる)ように処理しています。

　パラメータごとに描画領域を分割して分布を描画します。

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

# サブプロットの列数を指定
col_num = 2
row_num = np.ceil(param_num / col_num).astype('int')

# 確率軸の範囲を設定
u = 0.05
prob_max = np.max(prob_lt)
prob_max = np.ceil(prob_max /u)*u # u単位で切り上げ

# パラメータごとにポアソン分布を作図
fig, axes = plt.subplots(nrows=row_num, ncols=col_num, constrained_layout=True, 
                         figsize=(12, 8), facecolor='white')
fig.suptitle('Poisson distribution', fontsize=20)
for i in range(param_num):
    # サブプロットを抽出
    r = i // col_num
    c = i % col_num
    ax = axes[r, c]
    
    # 値を取得
    lmd = lambda_vals[i]  # パラメータ
    prob_vec = prob_lt[i] # 確率
    
    ax.bar(
        x=x_vec, height=prob_vec, 
        color=cmap(i%color_num)
    ) # バー
    #ax.set_xticks(ticks=x_vec) # x軸目盛
    ax.grid() # グリッド線
    ax.set_xlabel('x')
    ax.set_ylabel('probability')
    ax.set_title(f'$\\lambda = {lmd}$', loc='left')
    ax.set_ylim(ymin=0.0, ymax=prob_max) # 確率軸の描画範囲
    
# 不要なサブプロットを非表示
for c in range(c+1, col_num):
    axes[r, c].axis('off')
plt.show()