load hw5_p1

numexp = 500;
nsamples = [42 62 92 135 198 292 431 635 936 1380 2035 3000];
small_data = marginalize(data, [1 2 3 4]);

bayesian_score_diff = zeros(size(nsamples));
bic_score_diff = zeros(size(nsamples));
for i = 1:length(nsamples)
    score1 = zeros(1,numexp);
    score2 = zeros(1,numexp);
    for e = 1:numexp
        sample = subsample_data(small_data, nsamples(i));
        score1(e) = score_pst(pst_extended, sample, @log_bayesian_model_score) - ...
            score_pst(pst_original, sample, @log_bayesian_model_score);
        score2(e) = score_pst(pst_extended, sample, @bic_model_score) - ...
            score_pst(pst_original, sample, @bic_model_score);
    end
    bayesian_score_diff(i) = mean(score1);
    bic_score_diff(i) = mean(score2);
end

plot(nsamples, bayesian_score_diff./nsamples, '-k', ...
    nsamples, bic_score_diff./nsamples, '-.k', ...
    nsamples, zeros(length(nsamples)), ':k');
xlabel 'number of samples'
legend('log Bayesian score difference per sample','BIC score difference per sample',4)
title 'Comparison of the model selection decision based on the Bayesian and BIC scores'
    