2013-04-13

１番になるための戦略

２つの選択肢があって、帰結が0/1であるときの決断方法について、ここなどで、選択肢が連続値で帰結も連続値であるときのことをここで書いた(そこではある変数によって定まる関数に関して積分する話があって、それは、functional integral 汎関数積分と呼ぶらしい。基本的にはそれは解析的には解けないらしい)
今日は、選択肢は連続値で帰結は0/1である場合にする(問題を簡単にしよう)
選択肢はx=[0,1]と有限範囲を仮定し、そのうえで、xの値に対して成功率を与える関数があるとする。たとえばこんな感じ。なんらかの規則はあるが、きれいな関数ではない、というようなそんなもの

true.p <- function(x){
	ret <- x^2*0.6+0.2
	s <- which(x<0.3)
	ret[s] <- -(1-x[s])^2*(0.4-x[s])^3*(x[s]-0.8)*0.8+0.9
	#ret <- x
	ret <- (sin(x*4*pi)+1)*0.3 +(x-0.4)^2
	#ret <- (ret-min(ret))/(max(ret)-min(ret))
	ss <- sample(1:length(x),length(x)*0.5)
	ret[ss] <- sample(ret[ss])
	#ret <- runif(length(x))
	ret
}

x <- seq(from=0,to=1,length=100)
p <- true.p(x)
plot(x,p)

背後にこんな成功率が潜んでいるけれど、まったくなんの情報もないところからスタートすることにする
ある値xに関して、繰り返して帰結を観察すれば、その値xにおける成功率はベータ分布で考えても悪くないだろう
問題は、ある値xに繰り返して観察すると他の値xに関する情報が得られないこと
その点を解決するために、xの範囲を適当にビンに区切って、そのビンごとに成功率が一定だという仮定をした上で、ベータ分布を想定してみることにする
それを進める上で問題になるのは、ビンの幅をどうすればよいのか、ということ
１次元の量的データの場合に分布を図示方法としてヒストグラムがある。ヒストグラムの場合には、ビンの幅をどうするかの方法があってRではhist()関数の中にビン幅決めのデフォルト方法も実装してある
同様にdensity()関数にもある
これらのビン幅決めルールを流用することにする
基本的には、このビン幅決めルールは、xの範囲に関して均等幅のビンを取る
そのビン幅決め戦略は、ビン幅がだんだん狭くなっていくのだが、ビンが狭くなってそのビンにまったく帰結情報がない状態になると、そのビンの成功確率の分布は[0,1]で均等分布になることに留意する。情報のあるビンはどこも成功しにくいようなそんな真の成功率分布があるときには、このような「情報なしのビン」は相対的に選ばれやすくなり、帰結情報のないビンはなくなりがちになる。逆に、「情報ありのビン」の真の成功率が、高めのときは、情報のないビンは選ばれにくくなるが、情報ありのビンのxが繰り返し選ばれるときも、hist(),density()関数的には、ビンをどんどん狭くするわけではないから、情報なしビンが増えることは防止される

# density()関数を使ってビン幅を決める
# 標本数が1のときはビン幅を大きくしてある
choose.bw <- function(x){
	if(length(x)==1){
		bw <- 2
	}else{
		d <- density(x)
		bw <- d$bw
	}
	bw
}

# モンテカルロで算出する
# n.trialオプションは、何回のモンテカルロ試行でそれを決めるか
# silentオプションはちょっと興味のあるプロットをしながらの処理にするかどうかの指定

# x is selected option value
# y is outcome value for x
# z is option value vector
choice.prob <- function(x,y,z=seq(from=0,to=1,length=100),n.trial=1000,silent=TRUE){
	if(length(x)==0){
		return(rep(1/length(z),length(z)))
	}
	bw <- choose.bw(x)
	if(!silent){
		if(length(which(y==1))>2){
			#plot(hist(x[which(y==1)]))
			plot(density(cbind(x,y)))
		}
		
	}
	m <- matrix(0,length(z),2)
	for(i in 1:length(z)){
		st <- z[i]-bw/2
		end <- z[i]+bw/2
		selected <- which(x>=st & x< end)
		m[i,1] <- length(which(y[selected]==0))
		m[i,2] <- length(which(y[selected]==1))
	}
	R <- matrix(0,length(z),n.trial)

	for(i in 1:length(z)){
		R[i,] <- rbeta(n.trial,m[i,2]+1,m[i,1]+1)
	}
	ranks <- apply(R,2,rank,ties.method="max")
	top.only <- ranks
	top.only[which(top.only != length(z))] <- 0
	top.only[which(top.only == length(z))] <- 1
	s.top.only <- apply (top.only,2,sum)
	s2.top.only <- t(t(top.only)/s.top.only)
	s3.top.only <- apply(s2.top.only,1,sum)
	s3.top.only/sum(s3.top.only)
}

# 順番にn.iter回、試してみる
n.iter <- 1000
# 選択されたxを第１列に、その帰結を第２列に格納する
x.history <- matrix(0,n.iter,2)
#z <- seq(from=0,to=1,length=100)
# x.historyの第１行は、初期選択(ランダムに選ぼう)
tmp.x <- runif(1)
tmp.p <- true.p(tmp.x)
x.history[1,1] <- tmp.x
x.history[1,2] <- sample(1:0,1,prob=c(tmp.p,1-tmp.p))
# ２回目以降は、それまでの情報を活用する
for(i in 2:n.iter){
	# 選択肢は網羅的にやってもいいけれど、有限にしておく
	z <- runif(100)
	# これまでのデータに従って、選択肢ごとにそれが最高の成功率を示すかの相対的重みを返す
	# モンテカルロで算出する
	# n.trialオプションは、何回のモンテカルロ試行でそれを決めるか
	# silentオプションはちょっと興味のあるプロットをしながらの処理にするかどうかの指定
	tmp.prob <- choice.prob(x.history[1:(i-1),1],x.history[1:(i-1),2],z,n.trial=100,silent=TRUE)
	# オプションごとの相対的重みづけをプロット
	plot(z,tmp.prob,type="h")
	# 重みに応じて、選択する
	tmp.x <- sample(z,1,prob=tmp.prob)
	# その選択肢によって帰結が確率的に決まる
	# それを記録する
	x.history[i,1] <- tmp.x
	tmp.p <- true.p(tmp.x)
	x.history[i,2] <- sample(1:0,1,prob=c(tmp.p,1-tmp.p))
}

par(mfcol=c(2,2))
# 横軸が仮説。縦軸が帰結
plot(x.history)
# 真の成功率のプロット
plot(x,p)
# 選ばれた仮説の履歴
plot(x.history[,1])
# 帰結の履歴
#plot(x.history[,2])
# 成功の累積
# 赤のラインは、最高の成功率のときの予測ライン
# 緑のラインは、最低の成功率のときの予測ライン
# 青のラインは、平均の成功率のときの予測ライン
# 赤のラインの傾きに収束してくれれば、戦略としては成功
plot(cumsum(x.history[,2]))
abline(0,max(p),col=2)
abline(0,min(p),col=3)
abline(0,mean(p),col=4)
# 成功率の履歴
plot(true.p(x.history[,1]))