報酬が分散未知の正規分布に従うときのThompson Sampling

前回

kujira16.hateblo.jp

問題設定

行ったのはHondaらのAISTATS 2014の論文 [1] の中の実験の追試で,報酬が \mu=1, \sigma=3正規分布に従うアームと \mu=0, \sigma=0.3正規分布に従うアームがそれぞれ1本ずつあるという設定の多腕バンディット問題です。

原理

最初にパラメータ \alpha を -1/2, 0, 1/2 あたりの値に設定します。\alpha については後述します。

Thompson Samplingを始める前に,まずそれぞれのアームを \max(2,3-\lfloor 2\alpha \rfloor) 回引いておきます。このとき得られた報酬も個々のアームについて記録しておきます。

この後からは,個々のアーム i についてStudentの t 分布 \mathrm{StudentT}(\nu=n_i+2\alpha-1, \mu=\mu_i, \sigma=\sigma_i/\sqrt{n_i+2\alpha-1}) から乱数を発生させ,最も大きい乱数が得られたアームを選択します。ここで n_i はアーム i を引いた回数,\mu_i はアーム i から得られた報酬の標本平均,\sigma_i はアーム i から得られた報酬の標本標準偏差(Nで割るほう)です。

Studentの t 分布はどこから出てくるのかというと,データ D が得られたもとで平均 \mu標準偏差 \sigma が未知の正規分布のパラメータの分布 P(\mu,\sigma \mid D)ベイズの定理から計算し,これを \sigma について周辺化すると出てくるようです [2]。このとき,事前分布 P(\mu,\sigma) の設定方法には諸説あるらしく,一様分布を設定したものが \alpha=-1/2 に,reference priorを設定したものが \alpha=0 に,Jeffreys priorを設定したものが \alpha=1/2 にそれぞれ対応するようです。

前回の実験では事前分布は結果にほとんど影響を及ぼさなかったのですが,今回の実験では違っていて,設定した事前分布によって結果が大きく異なっています。このあたりの難しい話については追えていないので [1] を参照してください。

ソースコード

以下のプログラムはコマンドラインでシード値と事前分布を与えるとThompson Samplingに従ってアームを10^4回引くシミュレーションを行うものです。報酬の標準偏差を計算するときには,得られた報酬の2乗の和を持っておけば \mathbf{Var}[X]=\mathbf{E}[X^2]-\mathbf{E}[X]^2 を使って O(1) で計算できます。

package main

import (
    "bufio"
    "flag"
    "fmt"
    "math"
    "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 {
    gaussian *rng.GaussianGenerator
    studentT *rng.StudentTGenerator
}

func NewRng(seed int64) *Rng {
    return &Rng{
        gaussian: rng.NewGaussianGenerator(seed),
        studentT: rng.NewStudentTGenerator(seed),
    }
}

func (rng *Rng) Gaussian(mu, sigma float64) float64 {
    return rng.gaussian.Gaussian(mu, sigma)
}

func (rng *Rng) StudentT(nu int64, mu, sigma float64) float64 {
    return sigma*rng.studentT.Student(nu) + mu
}

func IntMax(a, b int) int {
    if a < b {
        return b
    }
    return a
}

func ArgMax(xs []float64) int {
    argmax := 0
    for i := 1; i < len(xs); i++ {
        if xs[argmax] < xs[i] {
            argmax = i
        }
    }
    return argmax
}

func gaussianBandit(rng *Rng, prior string) ([]float64, error) {
    K := 2
    mu := []float64{1, 0}
    sigma := []float64{3, 0.3}
    T := 10000

    var alpha float64
    switch prior {
    case "uniform":
        alpha = -0.5
    case "reference":
        alpha = 0
    case "jeffreys":
        alpha = 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

    N := make([]float64, K)
    sum := make([]float64, K)
    sum2 := make([]float64, K)

    C := IntMax(2, 3-int(math.Floor(2*alpha)))
    for k := 0; k < K; k++ {
        for j := 0; j < C; j++ {
            r := rng.Gaussian(mu[k], sigma[k])
            N[k]++
            sum[k] += r
            sum2[k] += r * r
            t := k*C + j
            regret[t+1] = regret[t] + mu[0] - mu[k]
            if progress {
                bar.Increment()
            }
        }
    }

    for t := K * C; t < T; t++ {
        muEstm := make([]float64, K)
        for i := 0; i < K; i++ {
            nu := int64(N[i] + 2*alpha - 1)
            mu := sum[i] / N[i]
            sigma := math.Sqrt((sum2[i]/N[i] - math.Pow(sum[i]/N[i], 2)) / float64(N[i]+2*alpha-1))
            muEstm[i] = rng.StudentT(nu, mu, sigma)
        }
        k := ArgMax(muEstm)
        r := rng.Gaussian(mu[k], sigma[k])
        N[k]++
        sum[k] += r
        sum2[k] += r * r
        regret[t+1] = regret[t] + mu[0] - mu[k]
        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 := gaussianBandit(rng, prior)
    if err != nil {
        panic(err)
    }
    bw := bufio.NewWriter(w)
    defer bw.Flush()
    for i, x := range regret[1:] {
        fmt.Fprintln(bw, i+1, x)
    }
}

シミュレーションは論文 [1] と同じようにシード値を変えて20000回行いました。事前分布とシード値の組合せを変更して実験を行う必要があるので,GNU Parallelがたいへん便利でした。

#!/bin/bash
set -ue
seq 1 20000 > seed.txt
printf "uniform\nreference\njeffreys" > prior.txt
cat prior.txt | xargs mkdir
parallel --bar "go run honda2014.go --seed {1} --prior {2}" :::: seed.txt :::: prior.txt

リグレットの変化のグラフを描画するためのPythonコードは以下のとおりです。シミュレーション結果のファイルサイズが大きいので,一度Pickleに書き出しました。

import seaborn
import scipy
import matplotlib.pyplot as plt
from tqdm import tqdm
from sklearn.externals import joblib

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

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

prefixes = ['uniform', 'reference', 'jeffreys']
for prefix in prefixes:    
    A = load_grid(f'{prefix}/{prefix}-seed', 1, 20000)
    joblib.dump(A, f'{prefix}.pkl')
    del A

plt.ylim(ymax=40)
plt.xscale('log')
plt.xlim(1, 10000)

fnames = ['uniform.pkl', 'reference.pkl', 'jeffreys.pkl']
labels = [r'$\alpha=-1/2$ (uniform)', r'$\alpha=0$ (reference)', r'$\alpha=1/2$ (Jeffreys)']
for fname, label in zip(fnames, labels):
    A = joblib.load(fname)
    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)
coef = 1/(1/2*scipy.log(1+1/(0.3**2)))
plt.plot(xs, coef*ys, color='k', label='asymptotic bound')
plt.xlabel('plays')
plt.ylabel('regret')
plt.legend(loc='upper left')
plt.savefig('honda2014.png', dpi=100)

結果

追試で欲しかった図は [1] p.4のFigure 1です。\alpha によって性能が違っていますが,あくまでこれは期待値です。

f:id:kujira16:20170503223128p:plain

参考文献

  • [1] Honda, J. and Takemura, A.: Optimality of Thompson Sampling for Gaussian Bandits Depends on Priors, PMLR Vol.33 (2014) [PMLR] [arXiv].
  • [2] Yang, R. and Berger, J. O.: A Catalog of Noninformative Priors (1998).
  • [3] 本多淳也,中村篤祥:バンディット問題の理論とアルゴリズム講談社(2016).

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

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

575