1/17 MIKUセミナー

MikuHatsune2011-01-17

自分のテーマとして、野球の続き。
・シミュレーションと確率過程での計算が、ずれているのか、同じ傾向(状況?経過?)を示しているのかどうか。
・異なる打撃確率を持つpeople人で1チームが構成されている場合に、どう確率過程を計算するか。
・次のテーマ
について。
 
野球のシミュレーション確率計算が合っているような合っていないようなよくわからない話。
1イニングあたりに打席に立つ人数を何%まで考えるかで、かなり正確さが違う。
セミナー前までは95%で考えていたが、1時間くらいPC酷使してなんとか99%まで計算させた。
シミュレーションも100000試合させた。
セミナーで提示したときよりもずいぶん近付いてきたぞ…?

ただ、これもよくわからないので、自分で振ったパラメータ(イニング数、アウト数、打率、塁数)を簡単にしてバグ取りしよう。
1イニングあたりの人数を切らない方法もなんか思いつけたらうれしい。
 
1チーム1人から、複数に拡張したい。
シミュレーションなら、PCにがんばってもらえば得点分布が得られる。これを確率計算で出すにはどうしようか。
とりあえず、1イニングに打席に立ちうる人数を用いて、誰が何回先頭打者になりうるかをシミュレーションしてみた。

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.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,length(vec))
  for(result in 1:length(vec)){
    count1to.[result]<- length(which(vec==result))
  }
  return(count1to.)
}

counter1<- rep(0,people)
Bnum<- rep(1:people,Ninning*min(which(cumsum(A)>cut)))
for(I in 1:10000){
  k<- sample(c(1:length(bboxesprob)),prob=bboxesprob,size=Ninning-1,replace=TRUE)
  j<- Bnum[c(1,cumsum(k)+1)]
  counter1<- counter1+countfunc(j)
}
plot(counter1/sum(counter1),type="o",xaxt="n",lwd=2,frame=F,
     xlab="No",ylab="First batter probability")
axis(1,1:people) 


10000試合したら、1番打者は10000試合の1イニング目は確実に先頭打者。
イニングが無限大になれば、どの打者も先頭打者になる確率は均一になる。
1イニングと2イニングでの先頭打者になりやすさのわずかな差が、妙な山を作っている。
ということでたぶん漸近線が\frac{1}{people}のところに近付く、と思う。
これで、先頭打者が1番のとき、2番のとき、…とで何点取るか考えよう、とか思ったけど
Rの特性として、sampleで文字列を発生させるのは超速いから、ここから考えるというのは効率がいい。
まあなんとかやってみたいと思う。
 
次のテーマとして、ある組織中の細胞を考えて、「2つの細胞の存在比率」をモデル化、シミュレーションしようと思う。・時間(変化)
・(今のところとりあえず)2つの細胞の相互作用
・相手の増殖速度に影響する
を考える。2つの細胞は、決められたある空間内に存在するとして、各々が増殖や相互作用しあうことで最終的に2つの細胞はどうなるか(どんな比率で存在するか)。
最初は、外部からの影響が無い状況でやって、うまくモデルができたら、例えば癌細胞と正常細胞を考えて、抗癌剤が聞いて癌細胞が死滅した、とか、逆に効き目がなくて正常細胞が死にました、とか
そうしたら正常細胞が死にすぎて、臓器が正常な機能を保てなくなり、個体の死に至りました、くらいにまで拡張できたら〜
と思う。
 
最初の図のスクリプト。起動させたらスペックによっては1時間くらいかかる。かも。
パラメータ。

library(MCMCpack)
Nbase<- 4
Nout<- 3
Ninning<- 9
Bprob<- matrix(c(1/2,1/8,1/8,1/8,1/8),1,5)
hitp<- 1-Bprob[1] 
Vnum<- 300
game<- 100000   
cut<- c(0.5,0.8,0.85,0.9,0.95,0.99)

自作関数を用意。

num.count<- function(vec){
  counter<- rep(0,Nbase)
  for(result in 1:Nbase){
    counter[result]<- length(which(vec==result))
  }
  return(c(Nout,counter))
}
score<- function(vec){
  runner<- rep(0,Nbase)
  for(k in 1:length(vec)){
      runner<- c(runner,c(1,rep(0,(vec[k] - 1))))
  }
  runner<- runner[1:(length(runner) - (Nbase - 1))]
  return(sum(runner))
}
Prob<- function(vec){
  P<- choose((Nout+length(vec)-1),(Nout-1))*prod((Bprob^num.count(vec)))
  return(c(score(vec),P))
}
library(expm)
one.game.score<-
        function(batter,out.count){   
          si  <- c(1);dou <- c(1,0);tri <- c(1,0,0);hr  <- c(1,0,0,0)
          out.count<- out.count
          total.score<- rep(0,9)
          batter.box<- 1
          for(t in 1:9){
            runner<- rep(0,9)
            out<- 0
            repeat{
              x<- sample(c(1,2,3,4,5),prob=batter[batter.box,],size=1)
              switch(x,
                     (out<- out + 1),
                     (runner<- append(runner,si)),
                     (runner<- append(runner,dou)),
                     (runner<- append(runner,tri)),
                     (runner<- append(runner,hr))
              )
              switch(batter.box,
                     (batter.box<- 2),(batter.box<- 3),(batter.box<- 4),
                     (batter.box<- 5),(batter.box<- 6),(batter.box<- 7),
                     (batter.box<- 8),(batter.box<- 9),(batter.box<- 1)
              )
              if(out>(out.count - 1)){
                break
              }
            }
            total.score[t]<- sum(runner[c(1:(length(runner)-3))])
          }
          return(sum(total.score))
        }
hist.data<-
function(n,batter,out.count){                            
  data1<- rep(0,n)    
  for(i in 1:n){
    data1[i]<- one.game.score(batter,out.count)
  }
  result.count<- rep(0,(max(data1) + 1))
  for(p in 0:max(data1)){
    result.count[p]<- length(which(data1 == p))
  }
  return(data1)                                 
}

ここからスタート。

R<- Rr<- list(0)
for(r in 1:length(cut)){
  R<- c(R,Rr)
}
for(cc in 1:length(cut)){
  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)
  }
  cut.off<- min(which(cumsum(A)>cut[cc]))
  Nbatter<- cut.off
  Result<- rep(0,(Nbatter-Nout+1))
    a<- p<- list(c(1:Nbase))
    Bpatterns<- expand.grid(p)
    PROBABILITY<- apply(Bpatterns,1,Prob)
    for(w in 0:max(PROBABILITY)){
      Result[w+1]<- Result[w+1]+sum(PROBABILITY[2,which(PROBABILITY[1,]==w)])
    }
  for(b in 1:(Nbatter-Nout-1)){
    p<- c(p,a)
    Bpatterns<- expand.grid(p)
    PROBABILITY<- apply(Bpatterns,1,Prob)
    for(w in 0:max(PROBABILITY)){
      Result[w+1]<- Result[w+1]+sum(PROBABILITY[2,which(PROBABILITY[1,]==w)])
    }
  }
  Result[1]<- Result[1]+Bprob[1]^Nout
  Result<- Result/sum(Result)
  M<- matrix(0,                                             
             nr=(length(Result)-1)*Ninning+length(Result),
             nc=(length(Result)-1)*Ninning+1)
  M[1:length(Result),1]<- Result 
  m<- matrix(M[1:ncol(M),1],nc=1)
  for(j in 2:ncol(M)){
    M[j:(length(Result)+j-1),j]<- Result
  }
  M<- head(M,n=ncol(M))
  Rinning<- M %^% (Ninning-1) %*% m        
  R[[cc]]<- Rinning
}

Bprob1<- matrix(Bprob,9,5,byrow=T)

D<- rep(0,game)
for(i in 1:game){
  D[i]<- one.game.score(Bprob1,Nout)
}
Q<- rep(0,max(D)+1)
for(t in 0:max(D)){
  Q[t+1]<- sum(D==t)
}
Q<- Q/length(D)
R[[length(R)]]<- Q

Mx<- c(0,length(R[[length(R)]]))
Mx<- c(0,60)    
My<- c(0,max(unlist(lapply(R,max))))
color<- c(2:length(R),1)
for(p in 1:length(R)){
  if(p==3){L<- 2}else{L<- 1}
    plot(R[[p]],col=color[p],lwd=3,lty=L,type="l",
         xlab="",ylab="",xlim=Mx,ylim=My,frame=F)
    par(new=T)
}
axis(1:2)
title("Calculation and Simulation")
legend(40,max(My),c(cut,"Simulation"),col=color[1:length(R)],lty=1,lwd=3)
text(45,0.08,"Probability is (1/2,1/8,1/8,1/8,1/8).")