強弱と高低

Python3をメインに

二項分布のシミュレータを作った (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)

出力結果はこんな感じ f:id:Accent:20160529140026p:plain

次は各座席に座る確率

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)

結果はこんな感じ。まああってるでしょ。 f:id:Accent:20160529140030p:plain

何席座るかを確率で表現する

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) で平均と分散をポアソン分布に従った乱数から求めている。そしてそれを描画しているという訳。

青の点線が理論曲線で赤が計算値。ほぼほぼ同じになってる。 f:id:Accent:20160529140034p:plain