遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

打ち切り型のゲーム

R ゲーム
  • 前の記事では何回で終了するかわからない打ち切りゲームについて考えた
  • 少し複雑にしよう
  • こちらで交叉組み換えの計算を行列に乗せていた
  • 行列に乗せると個々の計算を忘れられるのがよいのと、行列計算はコンピュータが(Rが)、軽い動作で実行してくれるから
  • そのメリットを享受するためには、少し我慢をして、複雑な場合わけをすることにする(推移行列を大きくする)ことにする
  • 題材は野球の得点分布(こちら)
  • 野球の得点分布を出すのも、「いつになったらスリーアウトになるかわからない」ので、打ち切り型のゲームである
  • 野球よりもう少し単純なモデルで考えた上で、野球の１イニングの得点分布を作るソースを書こう
  • C枚のカードがある(c_1,c_2,...,c_C)
  • C個の袋がある(b_1,b_2,...,b_C)
  • 今、必ずc_1のカードから始めるのだが、次のカードをどれにするかは、c_iカードを持っているときにはb_i番目の袋に入れてある札(0,1,2,...,Cの数字が書いてある)を引いて、札に書いてある数字のカードに持ち替えることにする。ただし、0の札を引いたら、カードは持ち替えない。ただし、b_i番目の袋に入っている札の枚数はそれぞれ違うものとする
  • 0のカードをn回引いたら、終了することにする。
  • C枚のカードには、それぞれ、減点数p_1,...,p_Cが決まっているので、終了時点では、終了までにN回、札を引いたとしたら、N-p_1を、このときの得点とすることにしている
  • さて、何点とる確率がどれくらいになるだろうか
  • i回の札を引いた時点で、C枚のうちのどのカードを持っているかは、推移行列Mを用いてM^i Sで表される。ただしSは、C枚のカードのどれを持っているかの確率ベクトル
  • また、この時点で0を引いた回数は0,1,2,...,n回のうちのいずれか。
  • 子供との「おわりー」と同じように考えよう
  • i回引いた時点で、0を引いた回数が0,1,2,...,n回のそれぞれの場合にベクトルSが求まる
  • これをすべて格納しつつ、どんどんと回数を増やしていき、0引き回数がn回までの確率を求め続けることは、ただの行列計算
  • さて、野球。
  • 掲載図は３２ｘ３２推移行列
    • i行j列成分は状態jから状態iへの推移確率
    • アウトが３より大きくなる場合も「起こり得る」ものとし、その上で、そのような状態については計算しないことにした「ちぎれた推移行列」であることに注意
    • アウトにならない場合(1,2,3,4塁打)はアウト数が同じ領域(８ｘ８スペース)での変化
    • アウトになる場合(0)は、アウト数を１つ増やしながら、ベースの埋まり方は変えない変化
    • アウト数が減ったり、２個以上増えたりすることはない様子もわかる

http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2010/baseballMatrix.png

library(MCMCpack)
# １チームの打者数
Nmember<-9
# 各打者の推移行列を納めるリスト
Mlist<-list()
# 推移行列をNmember個作る
for(s in 1:Nmember){
	# 各打者の打率はランダムに決める
	# b[1]はシングルヒット、b[5]はアウトの確率
	b<-c(rdirichlet(1,rep(1,5)))
	# ベースの埋まり方、(000),(100):1塁だけ埋まっている,(010),(001),(110),(101),(011),(111)の８通りに関して　８ｘ８推移行列を作る
	m<-matrix(c(rep(b[4],8),c(b[1],0,0,b[1],0,0,0,0),c(b[2],0,b[2],b[2],0,0,b[2],0),rep(b[3],8),
c(0,b[1],0,0,0,b[1],0,0),c(0,0,b[1],0,0,0,b[1],0),c(0,b[2],0,0,b[2],b[2],0,b[2]),c(0,0,0,0,b[1],0,0,b[1])),8,8,byrow=TRUE)
# ベースの埋まり方８通りとアウトカウントの場合の数0,1,2,3の両方を合わせた場合の数は８ｘ４＝３２なので、32x32推移行列を作る
# アウトになる場合とアウトにならない場合とで、推移行列は先に作った
# 8x8行列を当てはめて作れる
	M<-matrix(0,4*8,4*8)
	for(i in 1:8){
		for(j in 1:8){
			M[i,j]<-M[i+8,j+8]<-M[i+16,j+16]<-M[i+24,j+24]<-m[i,j]
		}
	}
	for(i in 1:(3*8)){
		M[i+8,i]<-b[5]
	}
	Mlist[[s]]<-M
}


# 初期状態は０アウト、ランナーなしに確率が１
v<-c(1,rep(0,7+8*3))
library(expm)
# 最大打者数
TrialMax<-100
# 第i番目の打者のときの状態(アウト数ｘベースの埋まり方の３２の場合ごとの確率を納める行列
vs<-matrix(0,TrialMax,length(v))
# 初期状態を格納
vs[1,]<-v
# 打順をまわすためのカウンター
counter<-1
for(i in 2:TrialMax){
	#vs[i,]<-(M%^%(i-1))%*%v
	# 行列掛け算のみを続けると
	# スリーアウトになった後、次の打者がアウトにならない場合の分が多くなるので
	# 次打席の結果の計算では、「スリーアウトの状態の確率を０」にして進める
	effectivevs<-vs[i-1,]
	effectivevs[25:32]<-0
	vs[i,]<-Mlist[[counter]]%*%effectivevs[]
	# 次の打者、一巡したら１に戻す
	counter<-counter+1
	if(counter==Nmember){
		counter<-1
	}
}
# 32状態の確率の推移
matplot(vs,type="l")
# 得点計算には、「スリーアウト」の８ベース状態のみが必要なので取り出す
ZanruiProb<-vs[,25:32]
# 8状態の残塁人数
ZanruiPenalty<-c(0,1,1,1,2,2,2,3)
# 得点は、送り出した打者数から残塁人数を引いたものなので
# 送り出した打者数をベース状態のそれぞれに当てはめた行列を作る
PointMatrix<-matrix(rep(1:TrialMax,8),ncol=8)
# 残塁分を引く。これが得点
PointMatrix<-t(t(PointMatrix)-ZanruiPenalty)

# 得点を納めた行列と確率の行列ができたので
# 得点別に確率を足し合わせる
Points<-c(PointMatrix)
Probs<-c(ZanruiProb)
maxPoint<-max(Points)
ProbPerPoints<-0:maxPoint

for(i in 1:length(ProbPerPoints)){
	tmp<-which(Points==i)
	ProbPerPoints[i]<-sum(Probs[tmp])
}
# これが１イニングの得点分布
plot(0:maxPoint,ProbPerPoints,type="b")

# 複数のイニングをするときは、先頭打者(counter)をいじりつつ進めればよい