読者です 読者をやめる 読者になる 読者になる

StanでAizu Online Judgeの難易度・習熟度を推定したい(2:人工データによる実験)

シリーズ一覧

kujira16.hateblo.jp

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(相補累積分布関数)を使っています。また今回はやっていませんが,Y[i]=0 の問題と Y[i]=1 の問題でそれぞれデータをまとめればベクトル化することもできます。

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

これでは全然だめですw

とりあえず \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