【麻雀】安定段位のバイアスについて
はじめに
天鳳・雀魂といったネット麻雀の世界では、プレイヤーの実力を表す指標のひとつとして安定段位が利用されています。安定段位は長期的には何段で安定できるかを表す指標で、たとえば天鳳の安定段位が八段なら、鳳凰卓で数多く打つと八段で安定すること意味します。安定段位は「段位」として解釈できるというわかりやすさから、最も利用されている実力指標のひとつになっています。
さて、そんな安定段位ですが、実は、試合数が少ないと平均的には真の安定段位（長期的に収束する安定段位）よりも高めの数字が出ることがわかっています。具体的に、トップからラスまで全て25%のプレイヤーがいたとします。このプレイヤーの安定段位は理論的には七段ちょうどになるはずです。しかし、試合数を100戦から2000戦まで100戦刻みで5万回ずつシミュレーションした結果は下図になります。
青い線は安定段位の理論値で、ちょうど七段になっています。オレンジの線は戦績から計算した安定段位です。たとえば100試合の場合、「トップからラスまで確率25%まで取る設定で100試合行い、実現したシミュレーション結果から安定段位を計算する」ということを10万回行い、その結果を平均しています。それを100試合の場合、200試合の場合、300試合の場合、、、と最終的には2000試合の場合まで試して、結果を可視化しています。
ここからわかることは、データから計算された安定段位、つまり僕たちが普段やっているように戦績から普通に計算した安定段位は、理論的な安定段位より平均的には高めに出るということです。そして、その傾向は試合数が少ないほど強くなります。試合数が多くなるとデータから計算された安定段位は真の安定段位に近づいていきますが、2000戦だとまだ若干の乖離が残っています。
この記事では、安定段位のバイアスはどこから生まれているのかを調べ、それ補正する手法を提案しています。提案手法を用いることで、少ない試合数でもバイアスのかからない安定段位の推定が可能になります（ただし、取り除けるのはあくまで「バイアス」だけです。少ない試合数から計算された安定段位は運不運による「ブレ」があり、それは取り除くことができません）。
詳細は以降で見ていきますが、安定段位の計算式を既存の
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\widehat{安定段位} = 7 + \frac{1}{15}\paren{\frac{90\hp_1 + 45\hp_2}{\hp_4} - 135}
$$
から
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\widetilde{安定段位} = 7 + \frac{1}{15}\paren{\frac{90\hp_1 + 45\hp_2}{\hp_4} - \frac{90\hp_1\hp_4 + 45\hp_2\hp_4}{(N-1)\hp_4^2} - \frac{(90\hp_1 + 45\hp_2)(1 - \hp_4)}{(N-1)\hp_4^2} - 135}
$$
に変更することで、バイアスを取り除くことができます。
安定段位のバイアスはどこから来ている？
安定段位の定義
安定段位のバイアスについて考えるため、まずは安定段位の数式を確認します。安定段位は順位率$${(p_1, p_2, p_3, p_4)}$$をとった場合に、長期的には何段で安定するかを計算する指標です。たとえば、鳳南の順位点の配分はトップから順に$${(90, 45, 0, -135)}$$なので、順位点の期待値がちょうどゼロになる順位率、つまり
$$
90p_1 + 45p_2 + 0p_3 - 135p_4 = 0
$$
を満たすような順位率$${(p_1, p_2, p_3, p_4)}$$なら七段で安定することになります。八段だとラスったときのマイナスが15ポイント増えて150になるので、八段で安定するためには
$$
90p_1 + 45p_2 + 0p_3 - 150p_4 = 0
$$
を満たすような順位率をとる必要があります。段位がひとつ上がるとラスったときのマイナスが15ポイントずつ増えることに注意しながらこれをより一般化すると、$${D}$$段で安定するためには、
$$
90p_1 + 45p_2 + 0p_3 - (135 + 15(D - 7))p_4 = 0
$$
を満たすような順位率をとる必要があると考えられます。これを$${D}$$について解いたものが安定段位の計算式です。
$$
\newcommand{\paren}[1]{\left(#1\right)}
安定段位 = D = 7 + \frac{1}{15}\paren{\frac{90p_1 + 45p_2}{p_4} - 135}
$$
3着率の部分$${0p_3}$$は0なので表記を単純にするため消しました。さらに約分して式変形していくこともできますが、順位点の数字がそのまま残っていたほうが解釈が容易だと思うのでここまでにしておきます。また、ラス率$${p_4}$$が0のときは計算できませんが、現実的にはそういう状況はないので気にしなくていいでしょう。
データから安定段位を推定するには？
さて、理論的な安定段位は前述の式で求めることができますが、実際はプレイヤーの真の順位率を知ることはできません。なので、データ（対戦の結果）から順位率を推定することになります。推定は単純で、$${N}$$回打った結果が1着から順に$${(N_1, N_2, N_3, N_4)}$$回なら、順位率は
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\paren{\hp_1, \hp_2, \hp_3, \hp_4} = \paren{\frac{N_1}{N}, \frac{N_2}{N}, \frac{N_3}{N}, \frac{N_4}{N}}
$$
で推定します。つまり、単純に割合を使います。1000回打ってトップから300回、250回、250回、200回なら30%、25%、25%、20%という感じです。当たり前ですね。
普段、安定段位を求めるときは、このデータから推定された順位率を使って安定段を推定しています。
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\widehat{安定段位} = 7 + \frac{1}{15}\paren{\frac{90\hp_1 + 45\hp_2}{\hp_4} - 135}
$$
安定段位の推定量は不偏性を満たさない
それでは、この安定段位の推定はうまくいくのでしょうか？「はじめに」で見たとおり、この推定は（特に試合数の少ないときに）平均的には安定段位を高めに見積もるという問題を抱えています。つまり、推定にバイアスがかかっています。
バイアスが掛かっていない状態、つまり推定量の期待値をとったときに真の値と一致する性質を不偏性と呼びます。たとえば、さきほど出てきた順位率の推定量$${(\hat{p}_1, \hat{p}_2, \hat{p}_3, \hat{p}_4)}$$は不偏性を持ちます。
$${i(=1, \dots, N)}$$試合目の着順を確率変数$${X_i}$$で表すとします。本当の雀力、つまり真の順位率は$${(p_1, p_2, p_3, p_4)}$$なので、$${k=1, 2 ,3, 4}$$に対して$${\Pr(X_i = k) = p_k}$$です（ちなみにこれはカテゴリカル分布と呼ばれています）。
真の順位率はわからないので、$${N}$$試合打った結果から順位率を推定します。
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\hp_k = \frac{N_k}{N} = \frac{1}{N}\sum_{i=1}^N \mathbb{1}\{X_i = k\}
$$
ここで、$${ \mathbb{1}\{X_i = k\} }$$は、$${X_i = k}$$なら1をとり、そうでないなら0をとることを意味します。推定量$${\hat{p}_k}$$は確率変数$${X_i}$$から計算されているので、これも確率変数であることに注意してください。
このとき、$${\hat{p}_k}$$の期待値をとると、真の値$${p_k}$$に一致する、つまり不偏性があることが知られています。
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\begin{align*}
\E{\hp_k}
&= \E{ \frac{1}{N}\sum_{i=1}^N \mathbb{1}\{X_i = k\} }\\
&= \frac{1}{N}\sum_{i=1}^N \E{\mathbb{1}\{X_i = k\}} \\
&= \frac{1}{N}\sum_{i=1}^N \paren{1p_k + \sum_{k'\neq k}0p_{k'}} \\
&= \frac{1}{N}\sum_{i=1}^N p_k\\
&=p_k
\end{align*}
$$
ここで、定数は期待値の外に出せる性質と、足し算の期待値は期待値の足し算になる性質を使っています。
若干話がそれますが、例として安定段位と同様に実力を表す指標として平均順位を考えてみます。
$$
平均順位 = 1p_1 + 2p_2 + 3p_3 + 4p_4
$$
平均順位の推定量はデータから順位率を計算した以下になります。
$$
\newcommand{\hp}[0]{\hat{p}}
\widehat{平均順位} = 1\hp_1 + 2\hp_2 + 3\hp_3 + 4\hp_4
$$
これの期待値をとると
$$
\newcommand{\hp}[0]{\hat{p}}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\begin{align*}
\E{\widehat{平均順位}} &= \E{1\hp_1 + 2\hp_2 + 3\hp_3 + 4\hp_4}\\
&= 1\E{\hp_1} + 2\E{\hp_2} + 3\E{\hp_3} + 4\E{\hp_4}\\
&= 1p_1 + 2p_2 + 3p_3 + 4p_4\\
&= 平均順位
\end{align*}
$$
となり、真の値と一致します。つまり、平均順位の推定量にバイアスはなく、不偏性をもつことがわかります。
一方で、安定段位の推定量$${\widehat{安定段位}}$$は、不偏性を満たしません。期待値をとると
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\begin{align*}
\E{\widehat{安定段位}}
&= \E{7 + \frac{1}{15}\paren{\frac{90\hp_1 + 45\hp_2}{\hp_4} - 135}}\\
&= 7 + \frac{1}{15}\paren{\E{\frac{90\hp_1 + 45\hp_2}{\hp_4}} - 135}\\
&\neq 7 + \frac{1}{15}\paren{\frac{\E{90\hp_1 + 45\hp_2}}{\E{\hp_4}} - 135}\\
&= 7 + \frac{1}{15}\paren{\frac{90p_1 + 45p_2}{p_4} - 135}\\
&= 安定段位
\end{align*}
$$
となり、真の安定段位とは一致しないことがわかります（実際、最初のグラフでも、真の安定段位と実測値から推定された安定段位には乖離がありました）。問題が起きているのは3行目です。基本的に、確率変数の割り算の期待値は、期待値のわり算とは一致しません。なので、ここの等式が成り立たず、安定段位の推定量の期待値は真の安定段位と一致しません。つまり、データから計算された安定段位にはバイアスがかかってしまうことがわかります。
安定段位の推定量のバイアスを見積もる
安定段位の推定量にはバイアスがかかることはわかりましたが、それはどの方向にどのくらいの大きさでバイアスがかかるのでしょうか？それがわかれば、なにか補正を行うことで、真の安定段位をよりうまく推定することができそうです。
確率変数の割り算の期待値を直接考えることは難しいので近似を行うことにします。まずはちょっと一般的なケースを考えます。確率変数$${X}$$と$${Y}$$があるとして、$${f(X, Y)}$$を$${(\mathbb{E}[X], \mathbb{E}[Y])}$$の周りでテイラー展開で2階近似すると
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[#1\right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[#1\right]}
\begin{align*}
f(X, Y)
&\approx f(\E{X}, \E{Y}) \\
&+ f_X\paren{\E{X}, \E{Y}}\paren{X - \E{X}} +f_Y\paren{\E{X}, \E{Y}}\paren{Y - \E{Y}} \\
&+\frac{1}{2}f_{XX}(\E{X})\paren{X - \E{X}}^2 \\
&+ f_{XY}\paren{\E{X}, \E{Y}}\paren{X - \E{X}}\paren{Y - \E{Y}}\\
&+ \frac{1}{2}f_{YY}(\E{Y})\paren{Y - \E{Y}}^2
\end{align*}
$$
なので、$${f(X, Y)}$$の期待値も以下で近似できることがわかります。
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[#1\right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[#1\right]}
\begin{align*}
\E{f(X, Y)}
&\approx f(\E{X}, \E{Y}) \\
&+\frac{1}{2}f_{XX}(\E{X})\Var{X}+f_{XY}\paren{\E{X}, \E{Y}}\Cov{X, Y}+\frac{1}{2}f_{YY}(\E{Y})\Var{Y}
\end{align*}
$$
$${\mathbb{E}[X - \mathbb{E}[X]] = 0}$$なので、1階微分の項は消えてしまいます。とりあえず一般的なケースを考えたのですが、今回は$${f(X, Y) = X / Y}$$のケースを考えます。割り算の微分を考えながら置き換えると
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[#1\right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[#1\right]}
\begin{align*}
\E{\frac{X}{Y}}
&\approx \frac{\E{X}}{\E{Y}} - \frac{\Cov{X, Y}}{\E{Y}^2}+\frac{\E{X}}{\E{Y}^3}\Var{Y}
\end{align*}
$$
となります。さらに安定段位について考えると、問題になっている
$$
\newcommand{\hp}[0]{\hat{p}}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\E{\frac{90\hp_1 + 45\hp_2}{\hp_4}}
$$
のバイアスを近似から知りたいので、$${X = 90\hat{p}_1 + 45\hat{p}_2}$$、$${Y = \hat{p}_4}$$で置き換えてあげればよさそうです。
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[#1\right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[#1\right]}
\newcommand{\hp}[0]{\hat{p}}
\begin{align*}
\E{\frac{90\hp_1 + 45\hp_2}{\hp_4}}
&\approx \frac{\E{90\hp_1 + 45\hp_2}}{\E{\hp_4}} - \frac{\Cov{90\hp_1 + 45\hp_2, \hp_4}}{\E{\hp_4}^2}+\frac{\E{90\hp_1 + 45\hp_2}}{\E{\hp_4}^3}\Var{\hp_4}
\end{align*}
$$
これで2階近似が導出できました。
安定段位の推定量を補正する
右辺第一項に注目すると
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\newcommand{\hp}[0]{\hat{p}}
\begin{align*}
\frac{\E{90\hp_1 + 45\hp_2}}{\E{\hp_4}} = \frac{90p_1 + 45p_2}{p_4}
\end{align*}
$$
で、これは本当に求めたい真の安定段位になっています。それに右辺の第二項と第三項が追加されてしまっているので、安定段位の推定量にはバイアスがかかってしまうことがわかります。逆に考えると、右辺の第二項と第三項を$${(90\hat{p}_1+45\hat{p}_2)/\hat{p}_4}$$から引いてあげればバイアスを取り除くことができそうです。
まず、右辺第二項は以下のように変形できます。
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[#1\right]}
\newcommand{\hp}[0]{\hat{p}}
\begin{align*}
-\frac{\Cov{90\hp_1 + 45\hp_2, \hp_4}}{\E{\hp_4}^2} = -\frac{90\Cov{\hp_1,\hp_4} + 45\Cov{\hp_2, \hp_4}}{\E{\hp_4}^2}
\end{align*}
$$
これをデータから推定する必要があります。着順は確率$${(p_1, p_2, p_3, p_4)}$$のカテゴリカル分布に従うので、カテゴリカル分布の共分散の性質を使うと
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\begin{align*}
-\frac{-90\hp_1\hp_4/(N-1) - 45\hp_2\hp_4/(N-1)}{\hp_4^2} = \frac{90\hp_1\hp_4 + 45\hp_2\hp_4}{(N-1)\hp_4^2}
\end{align*}
$$
でデータから推定できます。一応、共分散を不偏推定するために$${N}$$ではなく$${N-1}$$で割っています。
さらに、右辺第三項もデータから推定する必要があります。こちらもカテゴリカル分布の分散の性質から以下で推定できます。
$$
\newcommand{\hp}[0]{\hat{p}}
\begin{align*}
\frac{90\hp_1 + 45\hp_2}{\hp_4^3}\times\frac{\hp_4(1 - \hp_4)}{N-1} = \frac{(90\hp_1 + 45\hp_2)(1 - \hp_4)}{(N-1)\hp_4^2}
\end{align*}
$$
これらをまとめると、以下のような補正を行うことで、安定段位をより正確に推定できることがわかります。
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\begin{align*}
\widetilde{安定段位} = 7 + \frac{1}{15}\paren{\frac{90\hp_1 + 45\hp_2}{\hp_4} - \frac{90\hp_1\hp_4 + 45\hp_2\hp_4}{(N-1)\hp_4^2} - \frac{(90\hp_1 + 45\hp_2)(1 - \hp_4)}{(N-1)\hp_4^2} - 135}
\end{align*}
$$
補正済み安定段位のシミュレーション結果
数式だけだと途中でなにか間違えていることがよくあるで、シミュレーションでも確認していきましょう。最初のグラフに通常の補正なし安定段位だけでなく、補正済みの安定段位を加えた結果が下図になります。
試合数が少なくても、補正によってうまくバイアスを除去できていることが見て取れます。この手法がうまくいっていると言えそうですね！
おわりに
この記事では、安定段位の推定にはバイアスがかかることを示し、2階近似を用いてそれを補正する手法を提案しました。具体的には、安定段位の計算式を既存の
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\widehat{安定段位} = 7 + \frac{1}{15}\paren{\frac{90\hp_1 + 45\hp_2}{\hp_4} - 135}
$$
から
$$
\newcommand{\paren}[1]{\left(#1\right)}
\newcommand{\hp}[0]{\hat{p}}
\widetilde{安定段位} = 7 + \frac{1}{15}\paren{\frac{90\hp_1 + 45\hp_2}{\hp_4} - \frac{90\hp_1\hp_4 + 45\hp_2\hp_4}{(N-1)\hp_4^2} - \frac{(90\hp_1 + 45\hp_2)(1 - \hp_4)}{(N-1)\hp_4^2} - 135}
$$
に変更することで、バイアスを取り除くことができます。
ただし、取り除けるのはあくまで「バイアス」だけです。少ない試合数から計算された安定段位は運不運による「ブレ」があり、それは取り除くことができません。
実は、今回補正により取り除いたバイアスよりもこの「ブレ」のほうが影響は大きいです。ついでにシミュレーションをしてみたところ、100試合の結果から計算された安定段位の標準誤差は2.5程度になります。これはバイアスの0.4程度よりもはるかに大きい水準です。補正でバイアスを取り除いたとしても、少ない試合数の安定段位で実力を語ることは依然として難しいことは注意していただければと思います。
とはいえ、取り除けるバイアスは取り除いたほうがいいと思います。また、シミュレーションで確かめただけですが、補正を行うとどうやら標準偏差も少し小さくなりそうです（100試合だと2.5から2.3に）。なので、この補正は「ブレ」も少し抑えられそうです。
上記の補正式で安定段位を計算してもらえると、真の雀力により近い安定段位が計算できるので、ぜひ使ってみてもらえると嬉しいです！
最後になりますが、この記事で伝えたかったのはバイアスの補正手法そのものというよりは、数式やシミュレーションを使って麻雀について考えていく楽しさです。この記事の内容は間違っているかもしれませんし、良い手法があるかもしれません。教えていただきたいですし、議論したいです。ご質問・ご意見・ご指摘があれば、Twitter（@ippou1000）などなんでもいいので、どんどんコメントしてもらって、議論が深められたら嬉しいです！
シミュレーションコード
一応、シミュレーションに使ったPythonコードを乗っけておきます。コメントとかなくてすいません。こんな感じで、シミュレーションは短いコードで実現できるので、ぜひやってみてください。
import warnings
import japanize_matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
np.random.seed(42)
# いくつか関数を用意
def calc_antei_danni(p):
return 7 + ((90 * p[0] + 45 * p[1]) / p[3] - 135) / 15
def calc_adjusted_antei_danni(p, N):
term1 = (90 * p[0] + 45 * p[1]) / p[3]
term2 = (90 * p[0] * p[3] + 45 * p[1] * p[3]) / (N - 1) / p[3] ** 2
term3 = (90 * p[0] + 45 * p[1]) / (p[3] ** 3) * (p[3] * (1 - p[3]) / (N - 1))
return (term1 - term2 - term3 - 135) / 15 + 7
def simulate_antei_danni(p, N, M):
x = np.random.choice([1, 2, 3, 4], size=(M, N), p=p)
p_hat = np.row_stack([np.mean(x == k, axis=1) for k in [1, 2, 3, 4]])
return calc_antei_danni(p_hat), calc_adjusted_antei_danni(p_hat, N)
# シミュレーション
M = 100000
p = np.array([0.25, 0.25, 0.25, 0.25])
Ns = np.arange(100, 2001, 100)
antei_danni = np.ones((len(Ns), M, 3))
for i, N in enumerate(Ns):
antei_danni[i, :, 0] = calc_antei_danni(p)
antei_danni[i, :, [1, 2]] = simulate_antei_danni(p, N, M)
# データフレームにまとめる
df = pd.DataFrame(
np.column_stack([Ns, antei_danni.mean(axis=1), antei_danni.std(axis=1)[:, [1, 2]]]),
columns=["試合数", "理論値", "推定値", "推定値（補正済み）", "標準偏差", "標準偏差_補正済み"],
)
# 結果の可視化
df_longer = df.melt(
id_vars="試合数",
value_vars=["理論値", "推定値", "推定値（補正済み）"],
var_name="安定段位の種類",
value_name="安定段位",
)
fig, ax = plt.subplots()
sns.lineplot(x="試合数", y="安定段位", hue="安定段位の種類", data=df_longer, ax=ax)
ax.set(title="真の安定段位とデータから推定された安定段位の比較", xticks=np.arange(0, 2001, step=200))この記事が気に入ったらサポートをしてみませんか？