2017-04-21

StanでAizu Online Judgeの難易度・習熟度を推定したい（３：IRTモデルによる習熟度推定）

Stan 機械学習 Python

シリーズ一覧

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}$ はデータから推定します。なぜ正規分布にしないのかと疑問に思う方もいると思いますが，プログラミングコンテストの世界には外れ値的な人がいるので，正規分布よりも裾の長い分布のほうが適切だと思います。

秋葉拓哉氏の名前を見て「TopCoderでRedCoderらしいからヤバイ」みたいなこと書いてる人いるけど、RedCoderとかそんな次元じゃねーから。
レーティンググラフにするとこんな感じだから。化け物だから。 pic.twitter.com/eWYQ3ZYv5B
— chokudai(高橋直大)🍆🍡🌸 (@chokudai) 2017年3月31日

$\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 の弱情報事前分布をどうやって決めたのか気になると思いますが，最初に q や pf の推定を階層モデルなしで行って，その結果からざっくり決めてみました。いかにも間違った事前分布の決め方のような気がしますが，こうしないと計算が収束しなかったので…

そのあとの 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

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

	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に上げた図とは少し順位が入れ替わっているところもあります。

習熟度 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

*1:https://ja.wikipedia.org/wiki/%E9%A0%85%E7%9B%AE%E5%BF%9C%E7%AD%94%E7%90%86%E8%AB%96

2017-04-17

StanでAizu Online Judgeの難易度・習熟度を推定したい（２：人工データによる実験）

Stan Python 機械学習

シリーズ一覧

kujira16.hateblo.jp

目的

あるユーザがある問題に正答したというデータは得られますが，正答していない問題は，取り組んだけれども実力不足で解けなかったのか，そもそも取り組んでいないのかを区別することができません。データの生成過程についての仮説が正しかったとしても，パラメータの自由度が高すぎてパラメータ推定が行えないかもしれません。

そのため，データの生成過程についての仮説は正しいと仮定して，パラメータ推定が収束するかどうか確かめるために，人工データに対してパラメータが正しく推定できるか実験してみることにしました。

人工データの生成

問題に取り組む確率が $q=0.5$ ，平均的なパフォーマンスが $\mu_\textit{pf}=550$ ，パフォーマンスのばらつきが $\sigma_\textit{pf}=50$ のユーザが1人いて，[100, 1200] の一様分布に従う難易度の問題が50問ある状況を考えます。ユーザは，前述のモデルに従って問題を解き，正答したかどうかを表す情報 $Y$ と，個々の問題の難易度 $D$ が与えられた状況で $q$ , $\mu_\textit{pf}$ , $\sigma_\textit{pf}$ が推定できるかどうか実験します。

データは以下のように生成しました。

scipy.random.seed(0)
N = 50
lo = 100
hi = 1200
q = 0.5
pf_mu = 550
pf_sigma = 50
uniform = scipy.stats.uniform(loc=lo, scale=hi-lo)
D = []
for i in range(N):
    D.append(uniform.rvs())
Y = []
bern = scipy.stats.bernoulli(q)
norm = scipy.stats.norm(loc=pf_mu, scale=pf_sigma)
for i in range(N):
    y = (bern.rvs() == 1) and (norm.rvs() > D[i])
    Y.append(y)
df = pandas.DataFrame(dict(D=D, Y=Y))

パラメータ推定

Stanのモデルコードは以下のとおりです。ユーザが問題に挑戦したかどうかを表す変数 try は隠れ変数ですが，ハミルトニアン・モンテカルロ法（およびその拡張であるNo-U-Turn Sampler）では離散値を扱うことができないので log_sum_exp を使って周辺化しています。また $\mathrm{Normal}(\mu_\textit{pf}, \sigma_\textit{pf})$ に従う乱数が $D[i]$ を上回る or 下回る確率を表すためにCDF（累積分布関数）やCCDF（相補累積分布関数）を使っています。

data {
    int N;
    int<lower=0,upper=1> Y[N];
    real<lower=100,upper=1200> D[N];
}

parameters {
    real<lower=0,upper=1> q;
    real pf_mu;
    real<lower=0> pf_sigma;
}

model {
    for (i in 1:N) {
        if (Y[i] == 1) {
            target += bernoulli_lpmf(1 | q) + normal_lccdf(D[i] | pf_mu, pf_sigma);
        } else {
            target += log_sum_exp(
                bernoulli_lpmf(0 | q),
                bernoulli_lpmf(1 | q) + normal_lcdf(D[i] | pf_mu, pf_sigma)
            );
        }
    }
}

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

f:id:kujira16:20170417211642p:plain

f:id:kujira16:20170417211648p:plain

これでは全然だめですｗ

とりあえず $\mu_\textit{pf}$ に $\mathrm{Normal}(650, 550)$ , $\sigma_\textit{pf}$ に $\mathrm{Cauchy}(0, 100)$ という事前分布を設定することにしました。 $\mu_\textit{pf}$ は [100, 1200] の中にあるだろうという事前知識， $\sigma_\textit{pf}$ はたかだか100程度だろう（実力よりも難しい問題を解けるのはせいぜい+100くらいの難易度の問題までだろう）という事前知識を反映したものです。

data {
    int N;
    int<lower=0,upper=1> Y[N];
    real<lower=100,upper=1200> D[N];
}

parameters {
    real<lower=0,upper=1> q;
    real pf_mu;
    real<lower=0> pf_sigma;
}

model {
    pf_mu ~ normal(650, 550);
    pf_sigma ~ cauchy(0, 100);
    for (i in 1:N) {
        if (Y[i] == 1) {
            target += bernoulli_lpmf(1 | q) + normal_lccdf(D[i] | pf_mu, pf_sigma);
        } else {
            target += log_sum_exp(
                bernoulli_lpmf(0 | q),
                bernoulli_lpmf(1 | q) + normal_lcdf(D[i] | pf_mu, pf_sigma)
            );
        }
    }
}

結果は以下のようになりました。ばっちり収束しています。

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
q	0.658	0.003	0.146	0.359	0.557	0.662	0.761	0.931	1827.000	1.001
pf_mu	460.214	1.955	71.291	288.970	424.561	472.218	506.573	574.057	1330.000	1.001
pf_sigma	83.949	1.634	60.969	13.590	39.530	67.979	112.308	242.430	1392.000	1.000

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

f:id:kujira16:20170417210819p:plain

f:id:kujira16:20170417210827p:plain

f:id:kujira16:20170417210833p:plain

95%信用区間ですが， $q$ は [0.359, 0.931], $\mu_\textit{pf}$ は [289, 574], $\sigma_\textit{pf}$ は [13.6, 242] となりました。MAP推定値はそれぞれ0.674, 491, 43.6となりました。う〜ん…？

ソースコード

Jupyter Notebookは以下の場所に置いています。

https://gist.github.com/arosh/8638f2772e8b6c60d61acaa9a8064490

stanutilは私がStanを使うときのための便利関数群です。traceplotやMAP推定などが入っています。以下の場所に置いていますが，いろいろ試行錯誤している途中なのでbreaking changeが入る予定です。

https://gist.github.com/arosh/8871c22030cdebab6bebada3f3741bfc

続き

kujira16.hateblo.jp

2017-04-17

StanでAizu Online Judgeの難易度・習熟度を推定したい（１：モデル式）

Stan Python 機械学習

シリーズ一覧

kujira16.hateblo.jp

はじめに

Aizu Online Judge (AOJ) という競技プログラミングの練習サイトがあります。

AIZU ONLINE JUDGE: Programming Challenge

解けるか解けないかくらいのちょうど良い難易度の問題に取り組むことが一番成長につながるはずなので，それぞれ問題の難易度やそれぞれのユーザの習熟度を定量的に評価して推薦システムを作りたい気分になります。AOJの各ユーザのページにはこれまでの正答数やレーティングが表示されていますが，この2つの値はユーザの習熟度とは必ずしも一致しません。簡単な問題ばかり解いていれば正答数は伸ばせますし*1，レーティング*2は回答者が少なくて簡単な問題を中心に解いていればかなり大きく伸びるためです。そのため，問題の難易度やそれぞれのユーザの習熟度を計算する気の利いた方法が必要です。

得られる情報

AOJのAPIを利用することで，誰がどの問題に正答したかという情報を得ることができます。

問題の難易度の情報はAOJでは提供されていませんが，AOJ-ICPC というサイトで有志によって一部の問題に対して人手で難易度が付与されています。

難しいところ

正答していない問題は，取り組んだけれども実力不足で解けなかったのか，そもそも取り組んでいなかったのかを区別することができません。そのため，あるユーザが問題に取り組むかどうかを表す確率をモデルに組み込むことにしました。

モデル式

考えたモデルは以下のとおりです。

$\displaystyle q[i] \sim \mathrm{Beta}(a_0,b_0)$

$\displaystyle \mathit{try}[i,j] \sim \mathrm{Bernoulli}(q[i])$

$\displaystyle \mathit{pf}[i,j] \sim \mathrm{Normal}(\mu_\mathit{pf}[i], \sigma_\mathit{pf})$

$\displaystyle Y[i,j]=\begin{cases}1 & \text{if $\mathit{try}[i,j]=1$ and $\mathit{pf}[i,j] \geq D[j]$} \\ 0 & \text{if $\mathit{try}[i,j]=0$ or $\mathit{pf}[i,j] < D[j]$}\end{cases}$

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

$\mathit{pf}[i,j]$ は人 $i$ が問題 $j$ に取り組むときに発揮されるパフォーマンスです。人 $i$ が持っている平均的なパフォーマンス $\mu_\mathit{pf}[i]$ に標準偏差 $\sigma_\mathit{pf}$ のばらつきが乗ったものとしました。パフォーマンスのばらつき $\sigma_\mathit{pf}$ には個人差はないものとしています。

$Y[i,j]$ は人 $i$ が問題 $j$ に正答したかどうかを表す観測可能な変数です。問題に取り組む気分になり，かつ問題に取り組んだときにユーザが発揮したパフォーマンス $\mathit{pf}[i,j]$ が問題の難易度 $D[j]$ を上回った場合に1になります。問題に取り組む気分にならなかったり，ユーザが発揮したパフォーマンスが問題の難易度を下回ったりした場合には0になります。問題の難易度 $D[j]$ は一部の問題についてはAOJ-ICPCで有志によって人手で値が付与されています。