どのイニングに誰が先頭打者になりやすいか

野球の続き。チームの人数を複数にしたとき、誰が先頭打者になって、そこから何人が打席に立つかで、1イニングに取る得点が変わるだろう、ということを考えた。
この日は、全イニングの平均を出したので、1番バッターが多いのは当たり前だね。
シミュレーション的に、各イニングで誰が先頭打者になっているか、1000000回試行してみよう。

#team<- rdirichlet(people,c(5,rep(1,Nbase)))
team<- matrix(c(1/2,rep(1/8,Nbase)),nr=people,nc=5,byrow=T)
hitp<- 1-colMeans(team)[1]
Vnum<- 100   
cut<- 0.99
Ninning<- 50
people<- 9
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))))
cut.off<- min(which(cumsum(A)>cut))
cut.off
Nbatter<- cut.off
countfunc<- function(vec){
  count1to.<- rep(0,people)
  for(result in 1:people){
    count1to.[result]<- length(which(vec==result))
  }
  return(count1to.)
}
N<- 10000
counter1<- matrix(0,N,Ninning)
dimnames(counter1)<- list(1:N,1:Ninning)
names(dimnames(counter1))<- c("cycles","Inning")
Bnum<- rep(1:people,Ninning*min(which(cumsum(A)>cut)))
for(I in 1:N){
  k<- sample(c(1:length(bboxesprob)),prob=bboxesprob,size=Ninning-1,replace=TRUE)
  counter1[I,]<- Bnum[c(1,cumsum(k)+1)]
}
Zz<- apply(counter1,2,countfunc)
Z<- Zz/N
Z

Zの行は打者順、列はイニング数を表わしている。要素は確率。
Zzは生データ。
ここから、先頭打者のなりやすさに違いがあるか、χ自乗検定を使ってみたい。
よくあるサイコロ問題で、問題のサイコロが精密にできているかな、というやつ。
オレみたいな暇人が、600回もサイコロを独りさみしく振ったらしい。
精密にサイコロが作られているとしたら、全部100回ずつ出る…んじゃないかな?

これらが、100回（確率100/600）に近いか、遠いか、的な。
計算はwikiの通り分散的なやり方でできて、これをχ自乗分布と有意水準との比較で、有意差があるかどうかみる。
これをするのに、chisq.testという関数があるようだ。
調べたいベクトルx
確率p
に対してχ自乗検定を行う。

x<- c(92,98,105,103,101,101)
p<- rep(100,6)
chisq.test(x,p=p)
以下にエラー chisq.test(x, p = p) :  確率の総和は 1 でなければなりません。

確率を指定しないとダ〜メ。

chisq.test(x,p=p,rescale.p=TRUE)   # p=p/sum(p) も可。

        Chi-squared test for given probabilities

data:  x 
X-squared = 1.04, df = 5, p-value = 0.9593

これらはようわからんので後で見ておく。ここで欲しいのはp値。

chisq.test(x,p=p,rescale.p=TRUE)[3]   # p=p/sum(p) も可。
$p.value
[1] 0.9592754

とすると、p値だけ取ってこれる。形式はlist。
p値はかなり1に近い。相等精密に作られているっぽいよ。
独りさみしくサイコロを延々と振った甲斐があったね。
　
さ、やり方がわかったところで
どのイニングから、等確率で先頭打者になりやすそうかな。

p<- rep(1/people,people)
p.value<- rep(0,Ninning)
for(r in 1:Ninning){
  p.value[r]<- chisq.test(Zz[,r],p=p)[3]　　# pは確率でなくて頻度でもいい。
}                                           # そのときはrescale.p=TRUEを指定。
p.value<- unlist(p.value)
p.value

有意水準を0.001とかにしてみると

which(p.value<0.001)
[1] 1 2 3 4 5 6

6番目のイニングまでは、先頭打者になる確率は9人で等確率ではないようだ。
プロットしとくよ。

par(mfcol=c(1,3))
persp(Z,col=5,
      xlab="Batter.No",ylab="Inning.No",zlab="First batter probability",
      phi=20,theta=60)
plot(apply(Z,1,mean),type="o",xaxt="n",lwd=2,frame=F,
     xlab="No",ylab="First batter probability")
axis(1,1:people)
abline(h=1/people) 
plot(p.value,type="o",frame=F,ylim=c(0,1),xaxt="n",
     xlab="Ninning.No",ylab="p-value")
axis(1,1:Ninning)
axis(2)