コンテスト中に問題のdifficultyを推定した(かった)

～前回までのあらすじ～

コンテスト中に最終順位を推定したくなり、それをするにはコンテスト中にDiffを推定すればよいことが分かった。

keymoon.hatenablog.com

先行研究

コンテスト中にdiffを推定するということ自体は、すでにいくつかの先行研究がある。本項では、DEIM2020にて発表のあった「競技プログラミングコンテストにおけるタスクの難易度のリアルタイム推定*1」という論文で取られている手法の一つを紹介し、それの改善点について考察をする。

手法

まず、正解率pをレートとdifficultyの関係より導出される式 p=\frac{1}{1+6^{(d-r)/400}} より定めた。
そして、「時間T毎に各参加者が提出を行い、確率p_iで正答を得る。また、コンテスト参加者があるタスクに自分の解答を繰り返し提出し，正答であるという判定を得るまでにかかる時間が指数分布に従う。」というモデルを基に推定を行っていた。この場合の尤度関数L(d)は、指数分布の確率密度関数を用いて
\displaystyle L(d)=\prod_{i=1}^{N}{\frac{p_i}{T}\exp{\left( -\frac{p_i}{T}t_i \right)}}
と定められる。これの最大化については、尤度関数の対数を取り
\displaystyle \begin{eqnarray*} \log L(d)&=&\sum_{i=1}^{N}{\log(\frac{p_i}{T} \exp(-\frac{p_i}{T}t_i))} \\ &=&\sum_{i=1}^{N}{\left( \log\frac{p_i}{T}-\frac{p_i}{T}t_i \right)} \\ &=&\sum_{i=1}^{N}{\left( \log p_i-\log T \right)}-\sum_{i=1}^{N}{\frac{p_it_i}{T}} \\ &=&\left( \sum_{i=1}^{N}\log p_i \right) -N\log T -\sum_{i=1}^{N}{\frac{p_it_i}{T}} \cdots (1) \end{eqnarray*}
とした関数を最大化することと同じになる。これは、Tについて偏微分した
\displaystyle \begin{eqnarray*} \frac{\delta \log{L(d)}}{\delta T}=-\frac{N}{T}+\frac{1}{T^2}\sum_{i=1}^{N}{p_it_i} \end{eqnarray*}
が0となるようなTを離散的な各 d\ \mbox{s.t.}\ d \in ℤ, \rm{possibleMinDiff} \leq d \leq \rm{possibleMaxDiff}について求め、それによって\log L(d)を求めることで擬似的に実現できる。
そのようなTは
\displaystyle T=\frac{1}{N}\sum_{i=1}^{N}{p_it_i}
と表せるので、(1)式は
\displaystyle (1)=\left( \sum_{i=1}^{N}\log p_i \right)-N\log T-N \\
と表すことができる。

問題点

このモデルの問題点について考察する。
まず、各人がコンテストを通して解答できる時間の期待値e_iをT/p_iとしているが、それの根拠となる仮定である「時間T毎に提出して確率p_iで正答を得る」という点について。p_iを算出するのに用いているdifficultyの定義は「全体で正答できる確率が50%になるようなレート」であるが、この仮定の下では全体での正答確率がp_iとならないため、これは定義を満たしていない。
また、確率分布が以上の仮定に従うとすると離散的な指数分布を取ると考えられるが、その仮定とは別に指数分布となるモデルを考えているため、複数のモデルが並立していることになっている。
また、レート毎の解答時間の分布に指数分布を用いているが、レート毎の解答時間の分布は対数正規分布になることが観察と考察によって確かめられる。

pepsin-amylase.hatenablog.com

よって、本記事では対数正規分布を解答時間の分布として用いて推定を行ってみたいと思う。

正規分布による推定

difficultyの定義より、レートrの人が難易度dの問題を終了時刻t_{end}において正答できている確率は
\displaystyle \frac{1}{1+6^{(d-r)/400}}\cdots(2)
と表せる。
ここで、あるレートrにおいて正答を得られるまでの時間tは対数正規分布に依り、すべてのtについて対数スケールでの分散が一定であると仮定する。
そのため、あるレートrに対して対数スケールでの解答時間の中央値がe^\mu、同じく分散が\sigma^2のとき、時間tに対する確率密度関数は
\displaystyle \frac{1}{\sqrt{2\pi\sigma^2}t}\exp\left(-\frac{(\log t-\mu)^2}{2\sigma^2}\right)
となる。
ここで、正規分布は尺度s=\frac{2}{\pi}|\sigma|のロジスティック分布に近似できるため、そのような近似を用いて\muをsの式で表すことを考える。この\mu,sは、時間tに対する対数ロジスティック分布の累積分布関数\frac{1}{1+e^{(\mu-\log t)/s}}とレートあたりの正答率の式(2)より、
\displaystyle \begin{eqnarray*} \frac{1}{1+e^{(\mu-\log t_{end})/s}}&=&\frac{1}{1+6^{(d-r)/400}} \\ \Leftrightarrow e^{(\mu-\log t_{end})/s}&=&6^{(d-r)/400} \\ \Leftrightarrow (\mu-\log t_{end})/s&=&(d-r)/400\cdot\log 6=\frac{\log 6}{400} \cdot (d-r) \\ \Leftrightarrow \mu&=&\frac{s\log 6}{400} \cdot (d-r) + \log t_{end} \end{eqnarray*}
という式に従う。

以上より、\muを\sigmaに対する式として表すと、
\displaystyle \mu=\frac{|\sigma|\log 6}{200\pi} \cdot (d-r) + \log t_{end}
となる。ここで、簡便のためにc=\frac{200\pi}{\log 6}とする。

これを時間tに対する確率密度関数に代入すると、
\displaystyle \frac{1}{\sqrt{2\pi\sigma^2}t}\exp\left(-\frac{\left(\log t-\left(\frac{|\sigma|}{c} \cdot (d-r) + \log t_{end}\right)\right)^2}{2\sigma^2}\right) \\
となる。上式より、尤度関数L(d)を
\displaystyle L(d)=\prod_{i=1}^{N}\frac{1}{\sqrt{2\pi\sigma^2}t_i}\exp\left(-\frac{\left(\log t_i-\left(\frac{|\sigma|}{c} \cdot (d-r_i) + \log t_{end}\right)\right)^2}{2\sigma^2}\right) \\
と定める。対数を取ると、
\displaystyle \begin{eqnarray*} \log L(d)&=&\sum_{i=1}^{N}\log\left(\frac{1}{\sqrt{2\pi\sigma^2}t_i}\exp\left(-\frac{\left(\log t_i-\left(\frac{|\sigma|}{c} \cdot (d-r_i) + \log t_{end}\right)\right)^2}{2\sigma^2}\right)\right) \\ &=&\sum_{i=1}^{N}-\log(\sqrt{2\pi}|\sigma|)-\log t_i-\frac{\left(\log t_i-\left(\frac{|\sigma|}{c} \cdot (d-r_i) + \log t_{end}\right)\right)^2}{2\sigma^2} \\ &=&-N\log\sqrt{2\pi}-N\log|\sigma|-\sum_{i=1}^{N}\log t_i-\sum_{i=1}^{N}\frac{\left(\log t_i-\log t_{end}-\frac{|\sigma|}{c} \cdot (d-r_i)\right)^2}{2\sigma^2} \\ &=&-\left(N\log\sqrt{2\pi}+\sum_{i=1}^{N}\log t_i+N\log|\sigma|+\frac{{\displaystyle \sum_{i=1}^{N}}\left(\log t_i-\log t_{end}-\frac{|\sigma|}{c} \cdot (d-r_i)\right)^2}{2\sigma^2}\right) \\ \end{eqnarray*}
となる。\sigmaについて偏微分すると、
\displaystyle \begin{eqnarray*} \frac{\delta\log L(d)}{\delta\sigma}&=&-\frac{N}{\sigma}-\sum_{i=1}^{N}-\frac{(\log t_i-\log t_{end})(c(\log t_i-\log t_{end})-|\sigma|(d-r_i))}{c\sigma^3}\\ &=&-\frac{N}{\sigma}+\sum_{i=1}^{N}\frac{(\log t_i-\log t_{end})(c(\log t_i-\log t_{end})-|\sigma|(d-r_i))}{c\sigma^3}\\ \end{eqnarray*}
となり、上式の偏微分が0になるために\sigmaが満たすべき条件は

\displaystyle \begin{eqnarray*} -\frac{N}{\sigma}+\sum_{i=1}^{N}\frac{(\log t_i-\log t_{end})(c(\log t_i-\log t_{end})-|\sigma|(d-r_i))}{c\sigma^3}&=&0\\ \Leftrightarrow\sum_{i=1}^{N}\frac{(\log t_i-\log t_{end})(c(\log t_i-\log t_{end})-|\sigma|(d-r_i))}{\sigma^2}-cN&=&0\\ \Leftrightarrow\sum_{i=1}^{N}\frac{c(\log t_i-\log t_{end})^2-|\sigma|(d-r_i)(\log t_i-\log t_{end})}{\sigma^2}-cN&=&0\\ \Leftrightarrow\sum_{i=1}^{N}\frac{c(\log t_i-\log t_{end})^2}{\sigma^2}-\sum_{i=1}^{N}\frac{(d-r_i)(\log t_i-\log t_{end})}{|\sigma|}-cN&=&0\\ \Leftrightarrow\left(\sum_{i=1}^{N}c(\log t_i-\log t_{end})^2\right)\frac{1}{\sigma^2}-\displaystyle \left(\sum_{i=1}^{N}(d-r_i)(\log t_i-\log t_{end})\right)\frac{1}{|\sigma|}-cN&=&0\\ \Leftrightarrow cN\sigma^2+\left(\sum_{i=1}^{N}(d-r_i)(\log t_i-\log t_{end})\right)|\sigma|-\sum_{i=1}^{N}c(\log t_i-\log t_{end})^2&=&0\\ \end{eqnarray*}
である。これを満たすような\sigmaを離散的な各d\ \mbox{s.t.}\ d \in ℤ, \rm{minDiff} \leq d \leq \rm{maxDiff}について求め、それぞれについて\log L(d)を求めた後、これが最大値をとるdを擬似的な最尤推定量とする。

\sigmaについての項は\sigma^2と|\sigma|しか出てこないことから、実装においては非負と仮定しても結果は変化しないことを利用すると良い。

実装,性能評価

以上のモデルを既存手法とともにC#にて実装をし、比較をした。実装は以下のリンクよりgistにて閲覧できる。
AtCoderからのデータ取得周りは外部の自作ライブラリに依存してるのでこれ単体では動かないです · GitHub

ABC170 旧モデル 提案モデル difficulty
A -163 -1000 -3425.5
B 14 -792 -872.4
C 233 -414 -166
D 274 323 1033.5
E 625 805 1466.4
F -240 1261 1913.6
f:id:keymoon:20200804132256p:plain
ABC170の推定
ABC171 旧モデル 提案モデル difficulty
A -205 -985 -2837.3
B -77 -1000 -1340.5
C -261 119 466.3
D 38 115 422.3
E -1000 308 733.7
F 614 1035 1721.5
f:id:keymoon:20200804132358p:plain
ABC171の推定
ABC172 旧モデル 提案モデル difficulty
A -167 -1000 -7393.8
B -60 -1000 -1410.4
C 704 347 877.7
D 17 351 944.7
E 710 961 1877.9
F 553 1368 2158.5
f:id:keymoon:20200804132439p:plain
ABC172の推定

旧モデルは収束するのは低難度帯に限られるが収束は早い。一方、提案モデルでは難易度によらず推定を行えているが収束はあまりしているようには見えない。

また、提案モデルでは解がほぼ確実に下方向にズレているのが課題である。そのため、どの程度ズレているのかを確かめてみた。

f:id:keymoon:20200804150254p:plain
推定データと実データ

すると、新ABCに限ってサンプリングしたデータではあるが、実データと推定データが比例しているという傾向が全diffについて見られた。つまり、最終推定値に1次関数による補正をかけたらこのモデルはそこそこな精度になるということである。
ここまで綺麗なズレが観測されるとなると、ここには何らかのモデル自体の欠陥が絡んできているのではないか?と思えてしまう。そのため、ズレがどこから来ているのかということを今後調べていきたいと思っている。

また、ABC/AGC/ARCそれぞれについて同様に調べてみると、グラフのかたむきは0.72程で変わらずに切片のみが変化していることが分かった。

f:id:keymoon:20200804202956p:plain
AGC/ARC/ABCについての推定データと実データ

おわりに

尤度関数と偏微分の式変形がかなり美しい結果に落ち着いたので、かなり満足しています。ただ、実装してみた結果うまく行かなかったので、とりあえず進捗を書くだけ書いておいて本業(競プロ)をしようという不純なモチベーションより筆を取り始めました。しかし、その途中で気まぐれでデータを比較してみるとかなり綺麗なズレが観測できたのでもう一歩な気もしてきました。もう少し頑張ってみます…

*1:https://proceedings-of-deim.github.io/DEIM2020/papers/I4-2.pdf