はじめに

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

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

【前の内容】

zenn.dev

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
ポアソン分布の乱数生成
- サンプリング
- サンプルサイズの影響
  - サンプルサイズと形状の関係
参考文献
おわりに

ポアソン分布の乱数生成

　ポアソン分布(Poisson distribution)に従う乱数を生成して、グラフやアニメーションで確認します。この記事では、PythonのMatplotlibライブラリを利用して作図します。
　ポアソン分布については「ポアソン分布の定義式 - からっぽのしょこ」、グラフ作成については「【Python】ポアソン分布の作図 - からっぽのしょこ」、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 とします。

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

乱数の生成

　サンプルサイズ $N$ を指定して、ポアソン分布の乱数 $x_n$ を生成します。

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

# ポアソン分布の乱数を生成
x_n = np.random.poisson(lam=lmd, size=N)
print(x_n[:10])

[5 5 2 4 9 7 4 5 4 3]

　サンプルサイズ(データ数) $N$ を指定します。
　 $N$ 個のサンプル(乱数)の集合 $\mathbf{x} = \{x_1, \cdots, x_N\}$ を作成します。各要素は、単位時間における発生回数(非負の整数) $x_n \in \{0, 1, 2, \cdots\}$ です。
　ポアソン分布の乱数は、np.random.poisson() で生成できます。パラメータの引数 lamに $\lambda$ 、サンプルサイズの引数 size に $N$ を指定します。

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

乱数の集計

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

# x軸の範囲を指定
#x_max = 30
u = 5.0
x_max = np.max(x_n)
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[:10])

[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]

　集計や作図に用いるため、範囲を指定して、単位時間における発生回数(非負の整数) $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])

[ 1  8 14 17 20 17  8  9  4  1]
[0.01 0.08 0.14 0.17 0.2  0.17 0.08 0.09 0.04 0.01]

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

乱数の作図

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

# サンプルの度数を作図
plt.figure(figsize=(8, 6), facecolor='white')
plt.bar(
    x=x_vec, height=freq_vec, 
    color='#00A968', width=1.0
) # 度数
plt.xticks(ticks=x_vec) # x軸目盛
plt.grid() # グリッド線
plt.xlabel('x') # x軸ラベル
plt.ylabel('frequency') # y軸ラベル
plt.title(f'$N = {N}, \\lambda = {lmd}$', loc='left') # タイトル
plt.suptitle('Poisson distribution', fontsize=20) # 図タイトル
plt.show()

　棒グラフの描画関数 plt.bar() の縦軸の引数(第2引数) height に度数を指定してヒストグラムを作成します。

　ポアソン分布に従う確率を計算します。

# ポアソン分布を計算
prob_vec = poisson.pmf(k=x_vec, mu=lmd)

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

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

# サンプルの相対度数を作図
plt.figure(figsize=(8, 6), 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.xticks(ticks=x_vec) # x軸目盛
plt.grid() # グリッド線
plt.xlabel('x') # x軸ラベル
plt.ylabel('relative frequency') # y軸ラベル
plt.title(f'$N = {N}, \\lambda = {lmd}$', loc='left') # タイトル
plt.suptitle('Poisson distribution', fontsize=20) # 図タイトル
plt.legend(title='distribution') # 凡例
plt.show()

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

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

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

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

# サンプルの相対度数を作図
fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
fig.suptitle('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_xticks(ticks=x_vec) # x軸目盛
ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('relative frequency, probability') # 1軸ラベル
ax.set_title(f'$N = {N}, \\lambda = {lmd}$', loc='left')
ax.legend(title='distribution')
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() で対応させます。

　以上で、乱数生成に関する基本的な処理を確認しました。

スポンサードリンク

サンプルサイズの影響

　次は、サンプルサイズを増やしてグラフの変化をアニメーションで確認します。
　作図コードについては「Probability-Distribution/code/poisson/random_number.py at main · anemptyarchive/Probability-Distribution · GitHub」を参照してください。