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

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

Pythonを使った幾何分布の求め方

離散型確率分布には、成功確率がpの時のベルヌーイ試行を繰り返すとき、初めて成功するまでの試行回数に従う確率分布として幾何分布があります。この記事では、pythonを使った幾何分布の求め方について解析的に求める場合と、数値的に求める場合について紹介しています。

幾何分布とは

成功確率がpである独立なベルヌーイ試行を繰り返す時、初めて成功するまでの試行回数Xが従う確率分布を幾何分布といいます。


P(X = k) = (1-p)^{k-1} p

例えば、サイコロを2投目で初めて1が出る確率は以下のようになります。


P(X = 2) = (1-\frac{1}{6})^{2-1} \times\frac{1}{6} = \frac{5}{36} = 0.1388

幾何分布の期待値と分散

確率変数Xが成功確率pの幾何分布に従う時、Xの期待値と分散は以下のようになります。


E(X) = \frac{1}{p}\\
V(X) = \frac{1-p}{p^2}

例えば、けん玉を10回やった結果、4回は失敗する場合、次に失敗するまでの回数の期待値と分散は以下のようになります。


E(X) = \frac{1}{\frac{4}{10}} =\frac{5}{2} = 2.5\\
V(X) = \frac{1-\frac{4}{10}}{(\frac{4}{10})^2} = \frac{\frac{3}{5}}{\frac{4}{25}} = \frac{3}{5}\times\frac{25}{4} = \frac{75}{20} = 3.75

幾何分布を解析的に求める場合

幾何分布の確率を解析的に求める場合は、以下のように求めることができます。

p = 1/6
k = 2

px = p*(1 - p)**(k - 1)
print("P(X = 2) = ", px)

P(X = 2) = 0.1388888888888889

np.random.geometricを使って幾何分布を試行し数値的に確率を求める場合

np.random.geometricを使用すると、幾何分布を実際に試行することができます。以下の例では、確率が1/6の時の幾何分布を10000セット実行し、サイコロ2投目で1が出る回数を出力しています。そして、実際に実行したデータを使って幾何分布の確率を出しています。10000セットくらい実行したデータから求めた確率だと解析的に求めた確率に近くなっています。

import numpy as np
p = 1/6
k = 2

px = np.random.geometric(p,10000)
print("P(X = 2) = ", (px == k).sum() / 10000)

P(X = 2) = 0.1339

scipyのgeomを使って幾何分布を試行し数値的に確率を求める場合

scipyのgeomを使用することでも、幾何分布を実際に試行することができます。以下の例では、確率が1/6の時の幾何分布を10000セット実行し、サイコロ2投目で1が出る回数を出力しています。

from scipy.stats import geom
p = 1/6
k = 2

px = geom.rvs(p, size=10000)
print("P(X = 2) = ", (px == k).sum() / 10000)

P(X = 2) = 0.1383