2021年4月～9月に、最強のぷよぷよプレイヤー（以下、ぷよらー）を決める「ぷよぷよ最強リーグ」が開催されました。
その勝敗データを用いて各ぷよらーの「強さ」と「勝負ムラ」を階層ベイズモデルで推定します。
ぷよぷよ最強リーグの詳細はこちら → https://shikoukouketsu.com/puyopuyo-saikyou-league/
YouTubeにアーカイブがあるので試合に興味ある方はこちら → https://youtu.be/1V55egyb8p4 (個人的に一番熱い試合 momoken vs マッキー)
今回の分析に利用したコードはこちら → https://github.com/skbono/puyosaikyou2021
データの取得・加工
最強リーグは30本先取りで、結果はリーグ表にまとめられています。
これはseason1の結果ですが、season1 ~ season3の結果を用いて分析しました。
リーグ表から、1行1連戦の形式でデータを書き写します。
さらに、1行1試合で敗者・勝者ごとの列になるよう加工します。
↓ 完成系イメージ
season1の「ぴぽにあ vs ともくん」（30-28）を例にします。
ぴぽにあ選手が30本取得して勝利したので、W（勝者）列にぴぽにあ選手が、L（敗者）列にはともくん選手が30行分入ります。
逆に、ともくん選手は28本取得していたので、W列にともくん選手が、L列にはぴぽにあ選手が28行分入ります。
この試合に関して、合計58試合行われたので、データにも58行分記載されることになります。
データ加工からstan実行までRを利用しました。
### ライブラリの読み込み ### 必要に応じてinstall.packages("$パッケージ名$")でインストール
library(rstan)
library(bayesplot)
library(dplyr)
library(magrittr)
library(readr)
### データの読み込み（1行1連戦 形式） ###
dat <- read_csv("data/SaikyouPuyo2021.csv", locale = readr::locale(encoding = "CP932"))
dat %<>% select(P1, P2, score1, score2)
### 1行1連戦 → 1行1試合に加工 ###
W <- NULL
L <- NULL
for (i in 1:nrow(dat)) {
# Winner列
W_tmp1 <- rep(dat$P1[i], dat$score1[i])
W_tmp2 <- rep(dat$P2[i], dat$score2[i])
W <- c(W, W_tmp1, W_tmp2)
# Loser列
L_tmp1 <- rep(dat$P2[i], dat$score1[i])
L_tmp2 <- rep(dat$P1[i], dat$score2[i])
L <- c(L, L_tmp1, L_tmp2)
}
puyoWL <- data.frame(L, W)
このような1行１試合形式をインプットとし、後ほど紹介するモデルを回します。
最強リーグはseason1 ~ season3まで計3回、seasonごとに選手の入れ替えもあり、計8人が最強リーグに参加しており、合計3898試合となりました。
モデルの作成
「StanとRでベイズ統計モデリング」（通称：アヒル本）のP189にある、プロ棋士の強さを推定するモデルを参考にしました。
https://www.amazon.co.jp/dp/B07M8LWLS1?tag=maftracking142127-22&linkCode=ure&creative=6339
各ぷよらーの強さ $μ[n]$ は平均0・標準偏差$σ_μ$の正規分布に従います。
勝負ムラ $σ_{pf}[n]$ は事前分布であるガンマ分布(パラメータ α=β=10)に従うことを仮定しています。
各試合で発揮されるパフォーマンスは、各ぷよらーの強さ $μ[n]$ と勝負ムラ $σ_{pf}[n]$ の正規分布から生成されます。
これで各ぷよらーの強さを推定し、試合ごとに調子の波（勝負ムラ）があることを表現しています。
data{
int<lower=1> P;
int<lower=1> G;
int<lower=1, upper=P> LW[G, 2];
}
parameters{
ordered[2] performance[G];
vector[P] mu;
real<lower=0> s_mu;
vector<lower=0>[P] s_pf;
}
model{
for (g in 1:G)
for (i in 1:2)
performance[g, i] ~ normal(mu[LW[g, i]], s_pf[LW[g, i]]);
mu ~ normal(0, s_mu);
s_pf ~ gamma(10, 10);
}
stanファイルに記載したモデルは、アヒル本のP189から引用したものを流用しています。
このモデルの重要な部分はparameterブロックのordered[2] performance[G]で、敗者より勝者のperformanceが大きくなるよう大小関係の制約を与えています。
また、勝負ムラには弱情報事前分布としてgamma(10, 10)を設定しています。gamma(10, 10)は、平均1・標準偏差0.32の正規分布を近似しつつ、非負の値である勝負ムラを表現しています。
### stan用にリスト化 ###
P <- c(puyoWL$L, puyoWL$W) %>% unique() %>% length() #プレイヤー数
G <- nrow(puyoWL) #総試合数
data <- list(P=P, G=G, LW=puyoWL) #stanに渡すためにリストにまとめる
### MCMCの実行 ###
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
stanmodel <- stan_model(file='puyo.stan') #Stanファイルの指定
fit <- sampling(stanmodel,
data=data, pars=c('mu','s_mu','s_pf'),
chain = 4,
iter = 10000,
warmup = 5000,
thin = 3,
seed = 1)
### サンプリング結果 ###
# 「強さ」の事後分布の範囲
mcmc_intervals(fit,
pars = c("mu[2]", "mu[6]", "mu[4]", "mu[5]", "mu[1]", "mu[3]", "mu[8]", "mu[7]"),
prob = 0.5,
prob_outer = 0.95)
# 「勝負ムラ」の事後分布の範囲
mcmc_intervals(fit,
pars = c("s_pf[2]", "s_pf[6]", "s_pf[4]", "s_pf[5]", "s_pf[1]", "s_pf[3]", "s_pf[8]", "s_pf[7]"),
prob = 0.5,
prob_outer = 0.95)
R側ではデータを読み込み、必要なデータ「プレイヤー数・総試合数・勝敗データ」をリスト形式にまとめておきます。
あとはsampling関数を実行します。
ちなみに、iteration回数を指定しないと収束してくれなかったので、iteration=10000にしました。
結果
推定結果の前に、各ぷよらーの30本先取りの連戦での勝率を可視化しておきます。
(敬称略) 勝率はmomokenが0.88、マッキー0.80と圧倒的な勝率を納めています。次いで、ともくん、delta、ぴぽにあ、SAKIと続きます。
今度はstanで推定した結果を見てみます。
勝率の順位と若干異なります。勝率では「ともくん = delta > ぴぽにあ」とぴぽにあが若干劣る順位でしたが、推定結果では「ともくん ≧ ぴぽにあ ≧ delta」の順位になっており、3人にほとんど差はないようです。また、直接対決のないSAKI・ヨダソウマ・レインでは、SAKIが最も強い結果になっています。確かに、SAKIはseason2に出場し、ぴぽにあ・ともくんに勝っていたのは印象的でした。
勝負ムラも確認します。最も勝負ムラが小さいのはぴぽにあです。一方、勝負ムラが大きいのはヨダソウマ、レインでした。二人とも勝率では奮わなかったものの、他選手を追い詰めている試合も多々見られたことから、調子の波が大きいと推定されたようです。
最後に、最上級ぷよらーのmomokenとマッキーの「強さ」と「勝負ムラ」の結果を可視化してみます。
青線がmomokenで、赤線がマッキーです。わずかにmomokenの方が強いと推定されて、マッキーの方が少しだけ「勝負ムラ」が大きいです。
まとめ
階層ベイズモデルを用いて、最強リーグに出場しているぷよらーの「強さ」・「勝負ムラ」を推定しました。
「強さ」については、勝率通りにmomoken・マッキーが高く、
「勝負ムラ」は、ぴぽにあが最も安定しており、ヨダソウマ・レインが大きかったです。
今回の結果では、記述統計では把握しにくい「勝負ムラ」を推定できたり、直接対決していない選手同士の「強さ」を比較できたり楽しかったです。