離散値の確率を計算するのに便利な二項分布を、Pythonを使って実装してみようという趣旨。
離散型の確率分布として使われる二項分布は、例えば勝率が3割のチームが7戦した時に勝つ確率を求めるときなどに役立ちます。
この記事はこんな人におすすめ。
- 勝敗の確率を求めるシミュレーションの方法を知りたい。
- 二項分布を使って勝敗の確率を知りたい。
今回は、勝率2割5分のチームが6試合で3勝する確率を求めてみましょう。
勝敗の確率を求めるシミュレーション
ではまずシミュレーションで、6試合を4回行った場合の勝敗は以下のような結果になります。
ランダムで0-1の値を出し、0.25以下の場合勝ち、それ以上の場合負けとします。
# 勝率25%で6試合を4回行った時の勝敗を求める。
import numpy as np
np.random.random([4, 6])<0.25 # 勝率25%で6試合を4回行った時の勝敗の結果
array([[False, False, False, True, False, True],
[False, False, False, True, False, True],
[False, True, True, True, False, False],
[ True, False, False, False, True, False]])
勝ちをTrue、負けFalse で表した結果を見てみると、6試合を4回行って3勝できるのは1回のみという結果がでます。
以下の記述で、6回試合をして4回行って3勝するのは約4分の1ぐらい、ということがわかります。
(np.random.random([4, 6])<0.25).sum(axis=1)==3 # 勝率25%で6試合を4回行って3勝する割合。
array([False, False, False, True])
次に勝率25%で6試合を100万回行ったときに3勝する確率を求めます。
n_traial = 1000000
n_battle = 6
win3 = ((np.random.random([n_traial, n_battle])<0.25).sum(axis=1)==int(n_battle/2)).sum()
win3/n_traial
0.131671
結果、勝率は0.132%と出ました。
二項分布の公式を使ってみる
では、同じことを2項分布の公式を使って計算してみます。
二項分布の公式は、以下のようになります。
そしてこの の部分は以下のように求められます。
まず 、つまり6試合行って3勝する組み合わせが、何通りあるかを計算します。
式に書くと、
これを実装すると以下の結果。
import math
n= 6
k= 3
math.factorial(n)/(math.factorial(k) * math.factorial(n - k)) # 標準ライブラリmath.factorialを使う。
20.0
つまり、6試合して3勝するパターンが20通りあるということですね。
では実際に の次の部分も含めて二項分布の公式を実装し、勝率25%で6試合行った場合、3勝する確率を求めてみます。
def binominal_dist(k, n, p):
c = (math.factorial(n)/(math.factorial(k) * math.factorial(n - k)))*((p**k)*(1-p)**(n-k))
return c
win6 = binominal_dist(3, 6, 0.25)
import pandas as pd
pd.DataFrame([[win6]], index=['6試合'], columns = ['勝率'])
勝率
6試合 0.131836
実行結果は、0.132%という結果が出ました。
かなり低いですね。
scipy.stats を使って実装する
同じことを、 scipy.stats を使って計算してみましょう。
scipy.statsを使うと3行で済むので、実用的ではありますね。
同じように勝率25%で、6試合して3勝する確率を求めてみます。
from scipy.stats import binom
win6 = binom.pmf(3,6,0.25)
勝率
6試合 0.131836
結果はやはり、0.132%という結果になっています。
二項分布の応用編
応用ということもないのですが、勝率25%で6試合して1勝する確率を求めてみます。
まず、6試合行って1勝もできない確率を求めましょう。
シミュレーションで、6試合を100万回試合をして1勝もできない確率を求めます。
n_traial = 1000000
n_battle = 6
win3 = ((np.random.random([n_traial, n_battle])<0.25).sum(axis=1)==int(0)).sum()
print(win3/n_traial)
0.17807
一勝もできない確率が0.178% なので、少なくとも1勝できる確率は 1-0.178% になりますね。
n_traial = 1000000
n_battle = 6
win3 = ((np.random.random([n_traial, n_battle])<0.25).sum(axis=1)==int(0)).sum()
print(1-(win3/n_traial))
0.822262
シミュレーションを実行してみると、必ず1勝できる確率は0.822%。
二項分布の公式を使って計算してみましょう。
import math
n= 6
k= 0
math.factorial(n)/(math.factorial(k) * math.factorial(n - k))
def binominal_dist(k, n, p):
c = (math.factorial(n)/(math.factorial(k) * math.factorial(n - k)))*((p**k)*(1-p)**(n-k))
return c
win6 = binominal_dist(0, 6, 0.25)
import pandas as pd
print(1-(pd.DataFrame([[win6]], index=['6番勝負'], columns = ['勝率'])))
勝率
6番勝負 0.822021
やはり結果は、0.822%。
念のためscipy.statsを使っても計算してみましょう。
from scipy.stats import binom
win6 = binom.pmf(0,6,0.25)
print(1-win6)
0.822021484375
やはり結果は0.822%となり、勝率25%のチームが6試合で少なくとも1勝する確率は0.822%という結果がでました。
コメント