2022年12月2日金曜日

サイコロ振って中心極限定理を確かめる

 中心極限定理とは、

中心極限定理は標本平均と母平均との誤差を論ずるものである。多くの場合、母集団の分布がどんな分布であっても、その誤差は標本の大きさを大きくしたとき近似的に正規分布に従う。wikipedia

とあるように標本を大きくしていくと正規分布に近づくという定理です。

サイコロを想定してpythonで表現してみます。

numpyのrandintでランダムに1〜6の目を出すサイコロを回数分出します。

とりあえず10000回くらい振ってみましょう。

#! /usr/bin/env python

# -*- coding:utf-8 -*-

import pandas as pd
import numpy as np


print("[dice]")
print()

# サイコロの目(1〜6)をランダムにnumpy配列で出力
pip = np.random.randint(1, 7, 10000)

# 作られた配列の先頭を表示
print(pip[:10])
print()
print(f"平均値: %f" % pip.mean())
print()

そしてmatplotlibのヒストグラムで分布を見てみましょう。

#ヒストグラムでプロット
import matplotlib.pyplot as plt
import japanize_matplotlib

plt.title("サイコロ")
plt.xlabel("出た目")
plt.ylabel("出た回数")
plt.hist(pip, bins=6, ec="k")
plt.show()

平均値: 3.502500
1~6までほぼ同じ回数出てるのが分かります。
期待値は(1+2+3+4+5+6) / 6 = 7/2 = 3.5なので平均値もほぼ合ってます。
サイコロ1つなら感覚的にもつかみやすいでしょう。

次に2個のサイコロを振ってみます。



import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
#最瀕値を出すために使用
import statistics

#振る回数
roll = 10000

#1つ目のサイコロ
pip = np.random.randint(1, 7, roll)
#2こ目のサイコロの
pip2 = np.random.randint(1, 7, roll)

print("[dice]")
print()
#先頭10回分の目
print(pip[:10])
print(pip2[:10])

sm = pip + pip2
print()
print(f"平均値: %f" % sm.mean())
print(f'最瀕値: %f' % statistics.mode(sm))

plt.title("サイコロ")
plt.xlabel("出た目")
plt.ylabel("出た回数")
plt.hist(sm, bins=11, ec="k")

plt.show()

平均値: 7.044900

最瀕値: 7.000000


サイコロを2つにすると足した7が一番多い山のような形になることが分かります。




振るサイコロ6個の目の合計ではどうでしょう?

サイコロ6つだと、出る目は6の6乗の46656とおりになるので、それ以上回数を振らないとまばらな結果になりそうです。

10万回くらい振ってみましょう

#! /usr/bin/env python

# -*- coding:utf-8 -*-

#
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
#最瀕値を出すために使用
import statistics

#振る回数
roll = 100000

# ndarrayの初期化
sm = np.zeros(roll)

print("[dice]")


for i in range(6):
    pip = np.random.randint(1, 7, roll)
    
    #出ための先頭10回
    print(pip[:10])
    sm += pip

print()
print(f"平均値: %f" % sm.mean())
print(f"最瀕値: %f" % statistics.mode(sm))
print()
#出ための合計の種類
print(np.unique(sm))


plt.figure(figsize=(10, 8))
plt.title("サイコロ")
plt.xlabel("出た目の合計")
plt.ylabel("出た回数")
plt.grid()
plt.hist(sm, bins=31, ec="k")

plt.show()

平均値: 21.013420

最瀕値: 21.000000



ほぼ正規分布に近づきました。

標本を大きくすれば正規分布に近づいていくという中心極限定理が感覚的につかめました。