階層ベイズモデルでぷよらーの強さを推定してみた 

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の結果を用いて分析しました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/e02ca0c2-5ddb-3860-3bd5-5d4bfdfe9d37.png)
リーグ表から、1行1連戦の形式でデータを書き写します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/00314bc0-0497-1084-41ca-07764da76dc6.png)



さらに、1行1試合で敗者・勝者ごとの列になるよう加工します。
↓ 完成系イメージ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/c4fcf676-d984-7886-1268-c6e4b25ae9fb.png)
season1の「ぴぽにあ vs ともくん」（30-28）を例にします。
ぴぽにあ選手が30本取得して勝利したので、W（勝者）列にぴぽにあ選手が、L（敗者）列にはともくん選手が30行分入ります。
逆に、ともくん選手は28本取得していたので、W列にともくん選手が、L列にはぴぽにあ選手が28行分入ります。
この試合に関して、合計58試合行われたので、データにも58行分記載されることになります。


データ加工からstan実行までRを利用しました。

```puyo.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
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/24cf3bc9-3904-24b0-5a0b-8bb3e14dbd31.png)

各ぷよらーの強さ $μ[n]$ は平均0・標準偏差$σ\_μ$の正規分布に従います。
勝負ムラ $σ\_{pf}[n]$ は事前分布であるガンマ分布(パラメータ α=β=10)に従うことを仮定しています。
各試合で発揮されるパフォーマンスは、各ぷよらーの強さ $μ[n]$ と勝負ムラ $σ\_{pf}[n]$ の正規分布から生成されます。

これで各ぷよらーの強さを推定し、試合ごとに調子の波（勝負ムラ）があることを表現しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/d9fbfaa8-d3cd-efa0-4ede-ba21634cd264.png)


```puyo.stan
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の正規分布を近似しつつ、非負の値である勝負ムラを表現しています。

```puyo.R
### 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本先取りの連戦での勝率を可視化しておきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/9dc9f838-a2f1-508d-71bf-72806efad499.png)
(敬称略) 勝率はmomokenが0.88、マッキー0.80と圧倒的な勝率を納めています。次いで、ともくん、delta、ぴぽにあ、SAKIと続きます。

今度はstanで推定した結果を見てみます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/444bea7e-5d48-8b24-e1de-0afdbe208ad6.png)


勝率の順位と若干異なります。勝率では「ともくん = delta > ぴぽにあ」とぴぽにあが若干劣る順位でしたが、推定結果では「ともくん ≧ ぴぽにあ ≧ delta」の順位になっており、3人にほとんど差はないようです。また、直接対決のないSAKI・ヨダソウマ・レインでは、SAKIが最も強い結果になっています。確かに、SAKIはseason2に出場し、ぴぽにあ・ともくんに勝っていたのは印象的でした。


勝負ムラも確認します。最も勝負ムラが小さいのはぴぽにあです。一方、勝負ムラが大きいのはヨダソウマ、レインでした。二人とも勝率では奮わなかったものの、他選手を追い詰めている試合も多々見られたことから、調子の波が大きいと推定されたようです。

最後に、最上級ぷよらーのmomokenとマッキーの「強さ」と「勝負ムラ」の結果を可視化してみます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/466700/cb3b4bb9-dbb0-a25e-43c6-bfa3436ad531.png)
青線がmomokenで、赤線がマッキーです。わずかにmomokenの方が強いと推定されて、マッキーの方が少しだけ「勝負ムラ」が大きいです。

# まとめ
階層ベイズモデルを用いて、最強リーグに出場しているぷよらーの「強さ」・「勝負ムラ」を推定しました。
「強さ」については、勝率通りにmomoken・マッキーが高く、
「勝負ムラ」は、ぴぽにあが最も安定しており、ヨダソウマ・レインが大きかったです。

今回の結果では、記述統計では把握しにくい「勝負ムラ」を推定できたり、直接対決していない選手同士の「強さ」を比較できたり楽しかったです。