アヒル本10章棋士の強さのモデル化をpymcでやってみる。

この章の例題をpymcで実装してみた。 pymcによる実装の課題は以下の通り。

順序付けを行うorderedが、pymcにはない（と思われる）
単純にMCMCを実行するとエラーで実行できない

1点目は、以前のドラゴンボールの戦闘力推定の実装を使った。

gaiasky.hatenablog.com

# データの読み込み
players = pd.read_csv('shogi/data-shogi-player-name.txt', sep='\t')
battles = pd.read_csv('shogi/data-shogi-player.txt', sep=',')

モデル化は、ドラゴンボールの戦闘力推定モデルとほぼ一緒。

with pm.Model() as model:
    # prior
    s_mu = pm.Normal('s_mu', mu=0, sd=100)
    mu = pm.Normal('mu', mu=0, sd=s_mu, shape=len(players))
    sigma = pm.Gamma('sigma', 10, 10, shape=len(players))

    
    loser_perf = pm.Normal('loser_perf', mu=mu[battles.Loser-1], sd=sigma[battles.Loser-1], shape=len(battles))
    winner_perf = pm.Normal('winner_perf', mu=mu[battles.Winner-1], sd=sigma[battles.Winner-1], shape=len(battles))
    
    diff = winner_perf-loser_perf
    p = pm.math.sigmoid(diff)
    z = pm.Bernoulli('z',p=p, observed=np.ones(len(battles)))

MCMCの実行では、ADVI()による初期化を行うことで、NUTSによる実行が可能となった。（ADVIによる初期化をしないと、エラーで落ちる）

with model:
    trace = pm.sample(2000, init='advi', n_init=25000)

Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average Loss = 4,742.8: 100%|██████████| 25000/25000 [05:40<00:00, 73.47it/s]
Finished [100%]: Average Loss = 4,743.1
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [winner_perf, loser_perf, sigma, mu, s_mu]
Sampling 2 chains: 100%|██████████| 5000/5000 [44:36<00:00,  2.09draws/s]
The gelman-rubin statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.

pm.traceplot(trace, varnames=['mu','sigma'])

f:id:gaiasky:20180914230747p:plain

推定結果を表示する際、indexを棋士名に変換する。

def rename_indexes(trace, varname, players):
    ret = pm.summary(trace, varnames=[varname])
    
    tab = {}
    for p in range(len(players)):
        tab.setdefault('%s__%d'%(varname,p), players[players.nid==p+1].kname.values[0])
    ret = ret.rename(index=tab)
    return ret

各棋士の強さの平均を表示する。

ret_mu = rename_indexes(trace, 'mu', players)
ret_mu.sort_values('mean', ascending=False).head()

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
羽生善治	2.959619	0.268438	0.015966	2.412412	3.475476	31.911283	1.073706
渡辺明	2.538573	0.294463	0.017209	1.984500	3.114538	80.892714	1.030923
行方尚史	2.116030	0.313367	0.017155	1.464525	2.698630	223.610725	1.019783
豊島将之	2.067839	0.286370	0.015999	1.542303	2.653113	117.756783	1.045872
郷田真隆	1.995864	0.285740	0.016104	1.472823	2.532897	195.203563	1.026106

各棋士の勝負ムラを表示する。

ret_sigma = rename_indexes(trace, 'sigma', players)
ret_sigma.sort_values('mean', ascending=False).head()

勝負ムラが大きい棋士。

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
藤森哲也	1.126341	0.387405	0.030618	0.353202	1.831846	65.198876	0.999863
森内俊之	1.125223	0.344858	0.026571	0.511552	1.838470	89.927761	1.000393
深浦康市	1.116453	0.337526	0.025087	0.499118	1.743568	112.326690	1.003547
神崎健二	1.113792	0.336298	0.024777	0.483929	1.786775	59.970475	1.000581
阿久津主税	1.113714	0.327681	0.022573	0.472958	1.722987	119.652929	1.023349

勝負ムラが小さい棋士。

ret_sigma.sort_values('mean', ascending=False).tail()

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
菅井竜也	0.906430	0.282381	0.022321	0.375220	1.422869	56.414318	1.013031
先崎学	0.904325	0.279065	0.019915	0.380872	1.418062	103.896750	0.999974
飯島栄治	0.899098	0.250615	0.018413	0.408726	1.358156	72.390061	1.004667
竹内雄悟	0.892474	0.282656	0.019096	0.392110	1.476494	128.417525	1.001682
北島忠雄	0.862707	0.272092	0.018809	0.332110	1.346679	132.050279	1.002114