[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
While preparing for the new semester I have started reimplementing standard NetLogo examples in R. The first is Fire model.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The simulation in R is presented here:
# Forest matrix trees encoding:< o:p>
# 3 – green, 2 – burning, 1 – burnt, 0 – no tree< o:p>
simulate <- function (forest.density = 0.6, size = 251) {< o:p>
forest <- matrix(3 * rbinom(size ^ 2, 1, forest.density),< o:p>
nrow = size)< o:p>
forest[1,] <- 2< o:p>
while (sum(forest == 2) > 0) {< o:p>
list.red <- which(forest == 2)< o:p>
# find spots touching burning trees< o:p>
ignite.x <- 1 + floor((list.red – 1) / size)< o:p>
ignite.x <- rep(ignite.x, each = 4) +< o:p>
rep(c(1, –1, 0,0), length(ignite.x))< o:p>
ignite.y <- 1 + ((list.red – 1) %% size)< o:p>
ignite.y <- rep(ignite.y, each = 4) +< o:p>
rep(c(0, 0, 1,–1), length(ignite.y))< o:p>
include <- (ignite.x > 0) & (ignite.x <= size) &< o:p>
(ignite.y > 0) & (ignite.y <= size)< o:p>
ignite.n <- (ignite.x[include] – 1) * size +< o:p>
ignite.y[include]< o:p>
# ignite green trees< o:p>
is.green <- (forest[ignite.n] == 3)< o:p>
forest[ignite.n[is.green]] <- 2< o:p>
# trees burn for one period< o:p>
forest[list.red] <- 1< o:p>
}< o:p>
return(1 – sum(forest[-1,] == 3) / sum(forest[-1,] > 0))< o:p>
}< o:p>
The only thing left was to check that the results produced by NetLogo and R implementations are identical. I have simulated both models 1000 times for forest density equal to 60%. Here is the plot comparing percent burned density estimates and Kolmogorov-Smirnov test p-value.
It can be concluded that for the parametrization I have selected models produce similar results.
To leave a comment for the author, please follow the link and comment on their blog: R snippets.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.