前回
問題設定
行ったのはHondaらのAISTATS 2014の論文 [1] の中の実験の追試で,報酬が , の正規分布に従うアームと , の正規分布に従うアームがそれぞれ1本ずつあるという設定の多腕バンディット問題です。
原理
最初にパラメータ を -1/2, 0, 1/2 あたりの値に設定します。 については後述します。
Thompson Samplingを始める前に,まずそれぞれのアームを 回引いておきます。このとき得られた報酬も個々のアームについて記録しておきます。
この後からは,個々のアーム についてStudentの 分布 から乱数を発生させ,最も大きい乱数が得られたアームを選択します。ここで はアーム を引いた回数, はアーム から得られた報酬の標本平均, はアーム から得られた報酬の標本標準偏差(Nで割るほう)です。
Studentの 分布はどこから出てくるのかというと,データ が得られたもとで平均 と標準偏差 が未知の正規分布のパラメータの分布 をベイズの定理から計算し,これを について周辺化すると出てくるようです [2]。このとき,事前分布 の設定方法には諸説あるらしく,一様分布を設定したものが に,reference priorを設定したものが に,Jeffreys priorを設定したものが にそれぞれ対応するようです。
前回の実験では事前分布は結果にほとんど影響を及ぼさなかったのですが,今回の実験では違っていて,設定した事前分布によって結果が大きく異なっています。このあたりの難しい話については追えていないので [1] を参照してください。
ソースコード
以下のプログラムはコマンドラインでシード値と事前分布を与えるとThompson Samplingに従ってアームを10^4回引くシミュレーションを行うものです。報酬の標準偏差を計算するときには,得られた報酬の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です。 によって性能が違っていますが,あくまでこれは期待値です。
参考文献
- [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).
バンディット問題の理論とアルゴリズム (機械学習プロフェッショナルシリーズ)
- 作者: 本多淳也,中村篤祥
- 出版社/メーカー: 講談社
- 発売日: 2016/08/25
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (1件) を見る
575
報酬が
— しょラー (@shora_kujira16) 2017年5月3日
分散未知の
ガウシアンhttps://t.co/2pNjdjOBMj