StanでAizu Online Judgeの難易度・習熟度を推定したい(3:IRTモデルによる習熟度推定)

シリーズ一覧

kujira16.hateblo.jp

kujira16.hateblo.jp

はじめに

前回までの記事を公開したところ,Twitterで「問題に取り組んだときの正答確率の部分を項目応答理論でモデリングしないのはなぜか」というコメントをいただきました。

…すいません,項目応答理論というものを知りませんでした。

指摘を頂いてから勉強したのですが,この方法でモデリングするほうが自然だと感じたので,これからは正答確率の部分を項目応答理論でモデリングしていきたいと思います。

モデル式

項目応答理論の1パラメータロジスティックモデルでは,習熟度 \theta の人が 難易度 b_i の問題に正答する確率 p_i(\theta) を以下のようにモデリングします*1

\displaystyle
p_i(\theta)=\mathrm{InvLogit}(\theta-b_i)

ここで \mathrm{InvLogit}(x) はロジスティック関数 1/(1+\exp(-x)) です。

項目応答理論では被験者の習熟度と問題の難易度を同時に推定しますが,今回の記事で使うデータでは一部の問題については難易度が既に付与されているという点が異なります。難易度のスケールによって習熟度のスケールも自動的に決まってしまうので,ロジスティック関数のスケールに自由度を与えるパラメータがあるほうが適切でしょう。そのため,今回は以下のようなモデル式を考えました。

\displaystyle
q[i] \sim \mathrm{Beta}(a_0,b_0)
\displaystyle
\mathit{try}[i,j] \sim \mathrm{Bernoulli}(q[i])
\displaystyle
\mathit{pf}[i] \sim \mathrm{Cauchy}(\mu_\mathit{pf},\sigma_\mathit{pf})
\displaystyle
\mathit{solve}[i,j] \sim \mathrm{Bernoulli}(\mathrm{InvLogit}(\gamma(\mathit{pf}[i]-D[j]))
\displaystyle
Y[i,j]=
\begin{cases}1 & \text{if $\mathit{try}[i,j]=1$ and $\mathit{solve}[i,j]=1$} \\ 0 & \text{if $\mathit{try}[i,j]=0$ or $\mathit{solve}[i,j] =0$}\end{cases}

q[i] は人 i が個々の問題の取り組む確率です。事前分布としてベータ分布 \mathrm{Beta}(a_0,b_0) を設定します。a_0, b_0 はデータから推定します。

\mathit{try}[i,j] は人 i が問題 j に取り組むときに1,取り組まないときに0をとります。

\mathit{pf}[i] は人 i の習熟度で,コーシー分布 \mathrm{Cauchy}(\mu_\mathit{pf},\sigma_\mathit{pf}) に従って生成されると仮定します。位置母数 \mu_\mathit{pf} と 尺度母数 \sigma_\mathit{pf} はデータから推定します。なぜ正規分布にしないのかと疑問に思う方もいると思いますが,プログラミングコンテストの世界には外れ値的な人がいるので,正規分布よりも裾の長い分布のほうが適切だと思います。

\mathit{solve}[i,j] は人 i が問題 j を解けるとき1,解けないとき0をとります。\mathit{pf}[i] は人 i の習熟度,D[j] は問題の難易度です。1パラメータロジスティックモデルを真似たものですが,前述の通り,問題の難易度 D[j] は与えられており,D[j] のスケールから \mathit{pf}[i] のスケールが自動的に決まってしまうため,ロジスティック関数のスケールに自由度を与えるパラメータ \gamma を付けています。\mathit{pf}[i]=D[j] のとき問題を解ける確率は50%です。

Y[i,j] は 人 i が問題 j に正答したかどうかを表す観測可能な変数です。問題に取り組む気分になり,かつ問題に取り組んだときにユーザが問題を解けるならば Y[i,j]=1 になります。問題に取り組む気分にならなかったり,ユーザが問題を解けないならば0になります。

実験内容

問題の難易度とユーザの習熟度の同時推定は現状うまくいっていないので,AOJ-ICPCで難易度が付与されている問題のデータを使ってユーザの習熟度の推定をやった結果を紹介します。

作成したデータセットは以下のURLにあります。

https://github.com/arosh/performance-estimation/blob/master/data.csv

本当はたくさんのユーザの習熟度の推定を行いたかったのですが,問題をある程度たくさん解いている人でないと結果が収束しなかったので,AOJ-ICPCに掲載されている問題をたくさん解いている人上位200人の習熟度を推定します。

Stanのモデルコードは以下のようになります。

data {
    int N;
    int L;
    vector[N] D;
    int G[N, L];
}
parameters {
    vector<lower=0,upper=1>[L] q;
    real<lower=0> a0;
    real<lower=0> b0;
    vector[L] pf;
    real mu_pf;
    real<lower=0> sigma_pf;
    real<lower=0> gamma;
}
model {
    q ~ beta(a0, b0);
    a0 ~ cauchy(0, 2.0);
    b0 ~ cauchy(0, 0.64);
    pf ~ cauchy(mu_pf, sigma_pf);
    mu_pf ~ cauchy(0.455, 0.025);
    sigma_pf ~ cauchy(0, 0.14);
    gamma ~ cauchy(13.3, 0.03);
    for (i in 1:N) {
        for (j in 1:L) {
            if (G[i,j] == 1) {
                target += bernoulli_lpmf(1 | q[j]) + bernoulli_lpmf(1 | inv_logit(gamma * (pf[j] - D[i])));
            } else {
                target += log_sum_exp(
                    bernoulli_lpmf(0 | q[j]),
                    bernoulli_lpmf(1 | q[j]) + bernoulli_lpmf(0 | inv_logit(gamma * (pf[j] - D[i])))
                );
            }
        }
    }
}

a0, b0, pf, pf_mu, sigma_mu の弱情報事前分布をどうやって決めたのか気になると思いますが,最初に qpf の推定を階層モデルなしで行って,その結果からざっくり決めてみました。いかにも間違った事前分布の決め方のような気がしますが,こうしないと計算が収束しなかったので…

そのあとの target += ... と書いている部分ですが,シンプルに G[i,j] ~ bernoulli(q[j] * inv_logit(...)) とも書けてしまいそうです。しかしながら q[j] * inv_logit(...) の部分でアンダーフローしてしまうらしく,予備実験での結果が芳しくありませんでした。これを回避するために数値計算的に安定と思われる方法で冗長に書いています。

このStanのモデルコードをキックするPythonのコードは以下のようになります。stanutilは私がまとめている便利関数群です(第二回の記事を参照)。difficulty は [100, 1200] の値をとりますが,パラメータのスケールがなるべく同じになるようにしたほうが早く計算が収束するので1000で割ってスケールを調整しています。

import pandas
import pystan
import stanutil
L = 200
df = pandas.read_csv('data.csv', index_col=0)

model_code = '''

'''

data = {}
data['N'] = len(df)
data['L'] = L
data['D'] = df['difficulty'] / 1000
data['G'] = df.iloc[:, 1:L+1]
init = lambda: {
    'q': [0.7] * L,
    'a0': 1.6,
    'b0': 0.55,
    'pf': [0.455] * L,
    'mu_pf': 0.455,
    'sigma_pf': 0.115,
    'gamma': 13.3,
}
stan_model = stanutil.stan_cache(model_code)
fit = stan_model.sampling(data=data, seed=0, init=init)

結果

結果の分析に使ったJupyter Notebookは以下の場所に置いています。

https://github.com/arosh/performance-estimation/blob/master/plot.ipynb

サンプリングの結果は以下のとおりです。qpf は量が多いので,気になる人は ここ を見てください。

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a0 1.401 0.005 0.166 1.105 1.285 1.392 1.506 1.760 1100.000 1.003
b0 0.447 0.002 0.049 0.355 0.413 0.445 0.479 0.546 603.000 1.006
mu_pf 0.454 0.000 0.011 0.432 0.447 0.454 0.462 0.478 1741.000 1.002
sigma_pf 0.116 0.000 0.011 0.095 0.108 0.116 0.124 0.139 1456.000 1.001
gamma 13.213 0.005 0.119 12.906 13.161 13.256 13.296 13.349 641.000 1.007

traceplotは次のようになりました。

f:id:kujira16:20170421220849p:plain

f:id:kujira16:20170421220856p:plain

f:id:kujira16:20170421220903p:plain

f:id:kujira16:20170421220909p:plain

f:id:kujira16:20170421220914p:plain

もうTwitterで見た人もいるかと思いますが,pf[i] のサンプリング結果をアカウント名とともにプロットしたものが以下の図です。順番は pf[i] の分布のMAP推定値によって降順に並べています。箱ひげ図の箱の両端は50%ベイズ信頼区間,棒の両端は95%ベイズ信頼区間です。計算をやり直したのでTwitterに上げた図とは少し順位が入れ替わっているところもあります。

https://github.com/arosh/performance-estimation/raw/master/performance.png

習熟度 pf[i] と 問題に取り組む確率 q[i] の分布は次のようになりました。オレンジの線はpf_mu, pf_sigma, a0, b0 のMAP推定値を使ったコーシー分布とベータ分布です。

f:id:kujira16:20170421221942p:plain

f:id:kujira16:20170421221946p:plain

解釈

正答していない問題は,取り組んだけれども実力不足で解けなかったのか,そもそも取り組んでいないのかを区別することができない,という課題がある中で習熟度の推定を行っているので,AOJのsolved rankingよりもマシな指標になっているかどうかが気になるところです。以下の図は,正答した問題数(ただしAOJ-ICPCに掲載されている問題のみ)を横軸,推定した習熟度のMAP推定値を縦軸にとりプロットしたものです。これを眺めた感じでは,強い人として認知されている人はたくさん問題を解いていなくても習熟度が高くなっており,問題に取り組む確率 q[i] を導入した効果が現れていると考えられます。

f:id:kujira16:20170421221153p:plain

しかしながら,前の節に載せた箱ひげ図において中位以下の人の習熟度の推定結果は疑問に思われます。図の横軸の値×1000の難易度の問題を50%の確率で解くことができることになっているのですが,中位以下の人の名前を見てみると,もっと難しい難易度の問題を簡単に解いているように思えてならない人が多々見受けられます。このような結果になった原因として,問題の取り組む確率はどの問題でも独立同一であることを仮定していますが,実際にはこの仮定は成り立たず,簡単な問題から順番に埋めていくような取り組み方をしている人が多いためだと考えられます。

最後にロジスティック関数のスケールパラメータ gamma について言及しておきます。gamma のMAP推定値は13.3でした。問題を解ける確率と解けない確率のオッズ p/(1-p)\exp(\gamma(\mathit{pf}[i]-D[j])) であることから,p/(1-p)=\exp(13.3 \times 100 / 1000)p について解くと p=0.79 くらいになり,難易度が100上がると解ける確率が79%に減少するということが分かります。これはかなり疑わしい値で,実際には難易度が100上がると解ける確率は半分以下になるというのが私の感覚なのですがいかがでしょうか? gamma の事前分布にこの感覚を反映するべきかもしれません。大嘘を書いていました。p/(1-p)=\exp(13.3 \times -100 / 1000)p について解くと0.21くらいになり,難易度が100上がると解ける確率が21%に減少するということが分かります。これなら妥当な結果ですね。

追記

[twitter:@berobero11] さんが追加解析をしてくれました。ありがとうございます!

statmodeling.hatenablog.com