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

[Dan Davison]

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