[R] simulating a model

[Thomas Petzoldt]

Dear Rafael,

first of all, your simulation works, at least in a technical sense, so I 
don't understand what you mean with "can't simulate it properly".

Second, your SIR-based model is a quite different from the SIR models I 
know (e.g. http://en.wikipedia.org/wiki/SIR_Model).

The R code, however, seems to be technically correct, if compared with 
your system of equations. To help you solving your problem, we need more 
information, e.g. how the equations where derived, where the parameters 
come from, what is the process behind, and, most important, why do you 
think that the outcome is wrong.

In addition, I guess that your simulation time is too long, compared 
with the speed of the process. Try something like

times = c(from=0, to=1, by=0.01)


Thomas



Rafael Moral wrote:
> Dear useRs,
> 
> I have written an ecological model, based on the epidemiology SIR model.
> I've been trying to simulate it in R.
> However, I can't simulate it properly.
> Two guesses: my script isn't right; I'm not setting the parameters properly
> 
> I have uploaded an image to the model here:
> http://img24.imageshack.us/img24/743/imagemutr.jpg
> 
> The script I am using is as it follows:
> 
> require(simecol)
> 
> mod1 <- new("odeModel",
>   main = function(time, init, parms) {
>   x <- init
>   p <- parms
>     dx1 <- p["K"] - p["alpha"]*x[1]*x[2] - p["gamma"]*x[1]
>     dx2 <- x[1]*x[2]*(p["alpha"] - p["beta"])
>     dx3 <- p["beta"]*x[1]*x[2] + p["gamma"]*x[1]
>     list(c(dx1, dx2, dx3))
> },
> times = c(from=0, to=100, by=0.1),
> parms = c(K=100, alpha=0.3, gamma=0.5, beta=0.2),
> init = c(S=500, V=100, R=0),
> solver = "lsoda"
> )
> 
> plot(sim(mod1))
> 
> Thanks in advance!
> Rafael.
> 
> 
>       ____________________________________________________________________________________
> [[elided Yahoo spam]]
> 
> 	[[alternative HTML version deleted]]
> 
> 
> 
> ------------------------------------------------------------------------
> 
> ______________________________________________
> 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.