バグ取りその1

こちらの続き。
上のリンクの関数を用意しておく。パラメータは

Nbase<- 2
people<- 9
Nout<- 3
Ninning<- 3
sup<- 18 
cut<- 0.99
hit<- 0.05   

確率計算で出た結果を

datacal<- data

に格納する。
次に、シミュレーションする。

game<- 100000                         # 15分くらい。
K<- matrix(0,people,game)
system.time(
for(b in 1:people){
  for(g in 1:game){
    K[b,g]<- one.game.score(team.array[,,b],Nout,Nbase,people,Ninning)
    #K[b,g]<- one.game.score(team.array[,,b],Nout,Nbase,people,1)
  }
}
)
FU<- function(vec){
  counter<- rep(0,max(vec))
  for(result in 0:length(counter)){
    counter[result+1]<- length(which(vec==result))
  }
  return(counter/game)
}
datasimu<- apply(K,1,FU)

ori.datacal<- datacal
Lmax<-  max(unlist(lapply(datacal,length)))
# plot用にlengthをそろえておく。
  for(l in 1:people){
    datacal[[l]]<- append(datacal[[l]],rep(0,Lmax-length(datacal[[l]])+1))
  }
Xg<- min(max(unlist(lapply(datacal,length))),
         max(unlist(lapply(datasimu,length))))
Yg<- max(unlist(lapply(datacal,max)),
         unlist(lapply(datasimu,max)))
xlim<- c(0,Xg)
ylim<- c(0,Yg)
par(mfrow=c(3,3))
for(pl in 1:people){
  plot(0:(Xg-1),datacal[[pl]],xlim=xlim,ylim=ylim,
       col=1,lwd=4,type="l",lty=1,frame=FALSE,
       xlab="score",ylab="probability")
  par(new=TRUE)
  plot(0:(Xg-1),datasimu[[pl]][1:Xg],xlim=xlim,ylim=ylim,
       col=5,lwd=4,type="l",lty=2,frame=FALSE,
       xlab="score",ylab="probability")
  title(pl)
  legend(4,Yg,c("Cal","Simu"),col=c(1,5),lwd=3,lty=1:2)
}

ずいぶん近付いた。
PHが2番から4番のとき二峰性らしい。
2番のときはふたつの峰が同じ感じを表している。
1番のときの1点のところがなんかこうクワッと盛り上がって、それで右側の波がすわわわぁぁぁと小さくなったようにも見える。