麻雀は実力と運が勝敗を左右するゲームですが、運と実力の比率はどのくらいなのでしょうか？
7：3くらい？ 8：2くらいでしょうか？ あるいは100：1？
RとStanを用いて、この問いに答えてみたいと思います。
使用データ
データには、オンライン対戦麻雀【天鳳】の公開しているログ（https://tenhou.net/sc/raw/）による、プレイヤーの着順履歴を使用します。
サンプルサイズを絞るため、今回は2018年〜2020年の鳳凰卓４人打ち東風戦のデータを利用します。*1
さらにサンプルサイズを効率的に絞るため、 対戦履歴を「3年間で合計500ゲーム以上プレイしているプレイヤー4人による試合」に限ることとします。
すなわち、鳳凰卓常連 の4人による対戦のみに注目します。*2
整形したデータは以下のようになっています。
3年間で東風戦を500戦以上プレイした全226プレイヤーによる、計40657試合の履歴です。
モデルの説明
モデルの概要を説明します。
各プレイヤーの潜在的な強さがあります。
麻雀が完全実力ゲームであれば、この潜在的な強さに従って、この場合、１位から４位が「Aさん→Bさん→Cさん→Dさん」の順となります。
しかし麻雀は完全実力ゲームではないので、運要素が発生します。
運要素は、各プレイヤーに独立に、正規分布に従うものとして発生すると仮定します。
実際に各プレイヤーに正規分布乱数としての運要素が付加されます。
こうして実力要素と運要素の両方から、最終的な着順が決まります。
この例では、「Bさん→Aさん→Dさん→Cさん」の順位になりました。
最後に、いくつか条件を付け加えます。
まず、各プレイヤーの強さも、何らかの正規分布にしたがっていると仮定します。
このとき、「強さ」とは相対的なものなので、この正規分布の平均はゼロに固定できます。
また、先の例において、運要素、実力要素の数字のオーダーは不定です。
（運要素10に対して実力要素1と、運要素1に対して実力要素0.1、は同じことです）
ですので、仮定した運要素の正規分布の標準偏差を1に固定します。
そして「各プレイヤーの強さの正規分布」の標準偏差をXとおき、これを推定します。
すなわち以下のようなイメージです。
麻雀には実力でどうこうできる標準偏差Xの正規分布と、実力ではどうにもならない運としての標準偏差1の正規分布、二つの要素から勝敗が決定する、ということです。
これにより「麻雀は運 : 実力 = 1 : X である」と推定できます。
なお、実際のStanコードは、松浦健太郎氏による書籍の、第10章のものを使用しました。
結果
モデルの収束には12時間ほどかかりました。
実力の分布である正規分布の標準偏差Xは、中央値で0.0183と推定されました（下図）。
よって、麻雀の運と実力の比は、1:0.0183 ... すなわち、およそ50:1であると分かりました。
解釈
結果の解釈に際して、いくつか注意すべきと思われる点を列挙します。
「勝負ムラ」を「運要素」と捉えている（完全実力ゲームなど存在しない？）
将棋を例に挙げます。
将棋は運要素が極端に少ないことは、多くの人が同意する点であると思います。
ほとんど将棋のルールを知っているだけの私が、羽生善治さんと対戦しても、一度も勝てないでしょう。
しかし、藤井聡太さんと羽生善治さんならどうでしょうか？
今回使用したモデルに合わせて考えると、仮に藤井聡太さんの強さを49、羽生善治さんの強さを51（逆でも構いません）とすると、「完全実力ゲーム」なる将棋では、100回やって100回羽生善治さんが勝たねければなりません。
しかし、これはありそうもないことです。
この説明として、先に紹介した松浦健太郎さんの書籍では、「勝負ムラ」という単語が使われています。
つまり、その日（あるいはその対局）におけるパフォーマンスが、二人の強さ（49と51）を中心にして確率的に振る舞うと考えることで解決します。
これを運要素と捉えるかどうかは、何とも言い難い問題であるような気がしますが、
本モデルではこれを「運要素」として一纏めにしていることには注意しておきたいと思われます。
そもそもレベルの高いフィールドでの話をしている点
これは、先の将棋の例と繋がります。
本解析では、天鳳鳳凰卓という、一定の雀力がないと対戦権を得られないフィールドでの履歴を参照しました。
しかし実際には、特上卓・上級卓・一般卓という「鳳凰卓以下」のフィールドも存在します。
なにが言いたいかというと、レベルの高いフィールド（正確には、同程度のレベルの集合）では、「勝ったり負けたり」が発生しがちである...すなわち、「運要素」が大きくなるということです。
将棋においても、私が最上位層のプロ棋士に挑んでも一向に歯が立たないでしょうが、
最上位層同士の争いであれば、「勝ったり負けたり」が発生し、
「運要素（これを運要素と捉えるかは問題ですが）」が大きくなることが分かるでしょう。
この点で、本当はさまざまなレベルのプレイヤーが一纏めに対戦する場の履歴を解析すれば良かったのですが、
ネット麻雀ではレベルによる卓分けが明確になされている分、対処は難しかったです。
いろんなレベルの人たちが集まるフリー麻雀の対戦履歴が得られれば一番いいですね。
そのため、今回の解析結果の解釈としては、
「麻雀」というよりも「鳳凰卓東風戦」の運と実力比が50:1であると捉えておくのが無難です。
「個人ごとのパフォーマンスのブレの違い」を考慮していない
今回、運要素を全て「平均ゼロ・標準偏差1」に固定しましたが、実際には、プレイヤーごとにいくらかブレ幅に違いがあることが予想されます。
たとえば、将棋と同様に「その日その対局におけるコンディション」による違いや、
あるいは採用しがちな戦術によるブレ幅もあるでしょう（たとえば「トップラス麻雀」の人はブレ幅が大きいと言えます）
ここではそのようなプレイヤーごとの違いをモデル化していない（正確には、全て同じであるとしてモデル化している）ことにも、一応注意しておきたいところだと思います。
終わりに
以上です。
今回得られた結果から安定段位*3に関するシミュレーションもできたので、次はそれを書いてみます。
質問・コメント等あれば是非お願いします。
8/31 追記
思ったよりも多くの方に見られて緊張しているので、追記です。
このモデルのポイントは
・ 運要素に正規分布を仮定する
・ プレイヤーの実力に正規分布を仮定する
・ プレイヤーはそれぞれの実力に応じた「下駄」を履いた状態で
・ 正規分布乱数としての運要素が独立に付加され
・ 最終的に着順が決まる
というものです
この４つ目の仮定に関しては、実は現実的ではないかもしれません。
なぜかというと、麻雀は「誰かが点棒を失えば、誰かが点棒を得る」からです。
すなわち、各プレイヤーに作用している運要素は、実際には独立ではないことが考えられます。
この点を改良したモデルを作れないかなぁと考えているので、できたら追記します。
モデルの適合の良さについては言及していなかったので、さらっと書いておきます。
推定されたパラメータ0.0183を用いて正規分布から発生させた色々な強さを持つ1000プレイヤーについて、総当たり、一人およそ2000対戦させ、
1000人分の安定段位を集計した結果が以下です。
男冥利さん（http://otokomyouri.com/toppage.aspx）というサイトからとってきた、実際の鳳東プレイヤーの安定段位の分布が以下です。
実力の標準偏差を0.05と0.1でシミュレートすると以下のようになります。
横軸の値に注目してください。
実力の幅を過大評価しているせいで、現実には誕生し得ないほどの安定段位を持った仮想プレイヤーが爆誕しています。
なお、安定段位の低い方の端は適合が悪そうですが、これは現実にはプレイヤーが鳳凰卓での対戦権を失ってしまい、2000ゲームという対戦数を重ねられないことが原因と考えられます。
以上から、簡易的な評価ですが、このモデルはまだ改良の余地がありそうですが、適合のほどは悪くないと思います。
Rコード
######################### 2019-2020年分の鳳東データの前処理 ##############################
# データの読み込み
files2018 <- list.files("data/2009.2.20_2021.8.3鳳凰卓scc", pattern = "scc2018", full.names = TRUE)
files2019 <- list.files("data/2009.2.20_2021.8.3鳳凰卓scc", pattern = "scc2019", full.names = TRUE)
files2020 <- list.files("data/2009.2.20_2021.8.3鳳凰卓scc", pattern = "scc2020", full.names = TRUE)
d_orig2018 <- do.call(rbind, lapply(files2018, read_delim, delim=" | ", quote="",
col_names = c("X1","X2","X3","X4","X5","X6","X7","X8","X9","T1","T2","T3","T4")))
d_orig2019 <- do.call(rbind, lapply(files2019, read_delim, delim=" | ", quote="",
col_names = c("X1","X2","X3","X4","X5","X6","X7","X8","X9","T1","T2","T3","T4")))
d_orig2020 <- do.call(rbind, lapply(files2020, read_delim, delim=" | ", quote="",
col_names = c("X1","X2","X3","X4","X5","X6","X7","X8","X9","T1","T2","T3","T4")))
d_orig <- rbind(d_orig2018, d_orig2019, d_orig2020)
# 四鳳南喰赤 を選択
d <- filter(d_orig, X5 == "四鳳東喰赤－" | X5 == "四鳳東喰赤")[,10:13]
# 不要な文字列を削除
d$T1 <- str_remove(d$T1, pattern = "\\(\\+.*?\\)")
d$T2 <- str_remove(d$T2, pattern = "\\(\\+.*?\\)")
d$T2 <- str_remove(d$T2, pattern = "\\(\\-.*?\\)")
d$T2 <- str_remove(d$T2, pattern = "\\(0.0\\)")
d$T3 <- str_remove(d$T3, pattern = "\\(\\+.*?\\)")
d$T3 <- str_remove(d$T3, pattern = "\\(\\-.*?\\)")
d$T3 <- str_remove(d$T3, pattern = "\\(0.0\\)")
d$T4 <- str_remove(d$T4, pattern = "\\(\\-.*?\\)")
d$T1 <- str_remove(d$T1, "<br>")
d$T2 <- str_remove(d$T2, "<br>")
d$T3 <- str_remove(d$T3, "<br>")
d$T4 <- str_remove(d$T4, "<br>")
# 4,3,2,1着順に並び替え
d <- d[,c(4,3,2,1)]
#### データの把握 ####
# 全対戦数
nmatch <- nrow(d)
# 縦長データに変換
d_long <- pivot_longer(d, cols=everything(), names_to="order", names_prefix="T", values_to="player")
d_long$match <- rep(1:nrow(d), each=4)
# 全プレイヤーの総対戦数
M <- arrange(count(d_long, player), n)
colnames(d) <- c("forth", "third", "second", "first")
##### 対戦数少プレイヤーを省く #####
# 3年間で対戦数500以上のプレイヤー
sum(M$n>=500) # 226player
# 対戦数500以上のプレイヤー名を抽出
P500 <- M$player[which(M$n>=500)]
res.p500 <- is.element(d_long$player, P500) # 対戦数500以上のプレイヤーが含まれるか？
# 対戦数が500以上のプレイヤーが4人以上含まれる試合ナンバーを抽出
m500 <- data.frame(match = d_long$match[which(res.p500 == TRUE)])
tmp <- m500 %>%
group_by(match) %>%
filter(n()>3)
M500 <- unique(tmp$match)
# 対戦数500以上のプレイヤーが4人以上含まれる試合を選抜
d500 <- d[M500, ]
#' *対戦数500以上のプレイヤーが4人以上含まれる試合で解析*
# 個々人のID名に対応するナンバー
d500_long <- pivot_longer(d500, cols=everything(), names_to="order", names_prefix="T", values_to="player")
conv <- 1:length(unique(d500_long$player))
names(conv) <- unique(d500_long$player)
LtoW <- array(NA, dim=c(nrow(d500), ncol(d500)))
for(i in 1:nrow(d500)){
for(j in 1:ncol(d500)){
LtoW[i,j] <- conv[as.character(d500[i,j])]
}
}
nmatch500 <- nrow(d500)
data <- list(nind=max(conv), nmatch=nmatch500, LtoW=LtoW)
#' *RUN*
params <- c("s_mu", "mu")
tonpu500 <- stan(file="analysis1.stan", data=data, seed=1234, pars=params,
thin=3, iter=6000, warmup=3000, chain=3)
Stanコード
data{
int nind;
int nmatch;
int<lower=1, upper=nind> LtoW[nmatch, 4];
}
parameters{
real<lower=0> s_mu;
vector[nind] mu;
ordered[4] performance[nmatch];
}
model{
// 事前分布
s_mu ~ normal(0, 0.3);
for(m in 1:nmatch){
for(k in 1:4){
performance[m,k] ~ normal(mu[LtoW[m,k]], 1);
}
}
mu ~ normal(0, s_mu);
}
