からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

トップ > 確率分布 > 【Python】ポアソン分布の計算

【Python】ポアソン分布の計算

確率分布 確率分布-ポアソン分布(全般) 確率分布-ポアソン分布(全般)-ポアソン分布 Python

はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

ポアソン分布の計算

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

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

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

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

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


確率の計算

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

パラメータの設定

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

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

4.0

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

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

# 確率変数の値を指定
x = 2.0
print(x)

2.0

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

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

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

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

# 定義式により確率を計算
prob = np.exp(-lmd) * lmd**x / gamma(x + 1.0)
print(prob)

0.14652511110987343

　ポアソン分布は、次の式で定義されます。

\displaystyle \mathrm{Poisson}(x \mid \lambda) = \exp(-\lambda) \frac{\lambda^x}{x!}

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

　対数をとった定義式を実装して計算します。

# 対数をとった定義式により確率を計算
log_prob = -lmd + x * np.log(lmd) - loggamma(x + 1.0)
prob     = np.exp(log_prob)
print(log_prob)
print(prob)

-1.9205584583201643
0.14652511110987343

　対数をとったポアソン分布は、次の式になります。

\displaystyle \log \mathrm{Poisson}(x \mid \lambda) = - \lambda + x \log \lambda - \log x!

　自然対数は np.log()、対数をとったガンマ関数は sp.special.lgamma() で計算できます。
　指数と対数の関係 \exp(\log x) = x より、計算結果の指数をとると確率が得られます。

\displaystyle \mathrm{Poisson}(x \mid \lambda) = \exp \Bigl( \log \mathrm{Poisson}(x \mid \lambda) \Bigr)

　ガンマ関数は、階乗の計算なので値が大きくなり、発散することがあります。

# 発散する例
print(gamma(np.arange(170, 174)))

[4.26906801e+304 7.25741562e+306             inf             inf]

　対数ガンマ関数を使うことで、発散を回避して途中の計算を行えることもあります。

# 対策例
print(loggamma(np.arange(170, 174)))

[701.43726381 706.57306225 711.7147258  716.86222028]


クラスによる計算

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

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

0.14652511110987343

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

　logpmf() だと対数をとった確率を計算します。

# クラスにより確率を計算
log_prob = poisson.logpmf(k=x, mu=lmd)
prob = np.exp(log_prob)
print(log_prob)
print(prob)

-1.9205584583201643
0.14652511110987343

　計算結果の指数をとると確率が得られます。

　ここまでで、ポアソン分布の確率の計算を確認しました。続いて、ポアソン分布に関する各種の計算を確認していきます。

スポンサードリンク

統計量の計算

　次は、ポアソン分布の統計量を計算します。
　計算式については「ポアソン分布の平均と分散の導出：定義式を利用 - からっぽのしょこ」を参照してください。

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

　計算式を実装して期待値を計算します。

# 計算式により期待値を計算
E_x = lmd
print(E_x)

4.0

　ポアソン分布の期待値はパラメータ \lambda です。

\displaystyle \mathbb{E}[x] = \lambda

　分散を計算します。

# 計算式により分散を計算
V_x = lmd
print(V_x)

4.0

　ポアソン分布の分散もパラメータ \lambda です。

\displaystyle \mathbb{V}[x] = \lambda

　最頻値を計算します。

# 計算式により最頻値を計算
mode_x = np.floor(lmd)
print(mode_x)

4.0

　ポアソン分布の最頻値(モード)は、パラメータ \lambda 以下の最大の整数です。

\displaystyle \mathrm{mode}[x] = \lfloor \lambda \rfloor

　ある値以下の最大の整数は、np.floor() で小数点以下を切り捨てることで求められます。

クラスによる計算

　scipy のクラスメソッドを使って期待値を計算します。

# クラスにより期待値を計算
E_x = poisson.mean(mu=lmd)
print(E_x)

4.0

　ポアソン分布の期待値は、sp.stats.poisson.mean() で計算できます。確率メソッドと同様に引数を指定します。

　分散を計算します。

# クラスにより分散を計算
V_x = poisson.var(mu=lmd)
print(V_x)

4.0

　ポアソン分布の分散は、sp.stats.poisson.var() で計算できます。

スポンサードリンク

モーメントの計算

　最後に、ポアソン分布の歪度と尖度を計算します。
　詳しくはいつか書きたい。

　歪度を計算します。

# 歪度を計算
skew = 1.0 / np.sqrt(lmd)
print(skew)

0.5

　ポアソン分布の歪度は、次の式で計算できます。

\displaystyle \mathrm{Skewness} = \frac{\mathbb{E}[(x - \mu)^3]}{\sigma^3} = \frac{1}{\sqrt{\lambda}}

　ここで、 x の期待値 \mu = \mathbb{E}[x]、標準偏差 \sigma = \sqrt{\mathbb{E}[(x - \mu)^2]} です。
　ルート \sqrt{x} は、np.sqrt() で計算できます。

　尖度を計算します。

# 尖度を計算
kurt = 1.0 / lmd
print(kurt)

0.25

　ポアソン分布の尖度は、次の式で計算できます。

\displaystyle \mathrm{Kurtosis} = \frac{\mathbb{E}[(x - \mu)^4]}{\sigma^4} - 3 = \frac{1}{\lambda}


　この記事では、ポアソン分布の確率などの計算を確認しました。次の記事では、確率分布のグラフを作成します。

参考文献

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

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

Amazon

おわりに

　加筆修正しました。その際に「【Python】ポアソン分布の作図」から記事を分割しました。

　第2稿(前回の全面改修)ではPython版は手を付けられなかったようなので、第3稿(今回の全面改修)ではしっかりとやっていきたい、という方針で各種の分布の既存記事のうち加筆修正した1つ目の分布の1つ目の記事です。

【次の内容】

www.anarchive-beta.com