ベイズモデルで大相撲力士の強さを判定してみた

はじめに

以前、RのBradleyTerry2というパッケージで各力士の強さを判定しました。
既存のパッケージを使いとてもお手軽に力士の強さを数値化できましたが、大相撲好きとしてはやはり、大相撲の特徴を捉えた強さ判定モデルを自分で考えて数値化したいところです。

ということで、PyStanを使ってモデル化してみました。

各力士の強さ

モデル化

モデル化にあたって考慮すべき大相撲の特徴として、かなり実力差のある力士同士の取組でも、ちょこちょこ番狂わせが起こることです。
平成の大横綱である白鵬であっても、下位力士に金星を配給することはままあります。

アヒル本を参考に、各力士の強さをmu、勝負ムラをsigmaとし、1回の勝負で発揮する力はmu、sigmaをパラメータとするコーシー分布（正規分布より裾がかなり長い分布）から生成され、勝敗はその力の大小で決まると考えます。

model = '''
data {
    int N;
    int M;
    int<lower=1, upper=N> Id[M, 2];
}

parameters {
    vector[N−1] mu0;
    real<lower=0> s_mu;
    vector<lower=0>[N] sigma;
    ordered[2] performance[M];
}

transformed parameters {
    vector[N] mu;
    mu[1:N-1] = mu0;
    mu[N] = -sum(mu0);
}

model {
    mu ~ normal(0, s_mu);
    sigma ~ gamma(10, 10);
    for (i in 1:M)
        for (j in 1:2)
            performance[i, j] ~ cauchy(mu[Id[i, j]], sigma[Id[i, j]]);
}
'''

下の様に、取組毎の敗者と勝者の力士ID、その取組が行われた年をまとめたデータ（全期間における全取組数が10未満の力士の取組は除きました）を、

	loser	winner	year
0	1	40	2001
1	2	41	2001
2	3	21	2001
3	4	31	2001
4	5	20	2001

下の様に与えました。

N = loserwinner['loser'].max()
M = loserwinner.shape[0]
Id = list(loserwinner[['loser', 'winner']].values)

data = dict(
    N=N,
    M=M,
    Id=Id
)

stanmodel = StanModel(model_code=model)
fit = stanmodel.vb(data=data)

2001年以降のほぼ全ての取組を対象とすると30000行を超すデータ量であり、PyStan標準のNUTSではchains=3, iter=400, thin=1で1時間以上かかってしまいました。
ADVIでやってみたところ、5分もかからずにサンプリングが終了しました。

可視化

ベイズモデリングで男子プロテニスの強さを分析してみたを参考に可視化しました。

plt.figure(figsize=(40,10))
cmap = plt.cm.get_cmap('tab10')

for i in range(strength.shape[1]):
    g = plt.violinplot(strength.iloc[:, i], positions=[i], showmeans=False, showextrema=True, showmedians=True)
    c = cmap(i%10)
    for pc in g['bodies']:
        pc.set_facecolor(c)
    g['cbars'].set_edgecolor(c)
    g['cmaxes'].set_edgecolor(c)
    g['cmedians'].set_edgecolor(c)
    g['cmins'].set_edgecolor(c)

plt.xticks(list(range(len(rikishi_id.keys()))), rikishi_id.keys())
plt.xticks(rotation=90)

plt.xlabel('rikishi')
plt.ylabel('strength')
plt.savefig('output/rikishi_bayes_2001to2018_advi', dpi=200)
plt.show()

これが全力士の強さです。
見にくいので上位20人だけ見てみます。

白鵬、朝青龍が圧倒的ですね。
武蔵丸や貴乃花は引退前の2,3年分のデータしか対象としていないので、プロットが縦に広がっています。
現役横綱である稀勢の里と鶴竜は、最高位が大関であった栃東、把瑠都、魁皇より弱いと判定されています。

現役力士に絞って見てみます。

中央値が大きく、プロットが広がっている、阿武咲、御嶽海、貴景勝、北勝富士、阿炎、朝乃山は、若手注目力士であると言えそうです。

強さの推移（直近5年）

モデル化

年毎の強さの推移を知りたいので、1階差分のトレンド項を導入してみます。

model = '''
data {
    int N;
    int M;
    int L;
    int<lower=1, upper=N> Id[M, 3];
}

parameters {
    ordered[2] performance[M];
    matrix[L, N-1] mu0;
    real<lower=0> s_mu;
    real<lower=0> s_mu_diff;
    vector<lower=0>[N] sigma;
}

transformed parameters {
    matrix[L, N] mu;
    mu[, 1:N-1] = mu0;
    for (i in 1:L)
        mu[i, N] = -sum(mu0[i, ]);
}

model {
    mu[1, ] ~ normal(0, s_mu);
    for (k in 2:L)
        mu[k, ] ~ normal(mu[k-1, ], s_mu_diff);
    sigma ~ gamma(10, 10);
    for (i in 1:M)
        for (j in 1:2)
            performance[i, j] ~ cauchy(mu[Id[i, 3], Id[i, j]], sigma[Id[i, j]]);
}
'''

N = loserwinner['loser'].max()
M = loserwinner.shape[0]
L = loserwinner['year'].max()
Id = list(loserwinner[['loser', 'winner', 'year']].values)

data = dict(
    N=N,
    M=M,
    L=L,
    Id=Id
)

可視化

2018年時点、上位10人の中央値の推移です。

なかなかいい感じに判定できている気がします。

まとめ

白鵬、朝青龍は強すぎる。
稀勢の里、鶴竜はもう少し勝たないと、真の横綱として広く国民に認められなさそう。
阿武咲、御嶽海、貴景勝、北勝富士、阿炎、朝乃山の今後に期待。

コードは全てgithubにあります。プログラミング、相撲、ともに初心者なのでご教授いただけたらうれしいです。