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)
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
