データ分析系プログラマーのブログ

主にPythonを使ったデータ分析や機械学習をやっています。

Pythonを使ったポアソン分布の求め方

離散型確率分布のうち、結果が2つしかないような試行のことをベルヌーイ試行と言いますが、このうちnが大きく、pが小さい時をポアソン分布と言います。この記事では、pythonのライブラリを使ってポアソン分布の確率を解析的に求める場合と、数値的に求める場合について紹介しています。

ポアソン分布とは

ベルヌーイ試行のうちnが大きく、pが小さい時、「np = 一定(λ)」と考えることができ、ある期間にX回起きる確率分布をポアソン分布呼びます。


\lambda = np\\
e = 2.7182818 \dots\\
P(X = k) = \frac{e^{-\lambda}\lambda^k}{k!}

例えば、100個に1個の当たりのあるお菓子を10個買った時に、当たりのお菓子が2個含まれる場合の確率は以下のようになります。


\lambda = 10\times\frac{1}{100} = 0.1\\
e = 2.7182818 \dots\\
P(X = 2) = \frac{e^{-0.1}\times0.1^2}{2!} = 0.0045

ポアソン分布の期待値と分散

確率変数Xがポアソン分布に従う時、Xの期待値と分散は以下のようにλなります。


E(X) = \lambda\\
V(X) = \lambda

例えば、100個に1個の当たりのあるお菓子を10個買った時の期待値と分散はともに以下のようになります。


E(X) = \lambda = np = 10\times\frac{1}{100} = 0.1\\
V(X) = \lambda = 0.1

ポアソン分布を解析的に求める場合その1

ポアソン分布の確率を解析的に求める場合は、mathにネイピア数が用意されているので以下のように求めることができます。

import math

e = math.e
n = 10
p = 1/100
lambda1 = n*p
k = 2

px = (e**(-lambda1)*lambda1**k)/math.factorial(k)
print("P(X = 2) = ", px)

P(X = 2) = 0.004524187090179798

ポアソン分布を解析的に求める場合その2(poisson)

pythonのライブラリのscipyのpoissonを使用すると引数を指定するだけでポアソン分布の確率を求めることができます。

from scipy.stats import poisson

n = 10
p = 1/100
lambda1 = n*p
k = 2

px = poisson.pmf(k,lambda1)
print("P(X = 2) = ", px)

P(X = 2) = 0.004524187090179801

np.random.poissonを使ってポアソン分布を試行し数値的に確率を求める場合

np.random.poissonを使用すると、ポアソン分布を実際に試行することができます。以下の例では、ラムダが0.1の時のポアソン分布を500セット実行し、あたりの出る回数を出力しています。そして、実際に実行したデータを使ってポアソン分布の確率を出しています。100000セットくらい実行したデータから求めた確率だと解析的に求めた確率に近くなっています。

import numpy as np

n = 10
p = 1/100
lambda1 = n*p

s = np.random.poisson(lambda1, 500)

print(s)
[0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 0 0
 0 0 0 0 0 2 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 1
 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 1 0 1 0 0 0 0 2 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0
 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0
 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
s = sum(np.random.poisson(lambda1, 100000) == 2)/100000
print("P(X = 2) = ", s)

P(X = 2) = 0.00466

参考

統計学の時間 | 統計WEB