ベイズの定理と言えば、赤玉白玉の問題が定番ですがせっかく勉強するので競馬を題材にやってみました。
ベイズの定理の公式
よく見かけるけど、よう分からんという式ですね。。
事象Aを観測したあとで事象Bが発生する確率は
と表すことができます。
事象Aに関するデータを観測することで、事象Bの確率が P(B)から P(B|A) へ、より尤もらしい値へと更新される ということを意味しています。
ここで、
- P(B): 事前確率
- P(B|A): 事後確率
- P(A|B): 尤度（ゆうど）
と呼びます。
P(A) は実際には B が取りうるすべての値に対しての P(A | B) * P(B) に関する積分であり、
となります。
P(A)と書くとシンプルに見えますが、実際の計算は大変です。。
その他、
- P(A, B) = P(A | B) * P(B) = P(B | A) * P(A): 同時確率
- P(A): 周辺尤度
と言うようです。
また、より学問的な世界では A, Bではなく DもしくはX（データ）, θ（分布パラメータ） という表現をするようです。
ただ、θ という見慣れない記号を見るだけで急に難しく感じ引いてしまうのは私だけではないでしょう。。
問題設定
では問題設定してみます。
問題
- 厩舎は成績によってAクラス、Bクラス、Cクラスに分類される
- 厩舎に所属する馬が2頭連続で勝利した場合、それらの馬がAクラス厩舎所属である確率は？
せっかくなのでリアルな数字でやってみたいと思います。
netkeibaの調教師リーディングのページから 厩舎ごとの勝利数を取得し、勝利数の分布をプロットしてみます。
厩舎の成績データを取得するコードはこんな感じです。（python2）
# encoding: utf-8
from bs4 import BeautifulSoup
import urllib2
import chardet
import traceback
import pickle
def fetchPage(url):
req = urllib2.Request(url)
res = urllib2.urlopen(req)
body = res.read()
guess_enc = chardet.detect(body)
try:
unicode_html = body.decode(guess_enc['encoding'])
except UnicodeDecodeError:
print url, guess_enc
print traceback.format_exc()
unicode_html = None
if not unicode_html:
try:
unicode_html = body.decode('shift_jisx0213')
print 'Decoding in shift_jisx0213 is successful.'
except UnicodeDecodeError:
print traceback.format_exc()
unicode_html = body
return unicode_html
def getTrainerResult(html, offset=0):
results = []
soup = BeautifulSoup(html, "html.parser")
trainer_leading_table = soup.select("table.race_table_01 tr")
for i, row in enumerate(trainer_leading_table[2:]):
# print row
data = row.select("td")
name = unicode(data[1].string)
stable = data[2].string
win_count = int(data[4].string)
second_place_count = int(data[5].string)
third_place_count = int(data[6].string)
unplaced_count = int(data[7].string)
grade_race_win = int(data[9].string)
stakes_race_win = int(data[11].string)
not_stakes_race_win = int(data[13].string)
if stable != u'栗東' and stable != u'美浦':
continue
result = {'no': i + offset,
'name': name,
'win_count': win_count,
'second_place_count': second_place_count,
'third_place_count': third_place_count,
'unplaced_count': unplaced_count,
'grade_race_win': grade_race_win,
'stakes_race_win': stakes_race_win,
'not_stakes_race_win': not_stakes_race_win}
results.append(result)
return results
leading_results = []
for i in range(4):
page = i + 1
html = fetchPage("http://db.netkeiba.com/?pid=trainer_leading&year=2017&page=%d" % page)
# print(html)
results = getTrainerResult(html, offset=len(leading_results))
# print(results)
leading_results += results
print(leading_results)
print(len(leading_results))
with open("trainer_leading.pkl", "w") as f:
pickle.dump(leading_results, f)
成績のヒストグラムをプロットすると次のようになりました。（X軸：勝利数、Y軸：頻度）
このヒストグラムを見て次のように3つのカテゴリに分けます。
- Aクラス: 25勝以上 （35厩舎）
- Bクラス: 11勝以上25勝未満 （110厩舎）
- Cクラス: 10勝以下 （51厩舎）
それぞれのクラスでの平均勝率は計算して以下のようになりました。
- Aクラスの勝率: 11.8%
- Bクラスの勝率: 6.7%
- Cクラスの勝率: 3.3%
求めるものは 事後確率 P(Stable=A | win,win) となります。
ベイズの定理を使って事後確率を計算する
事前確率を計算する
まず、事前確率を求めます。
厩舎の分布から計算します。
P(Stable=A) = 35 / (35 + 110 + 51) = 18% P(Stable=B) = 110 / (35 + 110 + 51) = 56% P(Stable=C) = 51 / (35 + 110 + 51) = 26%
尤度を計算する
次に、尤度を計算します。
各クラスでの勝率は
P(win | Stable=A) = 11.8% P(win | Stable=B) = 6.7% P(win | Stable=C) = 3.3%
なので、
2勝するという確率は
P(win,win | Stable=A) = 0.118 * 0.118 = 0.013923999999999999 P(win,win | Stable=B) = 0.067 * 0.067 = 0.004489000000000001 P(win,win | Stable=C) = 0.033 * 0.033 = 0.0010890000000000001
となります。
P(win, win)を計算する
とりうる厩舎クラスに対して同時確率を積分します。
P(win,win) = P(Stable=A) * P(win,win | Stable=A)
+ P(Stable=B) * P(win,win | Stable=B)
+ P(Stable=C) * P(win,win | Stable=C)
= 0.005303300000000001
事後確率を計算する
ベイズの定理に従って事後確率を計算します。
P(Stable=A | win,win) = P(win,win | Stable=A) * P(Stable=A) / P(win,win)
= 0.013923999999999999 * 0.18 / 0.005303300000000001
= 0.47259630795919505
約47%となりました。
何も情報がない状態では ある馬が Aクラスの厩舎の所属馬であるかどうかは 厩舎数の分布から見積もった 11.8% ということになりますが 勝利というデータを観測することによってその確率は 47% に高まります。
1勝の場合
2勝ではなく1勝だったときの事後確率 P(win | Stable=A) は
P(win | Stable=A) = 11.8% = 0.118
P(win | Stable=B) = 6.7% = 0.067
P(win | Stable=C) = 3.3% = 0.033
P(win) = P(Stable=A) * P(win | Stable=A) + P(Stable=B) * P(win | Stable=B) + P(Stable=C) * P(win | Stable=C)
= 0.06734000000000001
P(Stable=A | win) = P(win | Stable=A) * P(Stable=A) / P(win)
= 0.2988416988416988
となり、約30% となります。
勝利が続くことでAクラス厩舎である確率がより高まることが分かります。（30% → 47%）
分母は無視することが多い（らしい）
P(win, win) は 積分した結果であり P(Stable=A) を変数として含まないので分母は無視して計算することが多いようです。
実際に計算してみます。
分子のみを計算し、比率を求めます。
P(win,win | Stable=A) * P(Stable=A) = 0.013923999999999999 * 0.18 = 0.00250632
P(win,win | Stable=B) * P(Stable=B) = 0.004489000000000001 * 0.56 = 0.0025138400000000007
P(win,win | Stable=C) * P(Stable=C) = 0.0010890000000000001 * 0.26 = 0.00028314000000000003
P(Stable=A | win,win) = P(win,win | Stable=A) * P(Stable=A) /
(P(win,win | Stable=A) * P(Stable=A)
+ P(win,win | Stable=B) * P(Stable=B)
+ P(win,win | Stable=C) * P(Stable=C))
= 0.47259630795919505
分母も使って計算したときと同じ47%になりました!
