廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

かんたんなバンディットアルゴリズムのシミュレーション

今日の川柳

報酬がベルヌーイ分布に従うときのトンプソンサンプリングをシミュレーションしてみました。

問題設定

スロットマシーンが5台あってそれぞれ「当たり」が出る確率が違うとする。

「当たり」の出やすいスロットマシーンをうまく選んで、得られる「当たり」の回数を最大化したい。

トンプソンサンプリング

https://web.stanford.edu/~bvr/pubs/TS_Tutorial.pdf

によればトンプソンサンプリングのアルゴリズムは以下で与えられる。

f:id:abrahamcow:20180713214858p:plain

r が報酬(あたり・はずれの二値)x が選択肢です。

ふーん、ベータ分布を事前分布として設定して、ベイズ更新しながら当たりの出やすいやつを選ぶのね。

これをみて思うのは「乱数使う必要ある?」ってことじゃないだろうか。単純に期待値の高いやつを選んでいっても同じことになるのでは?

はい、そういうアルゴリズムもあってそれはグリーディー(貪欲法)とかよばれるらしい。

貪欲法のアルゴリズムは以下です。

f:id:abrahamcow:20180713215324p:plain

以下では R によるシミュレーションで貪欲法とトンプソンサンプリングを比べます。

シミュレーション

トンプソンサンプリングをシミュレートするコードはこうなった。

simTS <- function(seed,Tim=100){
  set.seed(seed)
  p <- c(0.3,0.4,0.5,0.6,0.7)
  a <- rep(1,5)
  b <- rep(1,5)
  r <- integer(Tim)
  for(t in 1:Tim){
    theta <-rbeta(5,a,b)
    x <-which.max(theta)
    r[t] <-rbinom(1,1,p[x])
    a[x] <- a[x] + r[t]
    b[x] <- b[x] + 1-r[t]
  }
  list(r=r,a=a,b=b)
}

貪欲法はこう。

simGreedy <- function(seed,Tim=100){
  set.seed(seed)
  p <- c(0.3,0.4,0.5,0.6,0.7)
  a <- rep(1,5)
  b <- rep(1,5)
  r <- integer(Tim)
  for(t in 1:Tim){
    theta <-a/(a+b)
    x <-which.max(theta)
    r[t] <-rbinom(1,1,p[x])
    a[x] <- a[x] + r[t]
    b[x] <- b[x] + 1-r[t]
  }
  list(r=r,a=a,b=b)
}

5台のスロットマシーンの当たりの確率は (0.3,0.4,0.5,0.6,0.7) とした。

とりあえず1000回回す。

ベータ事前分布のハイパーパラメータは (1, 1) とした。

トンプソンサンプリングにより得られた報酬の系列をr1、貪欲法により得られた報酬の系列をr0とおく。

r1、r0の30期ごとの移動平均を取ったのが下の図です。

f:id:abrahamcow:20180713215743p:plain

require(ggsomestat)
require(tidyverse)
res1 <-simTS(1,1000)
res0 <-simGreedy(1,1000)
df <-data.frame(t=1:1000,r1=res1$r,r0=res0$r)
ggplot(df,aes(x=t))+
  stat_ma(aes(y=r1,colour="r1"),windowsize = 30)+
  stat_ma(aes(y=r0,colour="r0"),windowsize = 30)+
  labs(colour="",y="reward")

(ggsomestat は自作パッケージです。ggplot2 のためのいくつかの関数をパッケージ化しました - 廿TT)

r1のほうは0.7くらいを中心に推移しており、一番当たりの出やすいスロットマシーンをうまく選んでいそう。

一方でr0のほうは0.5くらいを中心に推移している。

なんでこうなるのか。事後分布の形をみてみよう。

トンプソンサンプリングによって得られた事後分布はこちら。

f:id:abrahamcow:20180713220301p:plain

cCyan <- "#00a0e9"
cMagenta <- "#e4007f"
cGreen <- "#009944"
cOrange <- "#f39800"
cLightBlue <- "#0068b7"
qcolours <- c(cCyan,cMagenta,cGreen,cOrange,cLightBlue)

curve(dbeta(x,res1$a[1],res1$b[1]),ylim=c(0,30),xlab = "p",ylab = "density",
      col=qcolours[1],lwd=2)
for(i in 2:5){
  curve(dbeta(x,res1$a[i],res1$b[i]),add = TRUE,n=1000,
        col=qcolours[i],lwd=2)
}
legend("topright",legend = 1:5,lty=1,col=qcolours)

(色使いは 統計データをggplot2で可視化するときの色 - joker8phoenix's diary を参考にしている。)

事後分布が急峻な形になっているのは、それだけ試行回数が多いことを意味する。

0.7のスロットマシーンを一番多く選びながら、一番確率の低い0.3のスロットマシーンをかわしつつ、他のスロットマシーンの成功確率はそこそこうまく推定されている。

貪欲法によって得られた事後分布はこちら。

f:id:abrahamcow:20180713220617p:plain

curve(dbeta(x,res0$a[1],res0$b[1]),ylim=c(0,30),xlab = "p",ylab = "density",
      col=qcolours[1],lwd=2)
for(i in 2:5){
  curve(dbeta(x,res0$a[i],res0$b[i]),add = TRUE,n=1000,
        col=qcolours[i],lwd=2)
}
legend("topright",legend = 1:5,lty=1,col=qcolours)

0.5のスロットマシーンばかりを選んでいる。おそらく、最初の難解かで0.5のやつの成功が多くて、それに固執してしまって他をうまく探索できなかったのだろう。

次に、短い期間(20期)での探索を並行して100個回した。

下図は平均とその標準誤差のプロット。

f:id:abrahamcow:20180713221142p:plain

r1_100 <-sapply(1:100, function(x)simTS(x,20)$r)
r0_100 <-sapply(1:100, function(x)simGreedy(x,20)$r)

df1_100 <-r1_100 %>% 
  as.data.frame() %>% 
  mutate(t=row_number()) %>% 
  gather(key,reward,-t) %>% 
  mutate(strategy="TS")
df0_100 <-r0_100 %>% 
  as.data.frame() %>% 
  mutate(t=row_number()) %>% 
  gather(key,reward,-t) %>% 
  mutate(strategy="Greedy")
df100 <-bind_rows(df0_100,df1_100)

ggplot(df100,aes(x=t,y=reward,colour=strategy))+
  stat_summary(position =position_dodge(0.5))

平均的には貪欲法でも0.7のやつを選ぶほうに収束しそう。

でも、トンプソンサンプリングのほうが序盤の成績がいい。これはぼくには意外だった。

感想

むかし囲碁のつよい友達が「○手先まで一応読むけど、結局どこに打つかは勘で決める」みたいなこと言ってたのを思い出した。うまく偶然の力を借りれるやつが強いのかもしれない、アカギみたいに。

  • おもしろいのでもっと複雑なモデルでやってみたい。
  • でもその場合「ベイズ更新」が使えなくなるので、どうするんだろう。
  • 本業の方にいかせないだろうか。たとえば多剤併用療法の最適化とか……。無理か。
  • 広告配信とかはもう全部これでいいんじゃないか。安易か。
  • 何らかの意味での最適性って保証されてるのか気になる。

天もアカギもいいけど麻雀マンガではこれが一番好き。ノーマーク爆牌党

以下の本は読んでません。

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

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

おもしろいのかな。