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

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]); 
}
-- 
View this message in context:
http://www.nabble.com/-help--simulation-of-a-simple-Marcov-Stochastic-process-for-population-genetics-tp19085705p19085705.html
Sent from the R help mailing list archive at Nabble.com.

On Thu, Aug 21, 2008 at 03:00:51AM -0700, z3mao wrote:
> 
> 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. "
In a situation like this, try using options(error=recover) to debug.
> I wonder where > the bug is. Thanks very much!
There are several bugs... The most serious is that your homemade
binomial random number generator is wrong. (For example, look at what
happens when it is given a probability parameter of 0: it returns 1
rather than 0. Your alleles aren't going to be lost from the
population very often!). So, if someone has set you the task of
simulating drift without using a built-in binomial RNG, then you'll
need to think through your RNG code again. But if you are free to do
what you want, then you should use the R function rbinom to generate
binomial RVs.

Here are comments on the other bugs with a cleaned up (but still
probabilistically wrong) version below.
> 
> generation<-function(i,N)
> {
>     m<-1;
## Don't initialise m here; it gets initialised in the for loop
> gen<-numeric(); 
##  gen <- rep(NA, 50) is better
>     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've already said that the above is wrong
>      i<-j; gen[m]<-j/(2*N);
m<-m+1; 

## The for loop deals with incrementing m; don't do it yourself!

>    }
>   plot(m, gen[m]); 
## You want plot(1:50, gen, type="l")

## You don't need semicolons at the end of lines in R!
> }
Here's a version of your code that corrects the other bugs, but
still has your incorrect binomial RNG code in it.

generation <- function(i,N)
{
    warning("binomial RNG code is wrong")
    mvals<- 1:50;
    gen<- numeric();
    for(m in mvals)
    {
        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(mvals, gen, type="l");
}

Dan



> -- 
> View this message in context:
http://www.nabble.com/-help--simulation-of-a-simple-Marcov-Stochastic-process-for-population-genetics-tp19085705p19085705.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
-- 
http://www.stats.ox.ac.uk/~davison

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)

Thank you so much!

-- 
View this message in context:
http://www.nabble.com/-help--simulation-of-a-simple-Marcov-Stochastic-process-for-population-genetics-tp19085705p19101425.html
Sent from the R help mailing list archive at Nabble.com.