背景やデータはしょラーさんの以下のブログ記事を読んでください。
kujira16.hateblo.jp
この記事ではAOJ-ICPCで付加された貴重な難易度の情報をフル活用して、問題の真の難易度の推定と、各ユーザの習熟度の推定を行います。
この問題の難しさは「解いていない問題が、スキップして取り組んでいないのか、解こうとしたけど解けなかったのか区別できない」という点にあります。そこで、元記事にもあったように問題をスキップする確率を導入してモデリングする必要があります。
とはいえ、まずはモデルのヒントになりそうなグラフを作成します。
以下では元記事にあわせて、難易度をdifficulty（StanコードではD）、習熟度をperformance（Stanコードではpf）と表現します。
データの分布の確認
difficultyの分布
横軸にdifficulty、縦軸に問題の数をとったヒストグラムは以下になります。山型の分布、途中から100刻みしかない、1000を超えると問題数が減ってくる、などが分かります。
問題の解かれた割合(%)
横軸にdifficulty、縦軸に今回データに含まれる全ユーザー（1000人）のうち何人がその問題を解いたかの割合(%)をとった散布図は以下になります。1つの点が1つの問題です。difficultyがあがると解かれる問題の割合がなんとなくですが指数関数的に落ちていきます。簡単な問題でも解いている人が少ない問題もあります。はじめはこの性質をモデルに取り込もうとしましたが、うまくいきませんでした。
解いた問題のdifficultyの平均と標準偏差
横軸にユーザが解いた問題のdifficultyの平均、縦軸に同じく標準偏差をとった散布図は以下になります。1つの点が1つのユーザです。簡単な問題をスキップしまくって難しい問題だけチャレンジする集団（図だと右下らへん）とかいるのかなと思ってグラフを作りましたが、思ったよりみなさん幅広く解いています。
このように背景知識から仮説を確かめていく過程でグラフを作ることは大変有効です。
このような問題はIRT（Item Response Theory, 項目応答理論）というロジスティック回帰の一種で扱うことが一般的です*1。しかし、個人的にはこのような対戦ムード（問題 vs ユーザ）があるような現象に関してはプロビット回帰を使うのがよいと思っています。理由は「ユーザが対戦相手の力量を上回る（勝利となる）確率」が、累積正規分布で気持ちよく表現できるからです。また、ある説明変数を変えた場合にオッズの観点で議論になりにくいと考えているためです。そこで、この記事ではプロビット回帰（すなわちしょラーさんの2つ目の記事にあるモデルとほぼ同じ）を使います。実用上はロジスティック回帰とプロビット回帰はそこまで差がないと思うのでどちらを使ってもよいと思います。
Stanコードは以下になりました。
- 2～3行目： データに含まれる問題数を
Q（ここでは608）、ユーザ数をN（ここでは1000）で宣言しました。
- 4行目： AOJ-ICPCで付加された難易度です。あとで
1000で割ってスケーリングして渡します。
- 5行目： 元記事では
Gに対応します。問題数×ユーザ数の2次元配列です。あるユーザがある問題を解いている場合に1、その他の場合に0となっています。
- 6行目： ユーザが解いたことがある問題のうち、最も難しかった問題の難易度です。後で使います。
- 16行目： 問題をスキップする確率です。ざっとデータを見るとユーザごとに大きく異なりそうだったので、ユーザごとに宣言します。以降で事前分布は設定しないので一様分布に従います。
- 20行目： 今回は最も難しい問題の難易度が1200で、それ以上のperformanceは正確に測定できない（とても大きな値を取るか分からない）はずなので、ユーザのperformanceは正規分布に従うとしました。この仮定があると、解いた問題数が少ないユーザがいても、
mu_pfあたりに推定してくれて推定が安定します。なお、正規分布の代わりに student_t(6, mu_pf, s_pf) のような少しだけ裾が長い分布を試してもほぼ同様の結果でした。
- 21～22行目： 一応階層モデルにしてあります。
- 15, 23行目：
sigmaはプロビット回帰で使用する累積正規分布の標準偏差です。対戦を扱うプロビット回帰においては、問題のdifficultyとユーザのperformanceに差がある場合に、どれぐらい勝負のアヤがあるかを表していると解釈できます。sigmaが小さいと少しでも差があると強い方が順当に勝つことが多く、sigmaが大きいと差が少々あっても確率的に弱い方が勝つことがあるといった具合です。IRTにおける「識別パラメータ」に相当します。本来はsigmaは問題ごとに推定できるとよいのですが、今回は推定が厳しかったので全問題で共通のsigmaとしました。
- 14, 24行目：
d_trueは問題の真の難易度です。元記事のようにDをそのまま使うことも考えられますが、一般的に人がつけたものはキリのいい数字に偏りやすく、また、誤差を含んでいると考えた方がよいでしょう。そこで、24行目では平均D・標準偏差0.1の正規分布に従うとしています。標準偏差0.1は元の難易度が100ぐらいはブレるかなと考えていることに相当します。なお、標準偏差0.05でも実行してみましたが、そこまで大きな違いはありませんでした。経験ではこのような変数を導入することで、Dのままだとどうしても矛盾してしまうようなところがフニャリと解消されて推定が安定化することが多いです。
- 28～29行目： 問題が解けた場合です。ユーザが問題をスキップしないで解いた（勝利した）確率になります。なお、
log1m(x)はlog(1-x)です。
- 31～32行目： 問題を解いていない場合です。ユーザが解いたことがある最も難しい問題の難易度より
D_range下回っている場合、簡単すぎてつまらないからスキップしているとみなします。
- 33～34行目： 問題を解いていない場合です。ユーザが解いたことがある最も難しい問題の難易度より
D_range上回っている場合、難しすぎてチャレンジしても解けないとみなします。
- 35～39行目： 問題を解いていない場合で、その間の難易度の問題はスキップしたかチャレンジして解けなかったの混合分布になります。
31～39行目の仮説はかなり大胆ですが、このようにデータを稼がないと、「解いていない問題が、スキップして取り組んでいないのか、解こうとしたけど解けなかったのか区別できない」という問題を打破できずにq_skipが不自然な値に収束してしまいます。結局ここでハマってトータル50個ぐらいモデルを試行錯誤しました。このようにある仮説に従ってデータを置き換えるのは統計モデリングでは常套手段で、「StanとRでベイズ統計モデリング (Wonderful R)」では5.3.3項、「予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)」では8.3節で扱っています。なお、if文の中のDをd_trueにすると、if文のなかにパラメータを含む推定が非常に厳しいモデルとなってしまい、今回だと推定できなくなります。
以下はキックするRコードです。
- 20～21行目：
DとD_maxはスケーリングして渡します。D_rangeは仮説および推定結果のsigmaとの兼ね合いになりますが、0.2～0.25ぐらいでsigmaより十分に大きい値となりましたので0.2としました（元の難易度で200）。
- 24行目： ここでは初期値を設定していませんが収束しました。ただし、chainによっては時間がかかる場合があったので初期値を
init=function() { list(pf=max_difficulty/1000, d_true=d_ori$difficulty/1000) } のように定めた方がよいかもしれません。
結果
Solvedの要素数が60万を超えていることもあり、推定に要した時間はおよそ13時間でした。もう少しデータ増えたら、自動変分ベイズであるADVI使った方がよさそうです（もしくはK年後のGPU化Stanを待つ）。
推定されたsigma
推定されたsigmaの値は中央値が0.077、95%ベイズ信頼区間が[0.074, 0.079]でした。2×標準偏差で考えると、元の難易度のスコアにしておおよそ2×77≒150ぐらい差があると、解ける・解けないがはっきりするという解釈になります。与えたD_rangeはそれよりも50ほど大きな値になっています。
ユーザランキング Top 50
推定した習熟度（performance, pf）のMCMCサンプルの中央値Top 50のユーザを示します。
凡例は元記事にあわせてあります。すなわち、ヒゲが95%ベイズ信頼区間、箱が50%ベイズ信頼区間、真ん中の印が中央値です。問題を多く解いている人はベイズ信頼区間がせまくなっているのが分かります。
また、このモデルでは難易度重視（難しい問題が解けるか）でユーザランキングが決まります。例えば、下から13番目のasi1024さんは問題を非常に多く解いているため、現時点でAOJ-ICPCのランキングでは1位です。しかしながら、今回のモデルですと、解いている問題が多い＝q_skipが小さい、それにもかかわらず1000より難易度が高い問題を比較的あまり解いていない＝スキップではなく解けなかった確率が高い、と解釈されてpfが若干小さくなります。すなわち実力ある人が簡単な問題を多めに解くと損になります。「昔から真面目に解いてきたけど、忙しくなって最近出た難しい問題は（解けるにもかかわらず）着手できていない」場合が損になるのを避けたい場合には、問題が発表された日時や、そのユーザのアクティブ日時などをモデルに組み込むと改善する可能性があります。
q_skip vs. pfの比較
横軸にq_skipの中央値、縦軸にpfの中央値をとった散布図は以下になります。Top 20だけラベルを付けました。
AOJ-ICPCによる難易度と推定された難易度の比較
横軸にD、縦軸にd_trueの中央値をとった散布図は以下になります。1つの点が1つの問題です。ちゃんとy=xの直線に載っているのでそこまで大きくは変わらないことが分かります。
AOJ-ICPCによる難易度と推定された難易度が大きく異なるTop 30
横軸・縦軸はひとつ前のグラフと同じです。中央値だけでなく、95%ベイズ信頼区間をヒゲで、50%ベイズ信頼区間で箱を表現しました。
例えば問題番号2710はBNFで数式を定義している問題で、AOJ-ICPCで付加された難易度は400ですが、このモデルでは500超えてるんじゃないかなと推定されています。同様に、問題番号1185のチョコレートの問題は1000という難易度が付加されていますが、このモデルではせいぜい800ぐらいかなと推定されています。
僕はドメイン知識がないので、有識者の解釈を聞いてみたいです。結果はsummaryだけgistにあります。
まとめ
- StanでAizu Online Judgeの難易度・習熟度を推定しました。
- ユーザの力量にあった難易度の問題を、easyモード・normalモード・hardモードなどで推薦できそう。
- 有志によって付加された難易度の情報は非常に貴重で、モデリングの重要な足掛かりとなりました。
- 問題をスキップしたか、チャレンジして解けなかったかをデータとして保持できれば、よりよいモデリングができそう。
- AOJだけでなく、問題が公開されていて自由にチャレンジできるタイプのプログラミング教育サイトに適用できそう。
Enjoy!
