TECH MEMO

技術的な内容に関する個人的なメモ

PyMCによるベイズ推論 〜世界カーリング選手権における先行/後攻別の勝率の分析〜

前回の記事で世界カーリング選手権における先行/後攻別の勝率を算出し、勝率の差が誤差ではないか確かめるために、検定(二項検定、χ2検定)を実施しました。

上記に対して、今回はより直感的な理解ができるベイズ推論によるアプローチで、勝率に差が存在するのか確かめていきたいと思います。
今回も対象はいつもと同様に2002年〜2016年までの15年間の男女世界選手権のデータです。
変数result_score_df(pandasのDataFrame型)に以下のような形式でデータが格納されている前提で進めます。

f:id:curlyst:20170729184315p:plain

なお、今回の内容を実施するにあたって、以下の書籍を参考にしました。

Pythonで体験するベイズ推論 PyMCによるMCMC入門

Pythonで体験するベイズ推論 PyMCによるMCMC入門

  • 作者: キャメロンデビッドソン=ピロン,玉木徹
  • 出版社/メーカー: 森北出版
  • 発売日: 2017/04/06
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログを見る

1. 先行/後攻別の勝率

まずはデータ全体における先行/後攻での勝率の差ですが、前回記事では二項検定を実施して、後攻の勝率が二項分布(p=0.5)から有意に偏っているか確かめました。その結果、偏りが存在し、後攻の勝率はp=0.5とは異なる(先行と後攻で勝率は異なる)となりました。

上記についてのベイズ推論によるアプローチでは、まず後攻の勝率に関する事前分布をモデリングします。勝率については何も確信がないので、ここでは一様分布に従うと仮定します。

# 後攻(hammer)チームの勝率(事前分布)
p_hammer = pm.Uniform('p_hammer', lower=0, upper=1)

次に実際の勝ち負けのデータ(観測データ)をモデルに組み込みます。試合の勝ち負けに関する確率変数には、二値変数についての確率分布であるベルヌーイ分布*1を用います。

# 観測データ
hammer_result_list = result_score_df[result_score_df["hammer"] == 1]["result"].tolist()
hammer_result_list = [True if x == 1 else False for x in hammer_result_list]

# 観測データを設定
obs = pm.Bernoulli("obs", p_hammer, value=hammer_result_list, observed=True)

これでモデルの設定は完了したので、以下を実行し、MCMC*2により勝率に関する事後分布から、データをサンプリングします。ここではサンプリングを10000回、バーンイン期間を1000としました。

mcmc = pm.MCMC([p_hammer, obs])
mcmc.sample(10000, 1000)

実行結果

[-----------------100%-----------------] 10000 of 10000 complete in 1.1 sec

サンプリングされた勝率のヒストグラムを生成した結果を以下に示します。赤色の縦の点線は0.5(勝率5割)を示しており、大部分のデータが勝率0.5以上に分布している(先行に比べて後攻の方が勝率が良い)ことがわかります。

f:id:curlyst:20170730200443p:plain

次に勝率がx割を上回る確率がどの程度かを算出してみます。

# 5割5分を超える確率
print("p > 0.55: %.2f" % (mcmc.trace("p_hammer")[:] - 0.55 >= 0).mean())
# 5割6分を超える確率
print("p > 0.56: %.2f" % (mcmc.trace("p_hammer")[:] - 0.56 >= 0).mean())
# 5割7分を超える確率
print("p > 0.57: %.2f" % (mcmc.trace("p_hammer")[:] - 0.57 >= 0).mean())
# 5割8分を超える確率
print("p > 0.58: %.2f" % (mcmc.trace("p_hammer")[:] - 0.58 >= 0).mean())

実行結果

p > 0.55: 0.94
p > 0.56: 0.77
p > 0.57: 0.43
p > 0.58: 0.16

後攻チームが勝率5割5分以上になる確率が94%となり、先行チームに比べて後攻チームが有利であることに対する確信度はかなり高いと考えられます。

2. 日本、カナダにおける先行/後攻別の勝率

それでは次に国別の先行/後攻における勝率について、日本とカナダを対象に調べていきます。前回の記事では、χ2検定を実施した結果、日本は先行/後攻が勝敗と関連はしておらず、カナダは関連しているという結果になりました。

ベイズ推論によるアプローチでは、先行での勝率、後攻での勝率、後攻での勝率ー後攻での勝率、の3つについて事後分布を求めます。先ほどと同様に、まず各勝率を一様分布によってモデリングします。

#--- 日本 ---
p_hammer_jp = pm.Uniform('p_hammer_jp', lower=0, upper=1)
p_not_hammer_jp = pm.Uniform('p_not_hammer_jp', lower=0, upper=1)

#--- カナダ ---
p_hammer_ca = pm.Uniform('p_hammer_ca', lower=0, upper=1)
p_not_hammer_ca = pm.Uniform('p_not_hammer_ca', lower=0, upper=1)

次に実際の勝敗データをモデルに組み込みます。こちらも先ほどと同様にベルヌーイ分布を用います。

# 日本
result_score_jp_df = result_score_df[result_score_df["team"] == "Japan"]
hammer_result_jp_list = result_score_jp_df[result_score_jp_df["hammer"] == 1]["result"].tolist()
hammer_result_jp_list = [True if x == 1 else False for x in hammer_result_jp_list]
not_hammer_result_jp_list = result_score_jp_df[result_score_jp_df["hammer"] == 0]["result"].tolist()
not_hammer_result_jp_list = [True if x == 1 else False for x in not_hammer_result_jp_list]

# カナダ
result_score_ca_df = result_score_df[result_score_df["team"] == "Canada"]
hammer_result_ca_list = result_score_ca_df[result_score_ca_df["hammer"] == 1]["result"].tolist()
hammer_result_ca_list = [True if x == 1 else False for x in hammer_result_ca_list]
not_hammer_result_ca_list = result_score_ca_df[result_score_ca_df["hammer"] == 0]["result"].tolist()
not_hammer_result_ca_list = [True if x == 1 else False for x in not_hammer_result_ca_list]

国別の勝敗に関するベイズ推論によるアプローチでは、先行/後攻の勝率だけでなく、それらの差分(後攻の勝率ー先行の勝率)の事後分布も推定するため、差分に関する確率変数を定義します。

# 日本
@pm.deterministic
def p_diff_jp(p_hammer_jp=p_hammer_jp, p_not_hammer_jp=p_not_hammer_jp):
    return p_hammer_jp - p_not_hammer_jp

# カナダ
@pm.deterministic
def p_diff_ca(p_hammer_ca=p_hammer_ca, p_not_hammer_ca=p_not_hammer_ca):
    return p_hammer_ca - p_not_hammer_ca

最後にMCMCにより事後分布からデータをサンプリングします。

# 日本
mcmc_jp = pm.MCMC([p_hammer_jp, p_not_hammer_jp, obs_hammer_jp, obs_not_hammer_jp, p_diff_jp])
mcmc_jp.sample(10000, 1000)

# カナダ
mcmc_ca = pm.MCMC([p_hammer_ca, p_not_hammer_ca, obs_hammer_ca, obs_not_hammer_ca, p_diff_ca])
mcmc_ca.sample(10000, 1000)

それでは日本、カナダそれぞれの勝率に関する事後分布を見ていきます。

まず日本ですが、後攻の勝率、先行の勝率、後攻と先行の差分、のそれぞれに関するヒストグラムを作成すると以下のようになります。先行と後攻で分布にあまり違いがなく、分布の範囲も重複している部分が多いため、差分についても0の周辺に分布した形になっています。先行の方が勝率が良いケース(差分が負)も多く存在し、日本は先行と後攻で差があまりないことが分かります。

f:id:curlyst:20170804164806p:plain

次にカナダにですが。カナダの勝率に関するヒストグラムを作成すると、以下のように後攻の勝率と先行の勝率の分布の範囲がずれており、後攻の方が0.1ほど大きくなっています。そのため、差分についても大半が0より大きい部分で分布しており、後攻の勝率の方が良い可能性が高いことを示しています。

f:id:curlyst:20170804165822p:plain

3. まとめ

今回は前回と同じデータを対象に、前回は検定によって分析していたところを、ベイズ推論によるアプローチで分析してみました。分析結果(結論)についてはどちらの方法でも変わらないですが、ベイズ推論によるアプローチでは、勝率の分布が実際にグラフとして可視化できるため、直感的に分かりやすいのではないかと思います。

ただ、難点として、ベイズ統計を知らない人にとっては「勝率の分布」とか言われても、中々理解するのは厳しいですかね...?