[R] [help] simulation of a simple Marcov Stochastic process for population genetics

[Ben Bolker]

z3mao <zhangyij <at> gmail.com> writes:

> 
> 
> Hi, this is my first time using R. I want to simulate the following process:
> "in a population of size N, there are i individuals bearing genotype A, the
> number of those bearing A is j in the next generation, which following a
> binominal distribution (choose j from 2*N, the p is i/2*N), to plot the
> probability of the next generations, my script is as follows. It cannot run
> successfully, declaring that the "ylim should be limited. " I wonder where
> the bug is. Thanks very much!
> 
> generation<-function(i,N)
> {
>     m<-1;gen<-numeric(); 
>     for(m in 1:50)
>   {
>    testp<-runif(1,0,1);
>     j<-0; sump<-0;
>     while(sump < testp)
>       {  sump<-sump+dbinom(j,2*N,i/(2*N));
>          j<-j+1;
>        }
>      i<-j; gen[m]<-j/(2*N); m<-m+1; 
>    }
>   plot(m, gen[m]); 
> }



## basic sim
N <- 50
ngen <- 500
nA <- numeric(ngen)
nA[1] <- 25
for (i in 2:ngen) {
  fixed <- nA[i-1]==0 || nA[i-1]==N
  if (fixed) break
  nA[i] <- rbinom(1,prob=nA[i-1]/N,size=N)
}

## wrap it in a function
binsim <- function(N,ngen,A0) {
  nA <- numeric(ngen)
  nA[1] <- A0
  for (i in 2:ngen) {
    fixed <- nA[i-1]==0 || nA[i-1]==N
    if (fixed) break
    nA[i] <- rbinom(1,prob=nA[i-1]/N,size=N)
  }
  nA
}

## replicate it
m <- replicate(200,binsim(50,100,25))
matplot(m,col="grey",type="l",lty=1,pch=1)
means <- rowMeans(m)
quants <- t(apply(m,1,quantile,c(0.025,0.975)))
lines(means,lwd=2)
lines(rowMeans(m)-,lwd=2)
matlines(quants,lty=2,lwd=2,col=1)