どのイニングに誰が先頭打者になりやすいかを確率計算で求める

MikuHatsune2011-01-22

どのイニングに誰が先頭打者になりやすいかの続き。
P_{n}:nイニング目で打者その1、その2、…がそれぞれ、そのイニングの先頭打者になる確率ベクトル
M_{n}:1イニングで打席に立ちうる打者人数の確率推移行列
として、
P_{n+1}=M_{n}P_{n}
の形にしたい。
打者9人なら、1イニング目で1番打者から始まる。
あるk+1イニング目で、1番打者が先頭打者になろうと思ったら、そのひとつ前のkイニング目で
1番打者から始まり、9人打席に立つ。
2番打者から始まり、8人打席に立つ。
3番打者から始まり、7人打席に立つ。



9番打者から始まり、1人打席に立つ。
となればよくて、それをまあいろいろやってみよう。
チーム内にはJ人いるとしよう。
p_{n}(k):あるイニングnで、k番打者が先頭打者になる確率
b(t):1イニングで打席にt人立ちうる確率。
ただしここでは、今までのプログラムで1チームの人数を越えるくらいの打者が打席に立っても、折りたたんでいる。
例えば、9人チームで1イニング最大14人まで打席に立つことを考えたら、10人目は結局打者1人しか打席に立ってない、11人なら2人、…、という解釈。
なのでここでは、アウト数3なのに打者1もしくは2人しか回らず、第2イニングで2番打者が先頭かよ〜ってなるけど、それは前の回で10人打席に立ったということ。
そうすると、たぶん
p_{n+1}(k)=\sum_{k=1}^{J}p_{n}(k)\hspace{5}b(J+1-k)
と書けそうで、以下ごり押し。
上のリンクから借用して

Nbase<- 4
Nout<- 3
people<- 9
Vnum<- 300
team<- matrix(c(1/2,rep(1/8,Nbase)),nr=people,nc=5,byrow=T)
cut<- 0.99
Ninning<- 50
A<- rep(0,Vnum)
for(v in 1:Vnum){
  A[v]<- choose((v-1),(Nout-1))*(hitp)^(v-Nout)*(1-hitp)^(Nout-1)*(1-hitp)
}
bboxesprob<- head(A,min(which(cumsum(A)>cut)))/
             sum(head(A,min(which(cumsum(A)>cut))))
# はみ出る部分を折りたたむ。
rest<- ceiling(length(bboxesprob)/people)*people-length(bboxesprob)
bbp<- c(bboxesprob,rep(0,rest))
bp<- matrix(bbp,nr=ceiling(length(bboxesprob)/people),byrow=T)
Pn<- colSums(bp)                                                                 # 折りたたみ完了
# 推移行列を作る。
h<- function(vec){
  v<- rep(0,length(vec))
  Ps<- which(vec!=0)
  num<- vec[Ps]
  v[Ps:1]<- Pn[num:(num-Ps+1)]
  if(Ps!=length(vec)){
    v[(Ps+1):length(v)]<- Pn[1:(length(v)-Ps)]
  }else{
    v<- Pn[1:length(v)]
  }
  return(v)
}
L<- apply(diag(people,people),2,h)                                               # 推移行列完成。
# イニングを増やす。
t<- c(1,rep(0,people-1))                                                         # 1イニング目
R<- matrix(0,people,Ninning)
R[,1]<- t
for(i in 2:Ninning){
  R[,i]<- L%*%R[,i-1]
}
# 特定のイニングだけ欲しかったら
library(expm)
c<- 9
(L%^%(c-1))%*%t
# 上のリンクのプログラムを走らせて、Zを作ったうえで
par(mfrow=c(2,2))
persp(Z,phi=20,theta=60,col=5,
      xlab="Batter.No",ylab="Inning.No",zlab="First batter probability",
      main="Simulation")
persp(R,phi=20,theta=60,col=2,
      xlab="Batter.No",ylab="Inning.No",zlab="First batter probability",
      main="Calculation")
persp(Z[,1:10],phi=20,theta=60,col=5,
      xlab="Batter.No",ylab="Inning.No",zlab="First batter probability",
      main="Simulation 10th inning")
persp(R[,1:10],phi=20,theta=60,col=2,
      xlab="Batter.No",ylab="Inning.No",zlab="First batter probability",
      main="Calculation 10th inning")
# ZとRがどれくらい似てそうか考える。
ZR<- Z[,2:ncol(Z)]/R[,2:ncol(R)]


ZとRを描いたら、まあいけてそう。
ZとRの比を取ったら、まあいけてそう。
イニングが増えたとき、Rは平均値(チーム人数の逆数)に近付いていて、まあいけてそう。
まあいけてそう。
ここで、適合度検定なるもの(前回のχ自乗検定だが、言葉の使い方がまだ怪しい)をしよう。
またまた上のリンクから拝借して

Rr<- ceiling(R*10000000)      # たぶん一千万回のシミュレーションとみていいのだろう。
p<- rep(1/people,people)
p.value<- rep(0,Ninning)
for(r in 1:Ninning){
  p.value[r]<- chisq.test(Rr[,r],p=p)[3]  
}                                          
p.value<- unlist(p.value)
p.value
plot(p.value,type="o",frame=F,ylim=c(0,1),xaxt="n",
     xlab="Ninning.No",ylab="p-value",main="p-value plot")
axis(1:2)

前回のシミュレーションに比べてかなりきれいな図になった。
8イニング目で 0.011だった。