2011-01-18

子供と遊びを始めるとき

R ゲーム

「(簡単な)ゲームをしよう！」
「んーと、じゃあ、○○ちゃんが３(k)回勝つまでね」
よくある光景である
さて、ゲームを何回やることになるのだろうか
○○ちゃんの勝率 $p$ とする
$N$ 回目に「おわりー」となるとする
$1,2,...,N$ 回までに「おわりー」とている確率は、 $N$ 回終わった段階で、○○ちゃんが、まだ $1,2,...,k-1$ 回しか勝っていない場合を除いたすべての場合の確率の和だから
$Q(N)=1-\sum_{i=0}^{k-1} _N C_i p^i(1-p)^{N-i}$ である
$N$ 回目に「おわりー」となるのは、 $N-1$ 回目で「おわりー」とならずに、 $N$ 回目で「おわりー」となった場合なので $Q(N-1)-Q(N)$ である

maxN<-20
n<-4
p<-0.3
# N回目が終了してもまだ「おわり」になっていない確率
ret<-rep(0,maxN)

for(i in 1:maxN){
	tmp<-0
	for(j in 0:(n-1)){
		tmp<-tmp+choose(i,j)*(1-p)^(i-j)*p^(j)
	}
	ret[i]<-1-tmp
}

# N回目で「おわり」になる確率は
P<-diff(ret)

par(mfcol=c(1,2))
plot(P,type="b",main="打ち切り回数ごとの確率")
plot(cumsum(P),type="b",main="打ち切り回数の累積確率")
par(mfcol=c(1,1))

少し複雑にしよう
C枚のカードがある( $c_1,c_2,...,c_C$ )
C個の袋がある( $b_1,b_2,...,b_C$ )
今、必ず $c_1$ のカードから始めるのだが、次のカードをどれにするかは、 $b_1$ 番目の袋に入れてある札(0,1,2,...,Cの数字が書いてある)を引いて、札に書いてある数字のカードに持ち替えることにする。ただし、0の札を引いたら、カードは持ち替えない。
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回までの確率を求め続けることは、ただの行列計算

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)をいじりつつ進めればよい