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

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

指数分布について

連続型確率分布の一つに、機械が故障してから次に故障するまでの期間や、災害が起こってから次に起こるまでの期間のように、次に何かが起こるまでの期間が従う分布の指数分布があります。この記事では、pythonを使った指数分布の求め方について紹介しています。

指数分布とは

ある期間に平均して$\lambda$(ラムダ)回起こる現象が、次に起こるまでの期間$X$が指数分布に従うとき、$X=x$となる確率密度関数$f(x)$は次の式で表されます。


f(x) = \left\{
\begin{array}{ll}
\lambda e^{-\lambda x} & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.

例えば、1時間に平均20人のお客さんが送るお店に、あるお客がきて次のお客が来るまでの時間が10分となる確率密度は以下のように計算します。この場合は、10分は$\frac{1}{6}$時間となります。


f(\frac{1}{6})=\lambda e^{-\lambda x} = 20\times e^{-20\times \frac{1}{6}} = 0.714

指数分布に従う期待値と分散

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


E(X) = \frac{1}{\lambda}

V(X) = \frac{1}{\lambda^2}

例えば、先ほどの例の場合、$\lambda$は20となるので、期待値と分散は以下になります。


E(X) = \frac{1}{\lambda} =\frac{1}{20}

V(X) = \frac{1}{\lambda^2} =\frac{1}{400}

指数分布に従う累積分布関数

積分布関数は、ある期間に平均して$\lambda$回起こる現象が次に起こるまでの期間を$X$としたとき、期間$X$が$x$以下となる確率のことを言い、次のようになります。


F(x) = P (X \leq x) = 1 - e^{-\lambda x}

例えば、1時間に平均20人が来るお店に、ある客が来てから次の客が来るまでの時間が10分以内である確率は以下のように計算します。


F(\frac{1}{6}) = P (X \leq \frac{1}{6}) = 1 - e^{-20\times \frac{1}{6}} = 1 - 0.03569 = 0.96431

指数分布に従う確率密度の求め方

Pythonを使って指数分布に従う確率密度を求める場合は、以下のようになります。

import math

e = math.e
lambda1 = 20
x = 1/6

fx = lambda1*e**(-lambda1*x)
print("f(1/6) = ", fx)

f(1/6) = 0.7134798669450483

指数分布に従う累積分布関数の求め方その1

Pythonを使って、指数分布に従う累積分布関数の確率を求める場合は、以下のように解析的に求める方法があります。

e = math.e
lambda1 = 20
x = 1/6

Fx = 1 - e**(-lambda1*x)
print("P(X <= 1/6) = ", Fx)

P(X <= 1/6) = 0.9643260066527476

指数分布に従う累積分布関数の求め方その2(np.random.exponential)

Numpyのnp.random.exponentialを使用すると、指数分布に従うデータを生成することができます。なので、以下のように10分以内(1/6)の確率は以下のように計算することができます。

import numpy as np
size = 100000
lambda1 = 20
x = 1/6
beta = 1/lambda1

Fx = sum(np.random.exponential(beta, size) <= x)/size
print("P(X <= 1/6) = ", Fx)

P(X <= 1/6) = 0.96435

指数分布に従う累積分布関数の求め方その3(expon.rvs)

scipyのexpon.rvsを使用すると、指数分布に従うデータを生成することができます。なので、以下のように10分以内(1/6)の確率は以下のように計算することができます。

from scipy.stats import expon
size = 100000
lambda1 = 20
x = 1/6
beta = 1/lambda1

Fx = sum(expon.rvs(scale=beta, size=size) <= x)/size
print("P(X <= 1/6) = ", Fx)

P(X <= 1/6) = 0.96468

参考

統計学の時間 | 統計WEB