打席に立ちうる人数を1イニングから1試合に拡張して考えたらPC終了のお知らせ

前回セミナーでの、シミュレーションと確率計算でバグを取りたい話の続き。
今までは
・1イニングあたりに打席に立つ打者の人数
・あるkイニングで先頭打者になる確率
を出して、それらの積的なもので重みをつけて1イニングあたりの得点確率を出していた。
そうするとシミュレーションと確率計算で結果がかなり異なるので、段階的に調査したら
　

　
これで(2)がおかしいのだろうと。
そこで、上記の方法がちょっとまずかったのだろう、ということで、いっそのこと1イニングではなく、1試合に打席に立ちうる打者人数を出してやろうという、先生手法を使うことにしてみた。
イニングリセットのアウト数がNout
1試合がNinning
で構成されるとすると、1試合に最小でNout*Ninning人が打席に立つ。
そこからNout*Ninning+1、Nout*Ninning+2、…、Nout*Ninning+k、…、人打席に立つ積算確率を、cutまで考えることにしよう。
すると計算量多すぎてPC死んだ。
　
嫁の家が壊れてしまうじゃないか!!
　

# パラメータ
library(MCMCpack)
library(gtools)
Nbase<- 2
people<- 3
Nout<- 2
Ninning<- 3
sup<- 15        # 1試合の上限人数。うまく取らないと死亡。
cut<- 0.99
hit<- 0.03      # ヒットパターンはすべて同じ確率。
# チームの作成
TEAM<- rbind(c(rep(0,Nbase),1),
             matrix(c(1-sum(rep(hit,Nbase)),rep(hit,Nbase)),
             people-1,Nbase+1,byrow=TRUE))
  b1<- mapply(rep,2:people,1:(people-1))
  b2<- mapply(rep,2:people,(people-1):1)
  b<-  diag(1,people)
  b[upper.tri(b)]<- unlist(b1); b[lower.tri(b)]<- unlist(b2)
  team.array<- array(0,c(people,Nbase+1,people))
    for(t in 1:people){
      team.array[,,t]<- TEAM[b[t,],]
    }

# プログラムに必要な関数を用意する。
# 走者パターンから得点を出す野球プログラム。
seq.score<- function(vec,batter,out.count,base,people,Inning){
  hitpattern<- c(1,lapply(mapply(rep,0,1:(base-1)),append,c(1),after=F))
  out.count<- out.count
  batter.box<- 1
  out<- 0
  score.prob<-c(0,1)
  runner<- rep(0,Nbase*Ninning)
  for(t in 1:length(vec)){
    x<- vec[t]
    score.prob[2]<- score.prob[2]*batter[batter.box,x+1] 
      if(x==0){
        out<- out+1
      }else{
        runner<- append(runner,hitpattern[[x]])
      }
    batter.box<- batter.box+1
      if(batter.box>people){
        batter.box<- 1
      }
      if(out>(out.count - 1)){
        runner<- runner[1:(length(runner)-(Nbase-1))]
        runner<- append(runner,rep(0,Nbase-1))
        out<- 0
      }
  }
  runner<- runner[1:(length(runner)-(Nbase-1))]
  score.prob[2]<- score.prob[2]*batter[batter.box,1]
  score.prob[1]<- sum(runner)
  return(score.prob)
}
# アウトとヒットになる打者を数え上げるプログラム。
CMB<- function(vec){
  hit.vec<- 1:max(out.cmb)
  hit.vec<- hit.vec[-vec]
  return(hit.vec)
}

# ここからスタート。
data<- lapply(mapply(rep,0,1:people),unique)
for(q in 1:people){
#for(q in 1:1){
  out.seq<- rep(team.array[,1,q],length=sup)
  hit.seq<- 1-out.seq
  A<- rep(0,length(out.seq))
    for(os in (Ninning*Nout):((sup)-1)){
      out.cmb<- combinations(os,Nout*Ninning-1)
        if(os==Ninning*Nout){
          hit.cmb<- as.matrix(apply(out.cmb,1,CMB))
        }else{
          hit.cmb<- t(as.matrix(apply(out.cmb,1,CMB)))
        }
          for(cmb.r in 1:nrow(out.cmb)){
            out.hit.prob<- prod(out.seq[out.cmb[cmb.r,]])*
                           prod(hit.seq[hit.cmb[cmb.r,]])*
                           out.seq[os+1]
            A[os+1]<- A[os+1]+out.hit.prob
          }
    }
  A[Nout*Ninning]<- prod(out.seq[1:(Nout*Ninning)])
  Nbatter<- sum(cumsum(A)<cut)
  GameResult<- rep(0,Nbatter-Nout*Ninning+1)
    for(batters in 1:(Nbatter-Nout*Ninning)){
      a<- p<- list(c(1:Nbase))
        if(batters != 1){
          for(cycles in 1:(batters-1)){
            p<- c(p,a)
          }
        }
    Bpatterns<- expand.grid(p)
    cmb<- combinations(Nout*Ninning+batters-1,Nout*Ninning-1) # Who is out(=0).
    Y<- matrix(0,nrow(cmb)*nrow(Bpatterns),ncol(cmb)+ncol(Bpatterns))
      for(i in 1:nrow(cmb)){
        for(t in 1:ncol(Bpatterns)){
          Y[(1+nrow(Bpatterns)*(i-1)):((nrow(Bpatterns)*i)),
            (1:ncol(Y))[-c(cmb[i,])][t]] <- Bpatterns[[t]]
        }
      }
    Y.apply.score<- apply(Y,1,seq.score,
                          team.array[,,q],Nout,Nbase,people,Ninning)
      for(sc in 0:max(Y.apply.score[1,])){
        GameResult[sc+1]<- GameResult[sc+1]+
                           sum(Y.apply.score[2,which(Y.apply.score[1,]==sc)])
      }
    }
  score.0<- head(rep(team.array[,1,q],floor(Nout*Ninning)+1),n=Nout*Ninning)
  GameResult[1]<- GameResult[1]+prod(score.0)
  GameResult<- GameResult/sum(GameResult)
  data[[q]]<- c(GameResult)
}