2016-08-25

ベルマン方程式

ベルマン方程式動的計画法マルチアームド・バンディット Gittins index

調べものをしましたが、日本語サイトの中でこちらが一番、わかりやすく参考にさせていただきました
ベルマン方程式というものがある。Wikipediaではベルマン方程式(リチャード・E・ベルマン)、もしくはBellman_equation
選択肢があるときに、どれを選ぶとよいかを迷いながら選ぶような状況があるとする
どれを選ぶかの基準は、「選ぶことで得られる利益の総和」を大きくしたい、という基準
この基準に照らして、もし最適な選択肢を取ることができるのであれば、その利益の総和はどうなるだろうか、というようなことを表した方程式
選択肢を選んだときに、どのくらいの利益が生じるかは確率的に決まる（ある時は利益が出るし、あるときは出ない、とか）ので、常に、どの選択肢を選んだものか…と悩ましい状態に置かれるが、それに関する式だということ
この論文では、複数の選択肢があって、選ぶたびに利益は0か1かになるような単純な選択課題であるが、そこに出てくる式(3.1)を理解したい
式は以下のような形をしている
- $V(x_{1,t},x_{2,t},...,x_{K,t})=$
- $max_{k} \{(profit_{k,t} +$
- $d(q_{k,t} \times V(x_{1,t},...,x_{k,t}+(1,0),...) + (1-q_{k,t}) \times V(x_{i,t},...,x_{k,t}+(0,1),...)\}$
各選択肢について、0回以上の試行結果の成否情報がある。それが $x_{i,t}$ 。成否の回数からなる長さ２の整数ベクトル。ただし、iは第i番目の選択肢。tは、これまでの全部でt試行がなされていることを意味する（複数の選択肢があるので、それぞれが実行された回数の総和になる）。
成否情報があったときに、最適選択をするとすると、その先で得られる利益の総和（の期待値）が $V(x_{1,t},...)$ 。全試行数がtのときのそれ
このような状態( $x_{i,t}$ )にあって、最善選択肢が第k番だとなっている、とする。
そうすると、全試行数がt+1に増えた段階では、第k番選択肢を選ぶことで得られる利益の期待値 $proit_{k,t}$ が利益となることはわかる。ただし、この追加利益は未知であるから、tの時点でのその推定値を使うことになる。二値成否の場合には、ベータ分布から知られる期待値を使うのがよいので $profit_{k,t} = \frac{s_{k,t}+1}{s_{k,t}+f_{k,t}+2}$ 、というように成否回数 $s,f$ と下駄に当たる、１、２を使った期待値を使っているのが、リンク先の論文
この第k番選択肢が選ばれることがわかっているものとして、t+1以降にどんな利益が続くのかは不明だが、それを $V$ で表してしまっている。ただし、このt+1以降で考慮するべきなのは、tにて第k番選択肢を選んで、成功する確率qと失敗する確率(1-q)とがあって、それぞれの場合について、起きる確率の重みを考慮して加える必要がある
また、このqについては、ベータ分布にて成功確率を推定するとすれば、 $q= \frac{s_{k,t}+1}{s_{k,t}+f_{k,t}+2}$
さて、こんな風なVに関する確率コミの漸化式ができたので、これを用いて、推定することを考える
全体で最適選択をしたら、どんな総利益が出るかを一発で出すことはできないけれど、漸化式をちまちま解いていくと最終的には答えがでるね、このように全体の答えを解くのに、ちまちまと部分部分で考えるのが動的計画法(こちら)
ちなみにdは、繰り返して実行するにしても、そのうち止めるという現実的な条件を入れるためのたで、継続する確率をdとしている。無限に続けるならd=1
こちらにベルヌーイ事象のGittins indexのGlazebrook Weberアルゴリズム（ベタなアルゴリズム）のMatlabコードがある。それのベタアルゴリズムの方のR化したもの（バグが潜んでいる可能性あり）。一応、答え合わせは合っている感じ…
ぐるぐるとループが３重に回っているので遅いです
0-1の値についてグリッド探索すると得られることが知られているらしい。~~ペイパーのAppendixの書いてあるのだが…、以下のソースがそれを実装しているのかどうか、ちょっと心もとないが、結果があっているのでよしとする~~
成功率は不明であるが、それがわかりさえすれば、繰り返して実施したときの総報酬期待値はその成功率に関して、1/(1-d)倍したものになる
なので、その成功率をグリッドで定め多ものに基づく総報酬期待値について、小さい方から、試してみる。で、試した値が、漸化式の値に近いことが解れば、その総報酬期待値を「当該観測成否回数」における総報酬期待値とする
漸化式の値に近いことは、「成功率グリッド値を小さい方から試して行って、初めて、漸化式の値を超えたあたり」という判定方法とする
そのようにしたときの成功率グリッド値(グリッドの幅を考慮して少し加減した値)をGittins indexとして登録する

my.gittins <- function(d,Mh,grid){
	G <- matrix(0,Mh-1,Mh-1)
	Gittins <- matrix(0,Mh-1,Mh-1)
	
	for(s in 1:(Mh-1)){
		G[s,Mh-s] <- s/Mh
	}
	ps <- seq(from = grid/2, to =1, by=grid)
	for(ii in 1:length(ps)){
		p <- ps[ii]
		# 1/(1-d) = 1 + d + d^2 + .....
		# p/(1-d)はdの確率で継続するときの…、全報酬の期待値
		known <- p/(1-d)
		for(i in (Mh-1):2){
			for(s in 1:(i-1)){
				experimental <- s/i * (1+d*G[s+1,i-s]) + (i-s)/i * d * G[s,i-s+1]
				if(Gittins[s,i-s] ==0 & (known > experimental)){
					Gittins[s,i-s] <- p - grid/2
				}
				G[s,i-s] <- max(known,experimental)
			}
		}
	}
	

	
	
	return(list(G=G,Gitting=Gittins))
}

Git09 <- my.gittins(0.9,200,0.001)
Git09[[2]][1:8,1:8]

全試行回数を有限に設定し、そこから遡る形で計算
結果

> Git09[[2]][1:8,1:8]
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,] 0.703 0.500 0.380 0.302 0.249 0.210 0.181 0.159
[2,] 0.800 0.635 0.516 0.434 0.372 0.324 0.287 0.257
[3,] 0.845 0.707 0.601 0.518 0.456 0.406 0.365 0.331
[4,] 0.872 0.754 0.658 0.581 0.518 0.468 0.426 0.390
[5,] 0.891 0.787 0.700 0.628 0.568 0.517 0.475 0.439
[6,] 0.904 0.811 0.732 0.664 0.607 0.558 0.516 0.479
[7,] 0.914 0.831 0.757 0.694 0.640 0.592 0.551 0.514
[8,] 0.922 0.846 0.778 0.719 0.667 0.621 0.581 0.545

Gittins indexとして論文に掲載されているもの

g1 <- c(0.7029,0.8001,0.8453,0.8723,0.8905,0.9039,0.9141,0.9221)
g2 <- c(0.5001,0.6346,0.7072,0.7539,0.7869,0.8115,0.8307,0.8461)
g3 <- c(0.3796,0.5163,0.6010,0.6579,0.6996,0.7318,0.7573,0.7782)
g4 <- c(0.3021,0.4342,0.5184,0.5809,0.6276,0.6642,0.6940,0.7187)
g5 <- c(0.2488,0.3702,0.4561,0.5179,0.5676,0.6071,0.6395,0.6666)
g6 <- c(0.2103,0.3245,0.4058,0.4677,0.5168,0.5581,0.5923,0.6212)
g7 <- c(0.1815,0.2871,0.3647,0.4257,0.4748,0.5156,0.5510,0.5811)
g8 <- c(0.1591,0.2569,0.3308,0.3900,0.4387,0.4795,0.5144,0.5454)

G.mat <- rbind(g1,g2,g3,g4,g5,g6,g7,g8)

> t(G.mat)
         g1     g2     g3     g4     g5     g6     g7     g8
[1,] 0.7029 0.5001 0.3796 0.3021 0.2488 0.2103 0.1815 0.1591
[2,] 0.8001 0.6346 0.5163 0.4342 0.3702 0.3245 0.2871 0.2569
[3,] 0.8453 0.7072 0.6010 0.5184 0.4561 0.4058 0.3647 0.3308
[4,] 0.8723 0.7539 0.6579 0.5809 0.5179 0.4677 0.4257 0.3900
[5,] 0.8905 0.7869 0.6996 0.6276 0.5676 0.5168 0.4748 0.4387
[6,] 0.9039 0.8115 0.7318 0.6642 0.6071 0.5581 0.5156 0.4795
[7,] 0.9141 0.8307 0.7573 0.6940 0.6395 0.5923 0.5510 0.5144
[8,] 0.9221 0.8461 0.7782 0.7187 0.6666 0.6212 0.5811 0.5454