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

シリーズ一覧

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