戦闘力のモデル化をpymcでやってみる。

pymc3のExampleには、ラグビーチームの攻撃力や防御力をモデル化する例がある。
https://docs.pymc.io/notebooks/rugby_analytics.html

こういったモデル化もあるのかと思っていると、世の中には、面白いことを考える人がいて、ドラゴンボールのキャラの戦闘力をモデル化した人がいる。

BUGS版（大元の記事）
http://statmodeling.hatenablog.com/entry/bugs-skill-performance

stan版
http://kazufusa1484.hatenablog.com/entry/2017/07/31/043000

BUGS版をベースにpymc3で実装してみた。大体大元の記事と同じ結果が得られたが、モデルの一部コードを少し変えないと、MCMCがエラーで実行できなかったので、記事に書くことにした。

まず大元記事と同じデータを作る。

# キャラのスペック
data = {
    'CID': [1,2,3,4,5],
    'name':['ベジータ','悟空','ナッパ','クリリン','栽培マン'],
    'power.level':[18000,8000,4000,1770,1200]
}
senshi = pd.DataFrame.from_dict(data)

# 試合結果
data = {
    'winner':[1,1,1,2,3,4],
    'loser':[2,3,5,3,4,5]
}
battles = pd.DataFrame.from_dict(data)

BUGSのstep関数の代わりにswitch関数を使ってみたが、これだとmcmcがエラーで落ちたので、 sigmoid関数を使った。あとは同じ（と思う）。

with pm.Model() as model:
    # prior
    skill = pm.Normal('skill', mu=0, sd=100, shape=len(senshi.CID))
    
    winner_perf = pm.Normal('winner_perf', mu=skill[battles.winner-1], sd=1, shape=len(battles))
    loser_perf = pm.Normal('loser_perf', mu=skill[battles.loser-1], sd=1, shape=len(battles))
    
    diff = winner_perf-loser_perf
    #p = pm.math.switch(diff, 1.0, 0.0)
    p = pm.math.sigmoid(diff)
    z = pm.Bernoulli('z',p=p, observed=np.ones(len(battles)))

with model:
    trace = pm.sample(1000, start=pm.find_MAP(), step=pm.NUTS())

logp = -38.675, ||grad|| = 0.00058994: 100%|██████████| 84/84 [00:00<00:00, 643.83it/s]  
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [loser_perf, winner_perf, skill]
Sampling 2 chains: 100%|██████████| 3000/3000 [03:08<00:00, 15.91draws/s]
The estimated number of effective samples is smaller than 200 for some parameters.

pm.traceplot(trace)

f:id:gaiasky:20180825194111p:plain

# 戦闘力部分の推定結果を抽出
result_df = pm.summary(trace, varnames=['skill'])

# パラメータ名をキャラ名に変換
senshi_dict = senshi.to_dict()['name']
senshi_dict2 = {}
for id,name in senshi_dict.items():
    senshi_dict2.setdefault('skill__%s'%(id), name)

result_df.rename(index=senshi_dict2, inplace=True)

真の戦闘力と推定戦闘力をviolinplotで表示する。

# violinplotに入れるデータの生成
data = []
for name in result_df.index:
    mu = result_df.loc[name, 'mean']
    sd = result_df.loc[name, 'sd']
    d = pm.Normal.dist(mu=mu,sd=sd).random(size=100)
    tmp = pd.DataFrame(data=d, columns=['power'])
    tmp['name'] = name
    ap = np.log10(senshi[senshi['name']==name]['power.level'])
    tmp['actual_power'] = np.round(ap.values[0], decimals=2)
    data.append(tmp)
data_df = pd.concat(data)

plt.figure(figsize=(12,9))
sns.violinplot(x='actual_power', y='power', data=data_df, hue='name')