第４回 データの分布と確率モデル
プログラム（Pythonスクリプト）の完全なリストを以下に上げる。スクリプトのファイルは圧縮ファイルfiles4.zipとしてまとめたので、クリックしてダウンロード解凍すればスクリプトファイルを開くことができる。プログラムの説明は連載本文を参照されたい。
リスト１ 幾何分布と指数分布
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
p = 0.1
for t in range(100):
if np.random.uniform() <
p:
print('1', end = '')
else:
print('0', end = '')
print()
def w_time():
t = 0
while True:
if
np.random.uniform() < p:
break
t +=
1
return t
for i in range(100):
print(w_time(), end = ' ')
n_samples = 10000
p = 0.1
data = []
for i in range(n_samples):
data.append(w_time())
plt.hist(data, density = True, bins =
20)
def exp_d(t, b):
return b * np.exp(-t * b)
b = p / (1.0 - p)
max_data = max(data)
t_range = [t for t in np.arange(0,
max_data, 0.1)]
probs = [exp_d(t, b) for t in t_range]
plt.plot(t_range, probs, linewidth =
3)
plt.show()
リスト２ 指数分布とポアッソン分布
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
def my_poisson(b):
c = 0
v =
np.random.exponential(1/b)
while True:
if v
> 1:
break
v +=
np.random.exponential(1/b)
c +=
1
return c
b = 10
data = []
for i in range(10000):
data.append(my_poisson(b))
data = np.array(data)
max_v = data.max()
x_rng = np.arange(-0.5, max_v + 1,
1.0)
plt.hist(data, density = True, bins =
x_rng)
x_range = np.arange(0, max_v)
ck_pssn = [scipy.stats.poisson.pmf(v,
b) for v in x_range]
plt.plot(x_range, ck_pssn, linewidth =
3)
plt.show()
リスト３ カテゴリカル分布
import numpy as np
import matplotlib.pyplot as plt
v = np.random.choice([1,2,3,4,5,6], p
= [1/6, 1/6, 1/6, 1/6, 1/6, 1/6], size = 10)
print(v)
v = np.random.choice(['A', 'B'], p =
[2/3, 1/3], size = 10)
print(v)
data = np.random.choice([1,2,3,4,5,6],
p = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6], size = 10000)
xindx = [1,2,3,4,5,6]
plt.hist(data, bins = [0.5 + v for v
in range(7)])
plt.show()
リスト４ 多項分布
import numpy as np
import matplotlib.pyplot as plt
v = np.random.multinomial(3, [1/6,
1/6, 1/6, 1/6, 1/6, 1/6])
print(v)
v = np.random.multinomial(3, [1/6,
1/6, 1/6, 1/6, 1/6, 1/6], size = 10)
print(v)
data = np.random.multinomial(3, [1/6,
1/6, 1/6, 1/6, 1/6, 1/6], size = 10000)
freqs = data.sum(axis = 0)
print(freqs)
plt.bar(range(1,7), freqs)
plt.show()