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}$

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

$\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で有志によって人手で値が付与されています。