麻雀の実力分布の推定

結果から実力が知りたいんですけど、どう考えれば良いですか？
というのがこのノートのテーマです。

基本的には麻雀を題材にしていますが、別にポーカーであってもスポーツや将棋であってもある程度同じ考え方で問題ありません。(一応)

はじめに

遠い昔に、同じようなテーマでnoteを書いてみたのですが、自分でもよく分からないまま突っ走ってしまっていました、、
ありがたい意見や自分の理解が追いついてきたおかげでどうすれば良いのかが見えてきて、改訂版を出したいな、とずっと思ってはいたのですがなかなか筆が進まず、、
良い感じの安定段位の新しい指標を導出することもできたので前半部分がよく分からない方もその部分だけでもよかったら見ていただけたらな、、と思います。
最後は、色んなルールで試行回数を変えるとどれくらい実力が出るのか？という点を掘り下げてみていきます。

前回は今見ると無茶苦茶やってるな、という感じなのですが今回はそこそこな出来に仕上がったかなと思います。
コードと使用したデータもあげておくので動かしたい人はぜひ動かしてみてください。天鳳以外のデータも取り上げてみようと思ったのですがデータを集めるのがちょっと面倒になったので、なんかデータを持っている人はぜひやってみて結果を教えていただきたいです。それかデータをくれると僕がやります。

面白かったり、参考になったなと思った方は感想をいただけると喜びます。
内容に関する間違いの指摘やご意見等もいただけたら大変幸いですので是非ともよろしくお願いします。

では少し長くなりますが、お付き合いよろしくお願いします。

なぜ、なにをどうしたいのか

麻雀をしている人だと、
「100試合だと平均着順が±0.25の範囲でブレる」
なんてことを聞いたこともあると思います。
しかしながら、私たちはブレた後の結果しか見ることができません。

では、平均着順2.25の結果があったときに、その人の本当の平均着順は2.0か2.5のどっちの確率が高いのでしょうか？
同じではなく、間違いなく平均着順2.5である確率の方が高いでしょう。

この違いは事前分布から来ています。
わかりやすく、この世には
真の平均着順2.0の人が1人
真の平均着順2.5の人が30人
それ以外はいないとしましょう。
さらに2.0の人は100%平均着順2.0が出せて、2.5の人は30%の確率で平均着順2.0が出せるとします。
この時、平均着順2.0の結果を見て、その人が本当の平均着順2.0の人である確率は

$$
\def\arraystretch{2}
\begin{array}{}
P &=& \frac{1}{平均着順が2.0になる人数} \\
   &=& \frac{1}{1+30\times0.3} \\
   &=& 10\%
\end{array}
$$

となるので2.0の人が100%の確率で平均着順2.0の結果を叩き出せても、結果から実力を推定すると90%の確率で平均順位2.5の人だろうとなるわけなんですね。
もちろんこれは極端な例なのですが、初めにおいた人数の分布が非常に重要であることがわかると思います。
結果から実力を語るには実力の分布が必要になってしまうのです。

しかし難しいのが、私たちは結果を見ることはできても、真の実力の分布を見ることが出来ません。なので、それを結果から推測しましょう、というのが今回したいことになります。
また事前分布を考慮した統計値は考慮しない統計値よりも良い指標となってくれることを示します。

事前分布のパラメータ推定

推定方法

この章は少し数式多めになります、、、
全然読み飛ばしてください、、、

まず事前分布をディリクレ分布でおきます。
パラメータ$${(\alpha_1,\alpha_2,…,\alpha_k)}$$
であるとき、着順確率$${(x_1,x_2,…,x_k)}$$に対して

$$
P(\boldsymbol{x};\boldsymbol{\alpha})=
\frac{\Gamma(\sum_{i=1}^k \alpha_i)}{\Gamma(\alpha_1)\Gamma(\alpha_2)\cdots\Gamma(\alpha_k)}
x_0^{\alpha_0-1}x_1^{\alpha_1-1}\cdots x_k^{\alpha_k-1}
$$

となります。
着順分布は対称性があるので$${\alpha=\alpha_1=…=\alpha_k}$$です。
$${k}$$は次元数で、例えば4人麻雀なら4です。
$${\Gamma}$$関数は階乗を一般化したような関数で
$${\Gamma(n+1)=n!}$$です。
$${P(A;B)}$$はBである時にAである確率です。

式だけだと形が想像しにくいですが、グラフはこんな感じです。

3次元の場合だと$${\boldsymbol{x}=(0.33,0.33,0.33)}$$を中心とする山だと思ってもらって大丈夫です。$${\alpha}$$は山の広がりを表していて、$${\alpha}$$が小さいほど広くなだらかな山になっています。

$${\alpha}$$は一様分布だと仮定して$${\alpha}$$を推定します。
すなわち
N人のそれぞれの着順回数のデータ$${Y}$$があったときに、最も確率が高くなる$${\alpha}$$を探します。

$$
\def\arraystretch{2.5}
\begin{array}{}
arg \max_\alpha P(\alpha;Y) &=&arg \max_\alpha  \frac{P(Y;\alpha)P(\alpha)}
{P(Y)}
\\
&=& arg\max_\alpha P(Y;\alpha)
\end{array}
$$

ベイズの定理ってやつです。$${P(\alpha)}$$は一様分布であり$${P(Y)}$$は$${\alpha}$$を含まないので、結局$${P(Y;\alpha)}$$を最大化する$${\alpha}$$を求めることになります。

今回各人のデータはそれぞれ独立であると考えます。
そうすると単純に各人の確率の掛け算となって、

$$
\def\arraystretch{3}
\begin{array}{}
P(Y;\alpha) &=& \prod_i^n P(\boldsymbol{y}^j;\alpha)\\

   &=& \prod_i^n \int\int\cdots\int_{\boldsymbol{x}}P(\boldsymbol{y}^j;\boldsymbol{x})P(\boldsymbol{x};\alpha)\\

   &=& \prod_i^n \int\int\cdots\int_{\boldsymbol{x}}
\frac{\Gamma(\sum_j^k y_j^i+1)}{\Gamma(y^i_0+1)\Gamma(y^i_1+1)\cdots \Gamma(y^i_k+1)}
x_0^{y^i_0}x_1^{y^i_1}\cdots x_k^{y^i_k}\cdot
\frac{\Gamma(k\alpha)}{\Gamma(\alpha)^k}x_0^{\alpha-1}x_1^{\alpha-1}\cdots x_k^{\alpha-1}\\

&=& \prod_i^n 
\frac{\Gamma(k\alpha)}{\Gamma(\alpha)^k}
\frac{\Gamma(\sum_j^k y_j^i+1)}{\prod_j^k \Gamma(y_j^i +1)}
\int\int\cdots\int_{\boldsymbol{x}}
x_0^{y_0^i+\alpha-1}x_1^{y_1^i+\alpha-1}\cdots x_k^{y_k^i+\alpha-1}\\

&=&
\frac{\Gamma(k\alpha)^n}{\Gamma(\alpha)^{kn}}
 \prod_i^n 
\frac{\Gamma(\sum_j^k y_j^i+1)}{\prod_j^k \Gamma(y_j^i +1)}
\frac{\prod_j^k \Gamma(y_j^i+\alpha)}{\Gamma(\sum_j^k y_j^i+k\alpha)}

\end{array}
$$

何が何やらという感じですが、何とか$${P(Y;\alpha)}$$を計算出来そうな形になりました。$${\Gamma}$$関数の性質を使えばもう少し簡略化できるような気もしなくもないけど、コンピュータで計算する前提なのでこれくらいなら問題ないでしょう。
後は、このままだと精度の部分で取り回しが悪いので対数を取って和の形にして最小化問題に変えてから数値計算でゴリゴリ$${\alpha}$$を求めましょう。

ランダムなデータを作ってみて、うまく動いているかの確認ができるコードをあげておきます。グラフ描画も含めると完了まで1~2分かかります。
実行ボタンポチッとするだけなので携帯でも動かせます。たぶん

Google Colaboratory

こんな結果が見れます。
ずれることもありますが、そこそこちゃんと動いてそうですね。

実際のデータでやってみる

ということでパラメータを推定できるようになったので、次にやるのは当然実際のデータで計算してみることですよね。

今回使ったデータです。天鳳の各ルールでの各プレイヤーの着順回数のデータで、json形式のファイルです。
[1着,2着,3着,4着](三麻は[1着,2着,3着])の回数で、確か2010年~2019年の期間です。
欲しい人はダウンロードしてみてください。

ということでやってみました。

$$
\def\arraystretch{2.5}
\begin{array}{ccc|ccc}
ルール& \alpha&データ数&ルール& \alpha &データ数\\ \hline
三特東&72.4&11182&三鳳東&150.0&631\\
三特南&64.6&45898&三鳳南&191.6&3115\\ \hline
四特東&37.5&56527&四鳳東&262.1&4611\\
四特南&51.0&79383&四鳳南&194.5&7874\\
\end{array}
$$

さて、もう一度$${\alpha}$$の意味をおさらいします。
$${\alpha}$$は山の鋭さを示しており、$${\alpha}$$が大きいほど分布の山は鋭く、すなわち実力の分布が集団の平均に集まっていることを表します。この実力分布はゲームの本質的な差の出やすさを含んでいますので、各ルールの集団のスキルの部分は簡単に計ることはできないです。集団のスキルが同じであると仮定するならば、$${\alpha}$$が大きいほど差の出にくいゲームであるとは言えます。
次元数が違うと、$${\alpha}$$の比較は意味がありません。異なる次元数での比較は少し後回しにして現状は同じ人数でのルールで$${\alpha}$$を比較してみます。

感覚的には、特<鳳,南<東の関係になるだろうかと思いましたが、南<東にはなっていない関係もありますね。
四特東<四特南はどうしてだろうか、、、
そういうこともあるんだろうか、、わかりません。
三鳳東はデータ数も少なく、各データでの対局数も少ないので信頼度が低いのかもしれません。

特<鳳の関係は全て満たしていますね。何だか安心です。

事前分布を考慮した統計値

期待値

さて、真に実力を表すものは何かというと着順の分布ではなく、期待値になります。ルールによって色々と変わると思いますが、期待値は着順の分布がわかっていれば簡単に計算することが出来ます。
しかしながら、何度も言うように私たちは結果しか見ることが出来ず、真の着順分布を見ることはできません。
平均着順の話と同じで、期待値も事前分布を考慮した期待値を考える必要があります。

獲得ポイントを$${d_1,d_2,\cdots,d_k }$$、結果の着順回数を$${y_1,y_2,\cdots,y_k }$$とすると、求めたいものは

$$
\begin{array}{}
E(\boldsymbol{x}^T\boldsymbol{d};\boldsymbol{y})&=&
\int_{\boldsymbol{x}} \boldsymbol{x}^T\boldsymbol{d}P(\boldsymbol{x};\boldsymbol{y)} \\
\end{array}
$$

事後分布$${P(\boldsymbol{x};\boldsymbol{y})}$$を切り離して求めます。邪魔くさいので、定数項をCとおくと

$$
\def\arraystretch{2.5}
\begin{array}{}
P(\boldsymbol{x};\boldsymbol{y})&=&\frac{P(\boldsymbol{y};\boldsymbol{x})P(\boldsymbol{x})}{P(\boldsymbol{y})}\\
&=& C\prod_i^k x_i^{y_i+\alpha-1}
\end{array}
$$

となるので、ディリクレ分布となります。
この時、$${E(\boldsymbol{x}^T\boldsymbol{d};\boldsymbol{y})}$$もいくつかのディリクレ分布の積分の和となることが分かるので、

$$
\def\arraystretch{2.5}
\begin{array}{}
\frac{\Gamma(\alpha_1)\Gamma(\alpha_2)\cdots\Gamma(\alpha_k)}{\Gamma(\sum_{i=1}^k \alpha_i)}&=&\int_{\boldsymbol{x}}\prod_i^k x_i^{y_i+\alpha-1}\\
\Gamma(n+1)&=&n\Gamma(n)
\end{array}
$$

これらを用いて整理していくと、

$$
E(\boldsymbol{x}^T\boldsymbol{d};\boldsymbol{y})=\frac{\sum_i^k (y_i+\alpha)d_i}{k\alpha+\sum_i^k y_i}
$$

だいぶスッキリした形にはなったけど、なんというか、はいそうですか、という感じですね。

安定段位

ここからは四鳳南のルールで話をします。
天鳳には安定段位というものがありますね。
安定しないやらバイアスがあるやらなどよく聞く話ですが、求め方としては
$${\boldsymbol{d}=(6,3,0,-(X+2))}$$として期待ポイントが0となる$${X}$$について解きます。そうすると

$$
X=\frac{6y_1+3y_3-2y_4}{y_4}
$$

となります。
この式を解く時に4着率で割り算をするという部分が、統計的には諸悪の根源ではあるのですが、目安段位としてわかりやすい指標が得られるという点は評価できる点です。

では、これを先ほど求めた期待値の式で$${X}$$を計算してみるとどうなるのでしょうか。
そうすると

$$
\widehat{X}=\frac{6y_1+3y_3-2y_4+7\alpha}{y_4+\alpha}
$$

となります。若干変わりましたね。
$${y}$$が$${\alpha}$$に対して十分大きい時、$${X}$$と$${\widehat{X}}$$は一致します。$${X}$$は回数を増やすと正しい値に近づいていくことが知られているのでこれは嬉しい点です。

$${y}$$が小さい時には値が異なります。以前の式では、$${y}$$同士の比率が同じであれば$${X}$$も同じでしたが、$${\widehat{X}}$$は異なります。

例を見てみます。
次の着順はどれも$${X=8.5}$$となる着順です。
先ほど求めた$${\alpha=194.5}$$を用いて$${\widehat{X}}$$を計算すると

$$
\def\arraystretch{1.5}
\begin{array}{cclccc}
\boldsymbol{y}&=&(2,3,-,2)&\widehat{X}&=&7.02\\
\boldsymbol{y}&=&(20,30,-,20)&\widehat{X}&=&7.14\\
\boldsymbol{y}&=&(200,300,-,200)&\widehat{X}&=&7.76\\
\boldsymbol{y}&=&(2000,3000,-,2000)&\widehat{X}&=&8.37\\
\end{array}
$$

次の着順はどれも$${X=5.5}$$となる着順です。

$$
\def\arraystretch{1.5}
\begin{array}{cclccc}
\boldsymbol{y}&=&(3,4,-,4)&\widehat{X}&=&6.97\\
\boldsymbol{y}&=&(30,40,-,40)&\widehat{X}&=&6.74\\
\boldsymbol{y}&=&(300,400,-,400)&\widehat{X}&=&6.00\\
\boldsymbol{y}&=&(3000,4000,-,4000)&\widehat{X}&=&5.57\\
\end{array}
$$

同じ比率である人同士で、勝っているならば回数を打っている人の方が強いだろう、負けているならば回数を打っている人の方が弱いだろう。
という感覚的にも正しいことを$${\widehat{X}}$$は定量的に評価してくれています。
もちろん対局数と順位比率が異なる場合にも正しく比較することが出来ますし、対局数を考慮することが出来るので対局数が少ないが安定段位が高く出た場合は信頼度がどのくらいなのか、といった無駄な心配もなくなります。4着回数が0であっても問題ありません。

事前分布を考慮した$${\widehat{X}}$$はより良い安定段位の指標と言えそうです。

ちなみに以下は四鳳南における$${\widehat{X}}$$のランキングtop10です。
対局数の制限なく強そうな人を見ることが出来ます。
新しいデータにアップデートして、他のルールでもまたやってみようかなと思います。

$$
\def\arraystretch{1.5}
\begin{array}{cll}
一蓮峻潤&\boldsymbol{y}=(1277, 1168, 1143, 969) &\widehat{X}=9.10\\
トトリ先生19歳&\boldsymbol{y}=(648, 576, 531, 475)&\widehat{X}=9.00\\
黒猫@ぺろぺろ☆&\boldsymbol{y}=(554, 532, 502, 416)&\widehat{X}=8.93\\
zeRo&\boldsymbol{y}=(2141, 1827, 1873, 1651)&\widehat{X}=8.88\\
夕凪&\boldsymbol{y}=(868, 943, 949, 706)&\widehat{X}=8.87\\
水希♪&\boldsymbol{y}=(436, 390, 373, 317)&\widehat{X}=8.82\\
ASAPIN&\boldsymbol{y}=(438, 411, 411, 324)&\widehat{X}=8.82\\
（≧▽≦）&\boldsymbol{y}=(679, 720, 652, 552)&\widehat{X}=8.70\\
のたわ&\boldsymbol{y}=(230, 244, 216, 167)&\widehat{X}=8.68\\
太くないお&\boldsymbol{y}=(1565, 1438, 1383, 1254)&\widehat{X}=8.67\\
\end{array}
$$

ゲーム毎の実力の反映度

データがあれば、事前分布のパラメータ$${\alpha}$$を推定することが出来るようになりました。しかしながら、別のルールやゲームとの関係性はまだわかっていない状態でした。
となれば、次はその点を考えたいです。
例えば、三麻何回と四麻何回なら実力の出やすさが同じくらい、あるいは競技を超えて、将棋何回とサッカー何回で同じくらい、と数字で言えたら面白いわけです。
となるとそもそも、実力の反映度とはどう定義するのがいいのでしょうか。
実力$${\boldsymbol{x}}$$、結果$${\boldsymbol{y}}$$、得点関数$${f}$$があった時に、
実力の反映度は
「実力で勝っている時に、結果でも勝てる確率」
$${P(f(\boldsymbol{y})>0;f(\boldsymbol{x})>0)}$$
と定義します。これはあらゆる競技で妥当な定義だと思います。
ちなみに、プレイヤーの事前情報がなく、$${f(\boldsymbol{x}),f(\boldsymbol{y})}$$の分布が原点に関して対称な場合は

$$
\def\arraystretch{2.5}
\begin{array}{}
P(f(\boldsymbol{y})>0;f(\boldsymbol{x})>0)
&=&\frac{P(f(\boldsymbol{x})>0;f(\boldsymbol{y})>0)P(f(\boldsymbol{y})>0)}{P(f(\boldsymbol{x})>0)}\\
&=&\frac{P(f(\boldsymbol{x})>0;f(\boldsymbol{y})>0)\cdot 0.5}{0.5}\\
&=&P(f(\boldsymbol{x})>0;f(\boldsymbol{y})>0)
\end{array}
$$

となり「結果で勝っている時に真の実力でも勝っている確率」と同じになります。

これらは試行回数$${N}$$によって変化するのが面白いどころです。
具体的に$${P(f(\boldsymbol{y})>0;f(\boldsymbol{x})>0)}$$を求めていきます。

$$
\def\arraystretch{2.5}
\begin{array}{}
P(f(\boldsymbol{y})>0;f(\boldsymbol{x})>0)&=&\frac{P(f(\boldsymbol{y})>0,f(\boldsymbol{x})>0)}{P(f(\boldsymbol{x})>0)=0.5}\\
&=&2P(f(\boldsymbol{y})>0,f(\boldsymbol{x})>0)\\
&=&2\int_{\boldsymbol{f(x)>0}}P(f(\boldsymbol{y})>0;\boldsymbol{x})p(\boldsymbol{x})
\end{array}
$$

$${f(\boldsymbol{x})=\boldsymbol{x}^T\boldsymbol{d}}$$であるとすると、
$${\boldsymbol{y}}$$は独立な$${N}$$回の試行によって決定されるので、
$${f(\boldsymbol{y});\boldsymbol{x}}$$は平均$${N\boldsymbol{x}^T\boldsymbol{d}}$$、分散$${N\boldsymbol{x}^T\boldsymbol{d^2}}$$の正規分布に従うと近似することが出来る。
よって、$${P(f(\boldsymbol{y})>0;\boldsymbol{x})}$$を正規分布の累積分布関数を用いて計算することが出来るので、積分区間に注意して数値積分でゴリおします。

この計算のコードを置いておきます。
良かったら好きなパラメータで試してみてください。

Google Colaboratory

ルールと試行回数による実力の関係

各ルールで$${N}$$を変えながら行うとどのようになるのかを見ていきます。
三鳳南,四鳳東,四鳳南の3つのルールで調べてみました。それぞれhn3,ht4,hn4です。

_revicedとなっているものは、三鳳南,四鳳東の時間は四鳳南の約1/2倍であることを考慮して、横軸方向を2倍に圧縮したものです。すなわち時間が四鳳南と揃っていることになります。

対局数で見ると、三鳳南>四鳳南>四鳳東の順で$${P}$$が大きいのが分かります。時間を揃えると、三鳳南>四鳳東>四鳳南の順で$${P}$$が大きいです。半荘1回よりも東風2回の方が$${P}$$は大きいんですね。

どのルールも5000半荘では、いわゆる95%ラインを超えることはできませんでした、、、

$${N}$$<300の部分を拡大してみるとこうなります。

大体100半荘するとで62%前後みたいです。
逆に38%の確率で実力で勝ってるのに負けるってことですね。
思ったより高いのか、低いのか、そんなもんなのか、、

ということで長くなりましたが、今回はこの辺までにしたいと思います。

次回は将棋、サッカー、競馬など麻雀以外の題材を分析してようかなと思っています。
また、今回は触れなかった直接対決の場合の確率についても考えたいと思います。

ここまで読んでいただいた方は、本当にありがとうございました！！