状態遷移

有限個の状態の離散的世代経過を状態間遷移関数でシミュレートすれば、遺伝的ドリフトを試せる。

#染色体N本中p本が着目アレル
#それがk倍になった後に、N本を抜き出し
#そこにi本の着目アレルがある確率

probDrift<-function(N,p,k,i){
 kN<-k*N
 kp<-k*p
 kNp<-kN-kp
 k1N<-(k-1)*N
 n11<-i
 n12<-kp-n11
 n21<-N-n11
 n22<-k1N-n12
 return(exp(lgamma(kp+1)+lgamma(kNp+1)+lgamma(k1N+1)+lgamma(N+1)-lgamma(kN+1)-lgamma(n11+1)-lgamma(n12+1)-lgamma(n21+1)-lgamma(n22+1)))
}

#染色体N本中p本が着目アレル
#それがk倍になった後に、N本を抜き出し
#そこに0-N本の着目アレルがある確率

probDrift02N<-function(N,p,k){
 ret<-rep(0,(N+1))
 kN<-k*N
 kp<-k*p
 kNp<-kN-kp
 k1N<-(k-1)*N
 commonLN<-lgamma(kp+1)+lgamma(kNp+1)+lgamma(k1N+1)+lgamma(N+1)-lgamma(kN+1)
 for(i in 0:N){
  if(kp>=i){
   n11<-i
   n12<-kp-n11
   n21<-N-n11
   n22<-k1N-n12
   if(n11>=0 && n12>=0 && n21>=0 && n22>=00){
    ret[i+1]<-exp(commonLN-lgamma(n11+1)-lgamma(n12+1)-lgamma(n21+1)-lgamma(n22+1))
   }

  }
 }
 return(ret)
}
probDrift02NInf<-function(N,p,k,infty=FALSE){
 if(infty){
  ret<-rep(0,(N+1))
  pr1<-p/N
  pr2<-1-pr1
  if(pr1==0){
   ret[N+1]<-1
  }else if(pr2==0){
   ret[1]<-1
  }else{
   for(i in 0:N){
    ret[i+1]<-exp(lgamma(N+1)-lgamma(i+1)-lgamma(N-i+1)+i*log(pr1)+(N-i)*log(pr2))
   }
  }
  return(ret)
 }else{
  ret<-rep(0,(N+1))
  kN<-k*N
  kp<-k*p
  kNp<-kN-kp
  k1N<-(k-1)*N
  commonLN<-lgamma(kp+1)+lgamma(kNp+1)+lgamma(k1N+1)+lgamma(N+1)-lgamma(kN+1)
  for(i in 0:N){
   if(kp>=i){
    n11<-i
    n12<-kp-n11
    n21<-N-n11
    n22<-k1N-n12
    if(n11>=0 && n12>=0 && n21>=0 && n22>=00){
     ret[i+1]<-exp(commonLN-lgamma(n11+1)-lgamma(n12+1)-lgamma(n21+1)-lgamma(n22+1))
    }

   }
  }
  return(ret)
 }

}
nextGenExProb2<-function(x,k,infty){
 N<-length(x)-1
 ret<-rep(0,length(x))
 for(i in 1:length(x)){
  p<-i-1
  tmpret<-rep(0,length(x))
  tmpret<-probDrift02NInf(N,p,k,infty=infty)
  ret<-ret+tmpret*x[i]
 }
 return(ret)
}

DriftSim3<-function(k=2,initial=1,np=20,ng=100,infty=FALSE){
 m<-matrix(rep(0,(np*2+1)*ng),nrow=ng)
 m[1,initial+1]<-1
 for(i in 2:ng){
  m[i,]<-nextGenExProb2(m[i-1,],k=k,infty=infty)

 }
 return(m)
}


out<-DriftSim3(np=20,initial=5,ng=10,infty=FALSE)

persp(out,theta=-50,phi=60,shade=0.5)