Python で指数分布の計算をする方法を調べてみました。
結果を記載します。
指数分布とは
そもそも、指数分布とは何かについて調べると以下の記事がヒットします。
個人的には、指数分布 Web 解析 なので左記キーワードを含めると以下の記事が見つかります。
書籍
統計学入門 - 東京大学出版会 には以下の内容が記載されています。引用します。
指数分布は
exponetial distributionは、(5.10) で述べたように確率密度関数
$$ f(x) = \lambda e^{-\lambda x} (x >= 0), 0 (x < 0) $$
で定義される連続型の分布である。幾何分布において $ q = e^{-\lambda} $ すなわち $ \lambda = -log q $ とおくと、形式上この確率分布の主要な部分 $ e^{-\lambda x} $ が得られるから、この確率分布は連続的な待ち時間分布の性質をもつことがわかる。
自然科学の統計学 - 東京大学出版会 には以下の内容が記載されています。引用します。
指数分布は、ポアソン分布に従って起きる事象の生起感覚を表現する分布である。指数分布は、さまざまなサービスの待ち時間を表現する分布として用いられる。
実世界での例
実世界で、指数分布となる事象には以下があります。
- サービスの待ち行列の分布
- システムの寿命
- 災害が発生するまでの年数
- 雨滴の大きさ
- Wikipedia へのページリクエスト間の時間
Web
指数分布で検索をかけると以下のサイトが見つかりました。
* 指数分布 - Wikipedia
指数分布の計算方法
Python では、以下の関数でできることがわかりました。
numpy.random.exponential を使う
パラメータmu(分布の平均値)、サンプル数 n の指数分布 を求めることができます。
numpy.random.exponential は指数分布に従う配列を返し、その値をもとにグラフ描画ができます。
関数に指定できる引数はなかなか柔軟で以下のような指定ができます。
引数 scale のみ
第一引数scale のみ指定ができます。 この場合、スカラー値が取得できます。
スカラー値とは以下のことだと思います。正直あまりよく理解できておりません。
スカラー (数学) - Wikipedia
import numpy as np
import pylab as plt
target = 500
beta = 1.0/target
Y = np.random.exponential(beta)
print(beta, Y)
0.002 0.00026848474693625724
引数 scale にはリストも指定可能
引数 scale には、リストもタプルも指定が可能です。
戻り値が毎回異なりますが、理由がよくわかっていません。
import numpy as np
import pylab as plt
target = 500
beta = [1.0/target, 2.0/target, 3.0/target]
Y = np.random.exponential(beta)
print(Y)
beta = (1.0/target, 2.0/target, 3.0/target)
Y = np.random.exponential(beta)
print(Y)
[ 0.00035789 0.00079873 0.0059433 ]
[ 0.00070439 0.00116696 0.01094726]
size を指定する
第二引数size を指定すると、指数分布に従う配列が取得できます。
scale はポアソン分布のラムダ(ある期間に平均してn回起こるのn) の逆数を示します。
これが一般的な指数分布に思われますので、グラフとして描画します。
import numpy as np
import pylab as plt
# target が ラムダ ある期間に平均して10回起こる
target = 10
beta = 1.0/target
# 要素数10000 ある期間に平均して10回 指数分布
Y = np.random.exponential(beta, 10000)
plt.hist(Y, normed=True, bins=200,lw=0,alpha=.8)
# Y の最大値をプロットしておく
plt.plot([0,max(Y)],[target,target],'r--')
plt.ylim(0,target*1.1)
plt.show()
### size に3要素のタプルを指定する
size には 3要素の タブル を指定することができます。
戻り値が3次元配列で返ります。
import numpy as np
import pylab as plt
# target が ラムダ ある期間に平均して10回起こる
target = 10
beta = 1.0/target
# (m=3, n=5, k=2) を指定する
Y = np.random.exponential(beta, (3, 5, 2))
print(Y)
[[[ 0.03303321 0.31016278]
[ 0.2045412 0.24228506]
[ 0.12344505 0.15427801]
[ 0.21190096 0.02285032]
[ 0.00356998 0.08352543]]
[[ 0.05553856 0.06446413]
[ 0.11237701 0.12704291]
[ 0.00912467 0.1277612 ]
[ 0.01687101 0.19313137]
[ 0.0858672 0.08564713]]
[[ 0.0119443 0.0315812 ]
[ 0.00289961 0.07208328]
[ 0.04683225 0.04912672]
[ 0.08552624 0.04571101]
[ 0.01842625 0.0141593 ]]]
scipy.stats.expon を使う
続いてscipy.stats.expon を使ってみます。
モーメントを求める
expon.stats(moments='mvsk') でモーメントを求めることができます。
これはscipy の他の分布関数と同様です。
歪度、尖度 等は、指数分布 - Wikipedia に取りうる値の記載があり参考になりました。
from scipy.stats import expon
mean, var, skew, kurt = expon.stats(moments='mvsk')
#平均1、分散1、歪度 2、尖度 6
print(mean, var, skew, kurt)
1.0 1.0 2.0 6.0
指数分布のグラフを描く
scipy.stats.expon — SciPy v1.1.0 Reference Guide を参考に、グラフを描きます。
from scipy.stats import expon
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1)
# 指定した区間をN等分した数列を生成する
# N 100 を増減させると、曲線の滑らかさが変わる
# ppf は percent point function の略
x = np.linspace(expon.ppf(0.01),expon.ppf(0.99), 100)
# expon は、Exponential distribution の expon だと思う
ax.plot(x, expon.pdf(x),'r-', lw=5, alpha=0.6, label='expon pdf')
# frozen > 固定した
rv = expon()
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
# ランダムな要素数1000 の指数分布に従う配列を取得
# loc (指数分布のx軸の位置)、scale (1 / lambda) の rvs に指定可能。 というか他の関数にも指定できる
r = expon.rvs(size=1000)
# ヒストグラム描画
ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)
ax.legend(loc='best', frameon=False)
plt.show()
参考
以下参考になりました。
- Manipulating the numpy.random.exponential distribution in Python - Stack Overflow
- scipy.stats.expon — SciPy v1.1.0 Reference Guide
正直、私には難しいかったですが、ポアソン分布との関係があることと、Web の待ち時間と関係があることがわかったのが良い経験になりました。
以上です。
コメント