Kei Minagawa's Blog

皆川圭(@keimina)のブログ、Pythonで試したことを書いていきます

三面体のサイコロを2000回投げたとき、それぞれの面の出る回数の分布

1. はじめに

(注意:一般的に三面体のサイコロというものは現実には存在しません。説明の便宜上、3つの面をもつ空想のサイコロのことを三面体のサイコロと表現しています)
Pandas 勉強会で、サンプルの状態が割とよく発生することなのか、稀にしか発生しないことなのかを知りたいと言う話がありました。それに関連するような、確率の問題を自分で作って、シミュレーションを行いました。結果、仮定がデータと紐づいているのであれば、割とよく発生することなのか、稀にしか発生しないことなのかがわかることがわかりました。一方で仮定がデータと紐づかないのであれば、仮定に基づいて議論を進めたところであまり意味のない議論になってしまうこともわかりました。

2. 問題

以下のように問題を定義しました。

2000個のサンプルがあったとして、 3つの領域 A, B, C の内どれかに分類されるとした時、観測したそのサンプルが領域 A に分類された個数、 B に分類された個数、 C に分類された個数の状態が、割とよくおこるものなのか、そうでないのか知りたい。

3. 問題を解くための仮定

まず、Aに分類されるか、B または C に分類されるかは一定の確率で決まることとします。議論を進めるため、その確率を以下のように定義します。

A の出る確率: 1/3
B の出る確率: 1/3
C の出る確率: 1/3

これは、簡単に言えば三面体のサイコロのそれぞれの面がでる確率と同じです。サイコロの方が直感的にイメージしやすいため、問題を3面体のサイコロをを2000回投げた時、それぞれの面のでる回数の分布を求める問題に置き換えます。A に分類されるか、B または C に分類されるかの確率が決まっていますので、この確率を前提として、シミュレーションを行います。

4. シミュレーションを行うコード

シミュレーションを行うコードは以下になります。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy.stats

# 1. 
samples = np.random.choice(['A', 'B', 'C'], p=[1/3, 1/3, 1/3], size=(10000, 2000))
# 2.
df = pd.DataFrame(samples)
# 3.
df = df.T.apply(lambda x: x.value_counts())
# 4.
plt.hist(df.loc['A', :], bins=50)
plt.show()

# 5.
mean = df.loc['A', :].mean()
sigma = df.loc['A', :].std()
print('平均:{}回\n標準偏差:{}'.format(mean, sigma))

# 6.
n, bins, patches = plt.hist(df.loc['A', :], bins=50, density=True) 
# 7.
y = scipy.stats.norm.pdf(loc=mean, scale=sigma, x=bins)
plt.plot(bins, y)
plt.show()

# 8.
print("'A'のでる回数xは...")
print("68% は {:.0f} < x < {:.0f}となる".format(mean - sigma, mean + sigma))
print("95% は {:.0f} < x < {:.0f}となる".format(mean - 2*sigma, mean + 2*sigma))
print("99.7% は{:.0f} < x < {:.0f}となる".format(mean - 3*sigma, mean + 3*sigma))

5. コードの説明

コードの説明は以下の通りです。

1. まず、 `np.random.choice` 関数で 'A', 'B', 'C' をそれぞれ 1/3 の確率で 2000個 出現する配列を作成します。 第2引数は第1引数で指定したリストの要素が出現する確率を指定しています 。今回はさらに、 それを10000回行ったこととし、 size 引数に (10000, 2000) を渡すことで 10000行 の行列を作成します。
2. 次に 1. で作成した numpy配列を DataFrame に変換します。
3. df を転置させたものの全ての列に対して、 `Series.value_counts` メソッドを適用します。
4. 3. で求めたものをヒストグラムで可視化します。形状が正規分布となっていることがわかります。
5. 平均と標準偏差を求めます。
6. `scipy.stats.norm.pdf` 関数のパラメータに先ほどの平均と標準偏差と 3. で表示したヒストグラムのバーのx座標の値を渡し、その母数での正規分布のx座標に対するy座標の値を得ます。
7. 3. で表示したヒストグラムを面積が1となるように正規化したグラフと、6. で求めた正規分布のグラフを重ねて表示します。
8. 正規分布の性質(平均値から±標準偏差の3倍までの面積は 0.997 、2倍までの面積は 0.95 、1倍までの面積は 0.68) から、範囲に対する確率を求め、表示します。

6. 出力結果

出力結果は以下になります。

平均:666.8196回
標準偏差:21.023074294814155
'A'のでる回数xは...
68% は 646 < x < 688となる
95% は 625 < x < 709となる
99.7% は604 < x < 730となる

f:id:keimina:20191005012451p:plain
ヒストグラム1

f:id:keimina:20191005012533p:plain
ヒストグラム1(正規化済)&正規分布

7. 考察・まとめ

問題を定義し、仮定を定義しました。三面体のサイコロ投げで面が出る回数について10000回シミュレーションを行い、その確率の分布を求めました。サイコロの目が出る回数を数えたら結局は正規分布になりました。今回行ったシミュレーションから 99.7% の割合で 'A' が出現する回数xは 604 < x < 730 の範囲となると推測しました。 この範囲外にあった場合、それは割とよく発生するものでない(異常な状態)ということがわかります。ただし、この推測が役に立つのは、あくまで仮定が成立する時の話で、この仮定とデータが何ら紐づかない場合は実際あまり意味のない議論になります。しかし、データと仮定が紐づかない場合でも、このような条件が議論する上で重要であるということがわかっているのであれば、何の指標も方向性も定まっていない状態より良いのではないかと思います。

以上です。