どのイニングに誰が先頭打者になりやすいかの続き。
:nイニング目で打者その1、その2、…がそれぞれ、そのイニングの先頭打者になる確率ベクトル
:1イニングで打席に立ちうる打者人数の確率推移行列
として、
の形にしたい。
打者9人なら、1イニング目で1番打者から始まる。
あるk+1イニング目で、1番打者が先頭打者になろうと思ったら、そのひとつ前のkイニング目で
1番打者から始まり、9人打席に立つ。
2番打者から始まり、8人打席に立つ。
3番打者から始まり、7人打席に立つ。
・
・
・
9番打者から始まり、1人打席に立つ。
となればよくて、それをまあいろいろやってみよう。
チーム内にはJ人いるとしよう。
:あるイニングnで、k番打者が先頭打者になる確率
:1イニングで打席にt人立ちうる確率。
ただしここでは、今までのプログラムで1チームの人数を越えるくらいの打者が打席に立っても、折りたたんでいる。
例えば、9人チームで1イニング最大14人まで打席に立つことを考えたら、10人目は結局打者1人しか打席に立ってない、11人なら2人、…、という解釈。
なのでここでは、アウト数3なのに打者1もしくは2人しか回らず、第2イニングで2番打者が先頭かよ〜ってなるけど、それは前の回で10人打席に立ったということ。
そうすると、たぶん
と書けそうで、以下ごり押し。
上のリンクから借用して
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だった。