二項分布のシミュレータを作った (Python3)
内容
1列6席の座席に等確率で着席する様子をシミュレーションする
まず座席に座る様子を
n席に着席する様子をM回実行するというもの。
numpyのrandom.randint(2, size=n)
というのは、0と1のintの乱数をsize=nでリストを生成するという関数。
とっても便利だった。forでサイズ分回すという必要がないからコードも読みやすくなっているはず。
あとは、 print_seat()
。これは、一列6席というのを守るために長いリストを要素数6に分割するというもの。
もちろんリストの最後の余り(席数が10とかだと二列目は4席だけ)も吸収できる。
import random import itertools import numpy as np def M_loop(n, M): for i in range(M): seat = simu1(n) print_seat(seat) def simu1(n): occupied = np.random.randint(2, size=n) return occupied def print_seat(seat_list): for a in itertools.zip_longest(*[iter(seat_list)]*6): print(a) print("=" * 40) if __name__ == '__main__': M_loop(10, 100)
出力結果はこんな感じ
次は各座席に座る確率
n=6の時に書く座席に座る確率を求める。あまり変わりはないが、座席のリストに1ずつ加算していき最後はループ回数で除算する。すると確率は出てくる。これはもちろん大凡0.5となる。
import random import itertools import numpy as np def M_loop(n, M): seat_all = np.zeros(6) for i in range(M): seat = simu1(n) seat_all += seat return seat_all / M def simu1(n): occupied = np.random.randint(2, size=n) return occupied def print_seat(seat_list): for a in itertools.zip_longest(*[iter(seat_list)]*6): print(a) print("=" * 40) if __name__ == '__main__': result = M_loop(6, 10000) print(result)
結果はこんな感じ。まああってるでしょ。
何席座るかを確率で表現する
6席ある中で、何席座る確率を求める。3席座る確率は、とか6席座る確率は?とかそんなもの。今度は、どの座席に座るかでなくて、何席座ったかを集めている。count[np.sum(simu1(n))]
で[0. 1. 1. 0. 1. 1.]という座席の座り方をしたときに、sumで4と返す。そこでcountの4番目のリストに1を加算して回数をカウントしている。
import random import itertools import numpy as np import matplotlib.pyplot as plt def M_loop(n, M): count = np.zeros(7) for i in range(M): count[np.sum(simu1(n))] += 1 return count def simu1(n): occupied = np.random.randint(2, size=n) return occupied def print_seat(seat_list): for a in itertools.zip_longest(*[iter(seat_list)]*6): print(a) print("=" * 40) if __name__ == '__main__': n = 6 M = 100000000 result = M_loop(n, M) p_k = result / M print(p_k) x = range(7) y = p_k plt.xlabel("seats") plt.ylabel("probability") plt.plot(x, y) plt.show()
平均と分散を求める
シミュレーションした結果が妥当かどうかを確認する方法として平均と分散を求める。この場合の平均は期待値と同義。 さっきのコードとほとんど変わらないんだけど、要するに期待値を計算している。np.dotで行列の内積を求めている。dotを使うのが一番速いらしい。どこかのページで見た。。。。分散は、np.varで求めれば十分。
import random import itertools import numpy as np def M_loop(n, M): count = np.zeros(7) for i in range(M): count[np.sum(simu1(n))] += 1 return count / M def simu1(n): occupied = np.random.randint(2, size=n) return occupied def print_seat(seat_list): for a in itertools.zip_longest(*[iter(seat_list)]*6): print(a) print("=" * 40) if __name__ == '__main__': n = 6 M = 100000000 result = M_loop(n, M) x = np.arange(7) mean = np.dot(result, x) print(mean) print(np.var(result))
ポアソン分布もやってみるか
生起確率を0.02としてポアソン分布を試してみる。これなら、a = npですぐ確認できるから楽。 ほとんど変わらないが、コインフリップで片方の生起確率を下げるとかどうするんだろうと思っていたところ、np.random.random()が(0,1]であることがわかる。するとこれにif文かませばいけるな、ということでできた。
import random import itertools import numpy as np import matplotlib.pyplot as plt from scipy.stats import norm def M_loop(n, M): count = np.zeros(n+1) for i in range(M): seat = simu1(n) count[np.sum(simu1(n))] += 1 return count / M def coin_flip(p): return 1 if np.random.random() < p else 0 def simu1(n): occupied = np.zeros(n+1) for i in range(n+1): occupied[i] += coin_flip(0.02) return occupied def print_seat(seat_list): for a in itertools.zip_longest(*[iter(seat_list)]*6): print(a) print("=" * 40) def poisson_curve(n, M): poi_sample = np.random.poisson(2, n+1) param = norm.fit(poi_sample) return param if __name__ == '__main__': n = 100 M = 100000000 result = M_loop(n, M) print(result) x = np.arange(n+1) mean = np.dot(result, x) print(mean) print(np.var(result)) param = poisson_curve(n, M) print(param) fitted = norm.pdf(x,loc=param[0], scale=param[1]) plt.plot(x,fitted, 'b--', ) plt.plot(x,result, 'r') plt.xlabel("seats") plt.ylabel("probability") plt.show()
numpyにはポアソン分布に従った乱数発生はあるんだけど、分布を出すものはなかった(見つけられなかった)。そこで近似曲線を引いている。 param = norm.fit(poi_sample)
で平均と分散をポアソン分布に従った乱数から求めている。そしてそれを描画しているという訳。
青の点線が理論曲線で赤が計算値。ほぼほぼ同じになってる。