[R] Newbie q. need some help understanding this code.

[sean kim]

dear all. 

Would someone be kind and willing to explain the code
below for a person who has never used R?  ( that is if
one has enough time and inclination) 

It implements gillepsie's stochastic algorithm for
Lotka Volterra model. 

What would help me tremendously is to see the
breakdown of the line by line code into plain english.


thanks for any insights or other comments. 

sean 



library(stepfun)

lv <- function(N=1000,cvec=c(1,0.005,0.6),x=c(50,100))
{
        m<-length(cvec)
        n<-length(x)
        xmat<-matrix(nrow=N+1,ncol=n)
        tvec<-vector("numeric",N)
        h<-vector("numeric",m)
        t<-0
        xmat[1,]<-x
        for (i in 1:N) {
                h[1]<-cvec[1]*x[1]
                h[2]<-cvec[2]*x[1]*x[2]
                h[3]<-cvec[3]*x[2]
                h0=sum(h)
                tp<-rexp(1,h0)
                t<-t+tp
                u<-runif(1,0,1)
                if ( u < h[1]/h0 ) {
                        x[1] <- x[1]+1
                } else if ( u < (h[1]+h[2])/h0 ) {
                        x[1] <- x[1]-1
                        x[2] <- x[2]+1
                } else {
                        x[2] <- x[2]-1
                }
                xmat[i+1,]<-x
                tvec[i]<-t
        }
       
list(stepfun(tvec,xmat[,1]),stepfun(tvec,xmat[,2]))
}



results <- lv(N=10000)


op<-par(mfrow=c(2,1))
plot(results[[1]],do.points=True)
plot(results[[2]],do.points=False)
par(op)