幾何分布、指数分布、ポアッソン分布のシミュレーション
時間軸上において生起確率が一様である事象が、観測値としては幾つかがまとまって生起したり、あるいはある期間生起が観測されなかったりと、見かけ上の群化が観測される。これを幾何分布のシミュレーションで確認する。幾何分布の時間の連続量版としての指数分布のシミュレーションも行う。また、指数分布に基づいて単位時間内の事象の生起頻度のシミュレーションを行い、ポアッソン分布であることを確認する。
幾何分布
整数値で表される各時点における事象の生起確率が
で一定であり、各時点での事象の生起は独立であるとする。
としてシミュレーションを行うと、図１の結果を得た。
図１
横軸が時点を表し、縦軸が累積頻度を表す。曲線は、一様に右上がりではなく、急激に上がっている区間、すなわち、事象がまとまって生起している区間と、水平線の長い区間、したがって事象の生起のない区間が認められる。各時点における生起確率は
で一定であるが、事象はまとまって生起が観測される区間と、生起が（あまり）観測されない区間がある。図１のグラフは、リスト１に示すスクリプトで描いた。
リスト１
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng()
p = 0.1
xcoord = [0]
ycoord = [0]
x = 0
y = 0
for i in range(1000):
if rng.binomial(1, p) ==
1: # 事象が生起した
x +=
1
xcoord.append(x)
ycoord.append(y)
y +=
2
# y軸方向に上に2シフト
xcoord.append(x)
ycoord.append(y)
else:
x +=
1
xcoord.append(x)
ycoord.append(y)
plt.figure(figsize=(12,5))
plt.plot(xcoord, ycoord)
plt.show()
ある時点
から考えて、
時点の間事象が生起せず
時点において事象が生起する確率は、各時点における事象の生起確率をとおくと
で表される。これは、幾何分布と呼ばれている。
時点を連続量とすると指数分布が得られる。
時点0から時間後に事象が生起する確率を
とおく。
このとき、事象が生起する時点を
とおき、時点までは事象が生起しなかった場合の時点から
において事象の生起する確率は
となる。確率
は、に依存せずからの経過時間
と
に依存して決まる。
とおいて、リスト２のスクリプトでシミュレーションを行うと、図２の結果を得た。
図２
横軸に時間をとり、縦軸で事象の発生を示している。事象が発生すると上方向に曲線が移動している。上への移動が密な区間と、ほぼ横に移動している区間が認められる。上への移動が密な区間では、それに対応して事象が密に生起していることが示されているが、縦方向の変化がなく横に移動している区間では事象の生起のないことが表されている。
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
rng = np.random.default_rng()
b = 10.0
xcoord = [0]
ycoord = [0]
x = 0
y = 0
xcoord.append(x)
ycoord.append(y)
while xcoord[-1] < 1000:
x +=
rng.exponential(scale=b)
# 時間x後に事象が発生
xcoord.append(x)
ycoord.append(y)
y += 2
# 縦軸方向にy移動
xcoord.append(x)
ycoord.append(y)
plt.figure(figsize=(12,5))
plt.plot(xcoord, ycoord)
plt.show()
事象の生起時間間隔が式（１）で与えられるとき、単位時間内の生起事象の回数はパラメータ
のポアッソン分布で与えられる。パラメータのポアッソン分布は次式で表される。
生起間隔が式（１）で与えられる事象が単位時間内に生起する回数について、
の値に対してリスト３のスクリプトによりシミュレーションした結果を図３に示す。
図３
シミュレーションの結果をヒストグラムで示すとともに、対応するポアッソン分布を棒グラフで示した。ヒストグラムと棒グラフはよく一致している。
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
rng = np.random.default_rng()
b = 0.1 # rateパラメータ
counts = []
for i in range(1_000_000):
c = 0
t = 0
while True:
vt =
rng.exponential(b) # 次の事象生起までの時間
if t
+ vt > 1.0:
# 現在の事象生起は単位時間内か
break
c +=
1
t +=
vt
counts.append(c)
rng_pssn = np.arange(np.min(counts),
np.max(counts)+1)
plt.hist(counts, density=True,
bins=rng_pssn, align='left',
color='g', alpha=0.5, label='Exponential')
pssn = ss.poisson.pmf(rng_pssn,
1/b) # Lambda = 1/b
plt.bar(rng_pssn, pssn, color='r',
alpha=0.5, label='Poisson')
plt.legend()
plt.show()
参考文献
Hogg, R. V., McKean, J. W., & Craig, A. T. (2005) Introduction to Mathematical Statistics, 6th ed. Pearson Education, Inc.