2013-02-23

奇病の治療法の確立、正単体空間の比較・対応付け、カテゴリカルな決断

奇病が発生した。この奇病にかかると、音楽と美術とスポーツに関する能力が破壊されるという。音楽は演奏できなくなるし鑑賞もできなくなる。美術は創作活動ができなくなるし鑑賞もできなくなる。スポーツは自らプレイすることができなくなるし、鑑賞もできなくなる。
今、３つの治療法、M,A,Sが開発されたと言うが、効果は不明だという。この３つの治療法は一つしか選べない(併用すると効果は無になることが分かっている)とする。
また、治療は即刻開始しなければ、無効であり、副作用はないことにする。
M,A,Sの効果はそれぞれ、音楽系、美術系・スポーツ系にのみ効果があり、他の系には全く無効である
また、効果の出方は４通りある。
- 「著効」
- 「半効」
- 「無効」
それぞれの効果をM1,M2,M3,A1,A2,A3,S1,S2,S3とする
効果は立場によって優劣が異なる。例を挙げる。
- ホタル界ではA1>A2=S1>S2>M1=M2=M3=A3=S3と評価された
- ウグイス界ではM1=S1>S2>M2=A1>A2>M3=A3=S3と評価された
- ちなみに、選択肢２つで、その帰結が「成功と失敗」の２つの場合には、X1=Y1>X2=Y2という優劣情報になっている。

# オプション数
k <- 3
# 各オプションの帰結の種類数
s <- sample(2:5,k,replace=TRUE)
# 各オプション３個にしてみる
s <- rep(3,k)
s.list <- list()
for(i in 1:k){
	s.list[[i]] <- 1:s[i]
}
# 各オプションごとに帰結の生起確率
true.prob <- list()
for(i in 1:k){
	true.prob[[i]] <- rdirichlet(1,rep(10,s[i]))
}
# ちょっと固定してみる
# オプション２とオプション３とが競り合うように
#true.prob[[1]] <- c(0.5,0.5)
#true.prob[[2]] <- c(0.45,0.55)
true.prob[[3]] <- c(0.5,0.1,0.4)
true.prob[[2]] <- c(0.0,0.5,0.5)
# 起こりうる帰結の組合せの網羅
cat.combin <- expand.grid(s.list)

# どの帰結を好むかの行列を作る
#preference <- sample(1:k,prod(s),replace=TRUE)
preference <- matrix(1,prod(s),k)
for(i in 1:prod(s)){
	tmp.s <- sample(0:(k-1),1)
	preference[i,sample(1:k,tmp.s)] <- 0
	preference[i,] <- preference[i,]/sum(preference[i,])
}
#preference <- matrix(c(1,1,0,1,1,0,1,1),ncol=2)
# ホタルのモデル
scores <- matrix(c(1,4,3,1,3,2,1,1,1),byrow=TRUE,3,3)
for(i in 1:length(preference[,1])){
	tmp <- c(scores[cat.combin[i,1],1],scores[cat.combin[i,2],2],scores[cat.combin[i,3],3])
	preference[i,] <- as.numeric(tmp==max(tmp))
}
# 和が１となるように標準化
for(i in 1:length(preference[,1])){
	preference[i,] <- preference[i,]/sum(preference[i,])
}
# 生起確率は全く不明なところからスタート
pre <- list()
for(i in 1:k){
	pre[[i]] <- sample(0:0,s[i],replace=TRUE)
}
# そのつど、確率的によさげなオプションを選んでは、その結果をhistoryに記録し
# 記録に応じて、次の選択を判断する
n.history <- 100
history <- matrix(0,n.history,sum(s))
# 結果からディリクレ乱数で「生起確率」を推定する
# 作成するディリクレ乱数の個数
n.iter <- 1000
#N.iter <- 1

for(h in 1:n.history){
	tmp <- c()
	for(j in 1:k){
			tmp <- c(tmp,pre[[j]])
	}
	history[h,] <- tmp
	chosen <- rep(0,k)
	#for(ii in 1:N.iter){
		this.time.chosen <- rep(0,k)
		rs <- list()
		for(j in 1:k){
			rs[[j]] <- rdirichlet(n.iter,pre[[j]]+1)
		}

		# 推定した生起確率とディリクレ乱数から
		# もっとも好みに沿った結果が得られるオプションを
		# 確率的に選ぶ
		for(i in 1:n.iter){
			tmp.rs <- list()
			for(j in 1:k){
				tmp.rs[[j]] <- rs[[j]][i,]
			}
			probs <- expand.grid(tmp.rs)
			tmp.prob <- apply(probs,1,prod)
			tmp2 <- preference * tmp.prob
			tmp3 <- apply(tmp2,2,sum)
			tmp.max <- which(tmp3 == max(tmp3))
			if(length(tmp.max)==1){
				choice <- tmp.max
			}else{
				choice <- sample(tmp.max,1)
			}
			#choice <- tmp.max[sample(seq(tmp.max),1)]
			this.time.chosen[choice] <- this.time.chosen[choice]+1
			#print(choice)
			#tmp <- matrix(NA,N.iter,k)
			#for(j in 1:k){
			#	tmp[,j] <- sample(1:s[j],replace=TRUE,prob=rs[[j]])
			#}
		}
		#tmp.max <- which(this.time.chosen==max(this.time.chosen))
		#if(length(tmp.max)==1){
		#	selected.chosen <- tmp.max
		#}else{
		#	selected.chosen <- sample(tmp.max,1)
		#}
		#selected.chosen <- tmp.max[sample(seq(tmp.max),1)]
		chosen[selected.chosen] <- chosen[selected.chosen]+1
	#}
	#tmp.max <- which(chosen==max(chosen))
	#if(length(tmp.max)==1){
	#	really.chosen <- tmp.max
	#}else{
	#	really.chosen <- sample(tmp.max,1)
	#}
	# 選んだオプションの結果は、隠された真の生起確率によって観察される
	really.chosen <- sample(1:k,1,prob=this.time.chosen/sum(this.time.chosen))
	#really.chosen <- tmp.max[sample(seq(tmp.max),1)]
	really.happen <- sample(1:s[really.chosen],1,prob=true.prob[[really.chosen]])
	pre[[really.chosen]][really.happen] <- pre[[really.chosen]][really.happen]+1

}

matplot(history,type="l")