報酬がベルヌーイ分布に従うときのThompson Sampling

報酬がベルヌーイ分布に従うときのThompson Samplingの追試をだいぶ前にやったのですが,あまり受けが良くなかったのでブログで供養をしていきます。

問題設定

Thompson Samplingは多腕バンディット問題に対する方策の1つです。多腕バンディット問題については良い解説記事やスライドがたくさんあるので,興味がある人はググってみてください。以下のようなデモもあります*1

https://arosh.github.io/multi-armed-bandit/

行ったのはChapelleらのNIPS 2011の論文 [1] の中の実験の1つの追試で,報酬が \lambda=0.5 のベルヌーイ分布に従うアームが1本,\lambda=0.4 のベルヌーイ分布に従うアームが99本あるという設定です。

原理

Thompson Sampling では個々のアーム i を引いたときに当たりが出た回数 N_\mathrm{S}(i),外れが出た回数 N_\mathrm{F}(i) を記録しておきます。アームを選ぶときには個々のアームについてベータ分布 \beta\left(N_\mathrm{S}(i)+1, N_\mathrm{F}(i)+1\right) から乱数を発生させ,最も大きい乱数が得られたアームを選択します。

これはベルヌーイ分布のパラメータ \lambdaベイズ推定して \lambda が最大のアームが i である確率をモンテカルロ法で計算するときに乱数サンプルを1個にした場合と考えることができるのですが,詳細は [2] を参照してください。事前分布を一様分布 \beta(1,1) にするかJeffreys分布 \beta(0.5,0.5) にするかは文献によって異なりますが,報酬がベルヌーイ分布に従う場合にはどちらでもあまり関係ないようです。

ソースコード

乱数の種を変更して複数回実験を行いました。以下のプログラムはコマンドラインでシード値を与えるとThompson Samplingに従ってアームを106回引くシミュレーションを行うものです。pure Pythonでは高速化することが難しいのでGoで書きました。乱数ライブラリの leesper/go_rngプログレスバー表示ライブラリの cheggaaa/pb を使っています。最初に go get してから使ってください。

package main

import (
    "flag"
    "fmt"
    "os"

    "github.com/leesper/go_rng"
    "gopkg.in/cheggaaa/pb.v1"
)

var (
    progress bool
    prior    string
    seed     int64
)

func init() {
    flag.BoolVar(&progress, "progress", false, "")
    flag.StringVar(&prior, "prior", "uniform", "")
    flag.Int64Var(&seed, "seed", 0, "seed value for random generater")
}

type Rng struct {
    bernoulli *rng.BernoulliGenerator
    beta      *rng.BetaGenerator
}

func NewRng(seed int64) *Rng {
    return &Rng{
        bernoulli: rng.NewBernoulliGenerator(seed),
        beta:      rng.NewBetaGenerator(seed),
    }
}

func (rng *Rng) Bernoulli(p float64) bool {
    return rng.bernoulli.Bernoulli_P(p)
}

func (rng *Rng) Beta(alpha, beta float64) float64 {
    return rng.beta.Beta(alpha, beta)
}

func bernoulliBandit(rng *Rng, prior string) ([]float64, error) {
    K := 100
    eps := 0.1
    maxLamb := 0.5
    T := 1000000
    lamb := make([]float64, K)
    lamb[0] = maxLamb
    for i := 1; i < K; i++ {
        lamb[i] = maxLamb - eps
    }
    var alpha, beta float64
    switch prior {
    case "uniform":
        alpha = 1
        beta = 1
    case "jeffreys", "reference":
        alpha = 0.5
        beta = 0.5
    default:
        return nil, fmt.Errorf("Unknown prior: %s", prior)
    }

    var bar *pb.ProgressBar
    if progress {
        bar = pb.New(T)
        bar.Output = os.Stderr
        bar.Start()
    }

    regret := make([]float64, T+1)
    regret[0] = 0
    S := make([]int, K)
    F := make([]int, K)

    for t := 0; t < T; t++ {
        lambEstm := make([]float64, K)
        for i := 0; i < K; i++ {
            lambEstm[i] = rng.Beta(float64(S[i])+alpha, float64(F[i])+beta)
        }
        argMax := 0
        for i := 1; i < K; i++ {
            if lambEstm[i] > lambEstm[argMax] {
                argMax = i
            }
        }
        b := rng.Bernoulli(lamb[argMax])
        if b {
            S[argMax]++
        } else {
            F[argMax]++
        }
        regret[t+1] = regret[t] + lamb[0] - lamb[argMax]
        if progress {
            bar.Increment()
        }
    }
    return regret, nil
}

func main() {
    flag.Parse()
    filename := fmt.Sprintf("%s/%s-seed%03d.txt", prior, prior, seed)
    //計算してからファイル書き込みに失敗したら悲しいので,先にオープンしておく
    w, err := os.Create(filename)
    if err != nil {
        panic(err)
    }
    defer w.Close()
    rng := NewRng(seed)
    regret, err := bernoulliBandit(rng, prior)
    if err != nil {
        panic(err)
    }
    for i, x := range regret[1:] {
        fmt.Fprintln(w, i+1, x)
    }
}

シミュレーションは重たいのでGNU Parallelで並列に実行しました。

kujira16.hateblo.jp

seq 1 200 | parallel --progress "go run chapelle2011.go --seed {} --prior uniform"

リグレットの変化のグラフを描画するためのPythonコードは以下のとおりです。

import seaborn
import scipy
import matplotlib.pyplot as plt
from tqdm import tqdm

seaborn.reset_orig()
seaborn.set_context('talk')
seaborn.set_color_codes()

def load_grid(prefix, st, en, size=1000000):
    A = scipy.empty((size, en - st + 1))
    for i in tqdm(range(st, en + 1)):
        filename = prefix + '{:0>3}'.format(i) + '.txt'
        A[:, i - st] = scipy.loadtxt(filename)[:, 1]
    return A

A = load_grid('uniform/uniform-seed', 1, 200)

label = 'thompson'
ys = scipy.mean(A, axis=1)
del A
xs = scipy.arange(1, ys.shape[0] + 1)
plt.plot(xs, ys, label=label)

ys = scipy.log(xs)
p = 0.5
q = 0.4
coef = (p-q)/(q*scipy.log(q/p)+(1-q)*scipy.log((1-q)/p)) * 99
plt.plot(xs, coef*ys - coef*scipy.log(100), color='k', label='asymptotic bound')

plt.xlabel('plays')
plt.ylabel('regret')
plt.legend(loc='upper left')

plt.xlim(xmin=100)
plt.ylim(0, 10000)
plt.xscale('log')

plt.savefig('chapelle2011.png', dpi=100)

結果

追試で欲しかった図は [1] p.3のFigure 1の右上の図です。asymptotic boundはリグレットの漸近下界で,カルバック・ライブラー・ダイバージェンスを使って式 (2) から計算しています。

f:id:kujira16:20170501140454p:plain

続き

kujira16.hateblo.jp

参考文献

  • [1] Chapelle, O. and Li, L.: An Empirical Evaluation of Thompson Sampling, in NIPS 2011 [PDF].
  • [2] 本多淳也,中村篤祥:バンディット問題の理論とアルゴリズム講談社(2016).

バンディット問題の理論とアルゴリズム (機械学習プロフェッショナルシリーズ)

バンディット問題の理論とアルゴリズム (機械学習プロフェッショナルシリーズ)

*1:無駄にReact+Reduxで作りました

StanでAizu Online Judgeの難易度・習熟度を推定したい(3:IRTモデルによる習熟度推定)

シリーズ一覧

kujira16.hateblo.jp

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

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

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

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

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

https://github.com/arosh/performance-estimation/raw/master/performance.png

習熟度 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%に減少するということが分かります。これなら妥当な結果ですね。

追記

@ さんが追加解析をしてくれました。ありがとうございます!

statmodeling.hatenablog.com

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(相補累積分布関数)を使っています。

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