python で幾何分布の計算を行う方法を調べてみました。結果を記載します。
幾何分布とは何か
統計学入門から
X回目で初めてOOが起こる確率が取る分布を幾何分布といいます。
統計学入門 (基礎統計学Ⅰ) | 東京大学教養学部統計学教室 |本 | 通販 | Amazon には、以下のように記載されています。
-
幾何分布とは
幾何分布は、時間を1,2,3,…離散的に考えるとき、初めてのSを得るまで待つ時間の長さの確率分布であり、離散的な待ち時間分布 waiting time distribution と呼ばれる。
-
現実での例
災害の到来、 失敗する回数
-
負の二項分布との関係
負の二項分布 の 公式は以下で、k = 1 の場合、幾何分布となる。
$$ f(x)= {k+x-1}C = 0, 1, 2 … $$ }p^kq^{x}, {x
wikipediaから
幾何分布 - Wikipedia に幾何分布の性質として、無記憶性の記載があります。
幾何分布の重要な性質として、無記憶性と呼ばれるものがある。 中略..
各種のギャンブルにおいて負けが続くと、しばしば「運がたまっている」とか「そろそろ勝ちが巡ってくる」といった考えに陥りがちである。しかし、試行の独立性を仮定する限りにおいては、この考えは誤謬であり、負けが続いているという情報は未来の確率に何の影響も与えないということが、無記憶性からわかる。 この逆、すなわち無記憶性を持つ離散型確率分布が幾何分布のみであることも、比較的容易に示される。
幾何分布の計算ができるライブラリについて
他の分布と同様に、numpy と scipy で計算ができます。
-
numpy.random.geometric
-
scipy.stats.geom
numpy.random.geometric を使う
numpy.random.geometric は 負の二項分布に従う、配列を取得できます。
* ヒストグラム描画
# ヒストグラム描画
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# グラフサイズの指定
plt.figure(figsize=(6,6),dpi=100)
#成功確率 0.05 試行回数 10000
x1 = np.random.geometric(0.05,10000)
count, bins, ignored = plt.hist(x1, 30, normed=True,histtype="bar",color="blue", alpha=0.3)
#成功確率 0.1 試行回数 10000
x2 = np.random.geometric(0.1,10000)
count, bins, ignored = plt.hist(x2, 30, normed=True,histtype="bar",color="yellow", alpha=0.3)
#成功確率 0.2 試行回数 10000
x3 = np.random.geometric(0.2,10000)
count, bins, ignored = plt.hist(x3, 30, normed=True,histtype="bar",color="green", alpha=0.3)
# 描画
plt.show()
- n回で成功する確率を求める
5回で成功する確率、5回以上で成功する確率 を求めます。
# ヒストグラム描画
import numpy as np
#成功確率 0.05 試行回数 10000
x1 = np.random.geometric(0.05,10000)
# 5回目で成功する確率
print((x1 == 5).sum() / 10000)
# 5回以上で成功する確率
print((x1 >= 5).sum() / 10000)
#成功確率 0.1 試行回数 10000
x2 = np.random.geometric(0.1,10000)
# 5回目で成功する確率
print((x2 == 5).sum() / 10000)
# 5回以上で成功する確率
print((x2 >= 5).sum() / 10000)
#成功確率 0.2 試行回数 10000
x3 = np.random.geometric(0.2,10000)
# 5回目で成功する確率
print((x3 == 5).sum() / 10000)
# 5回以上で成功する確率
print((x3 >= 5).sum() / 10000)
0.0398
0.8063
0.0625
0.6534
0.0855
0.4179
scipy.stats.geom を使う
scipy.stats.geom を使って、負の2項分布のグラフを描画してみます。
%matplotlib inline
from scipy.stats import geom
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1)
p = 0.1
# mements は 文字列 mvsk を指定すると、 m > 平均 v > 分散 > s > 歪度（わいど) k > 尖度(せんど) が返る
mean, var, skew, kurt = geom.stats(p, moments='mvsk')
print("平均", mean)
print("分散", var)
print("歪度", skew)
print("尖度", kurt)
# グラフ描画
x = np.arange(geom.ppf(0.01, p), geom.ppf(0.99, p))
ax.plot(x, geom.pmf(x, p), 'bo', ms=8, label='geom pmf')
ax.vlines(x, 0, geom.pmf(x, p), colors='b', lw=5, alpha=0.5)
rv = geom(p)
ax.vlines(x, 0, rv.pmf(x), colors='k', linestyles='-', lw=1, label='frozen pmf')
ax.legend(loc='best', frameon=False)
plt.show()
平均 10.0
分散 90.0
歪度 2.0027758514399734
尖度 6.011111111111111
- ランダムな配列を返す
scipy.stats.geom.rvsでnumpy.random.geometricと同様の幾何分布に従うランダムな配列を取得できます。
p = 0.1
from scipy.stats import geom
r = geom.rvs(p, size=10000)
# 5回目で成功する確率
print((r == 5).sum() / 10000)
0.067
参考
以下、参考にした記事になります。
以上です。
コメント