Barroso, Judit
2010-Nov-10 18:19 UTC
[R] Difficult doubt about choose distances randomly in a matrix with a probability of event
I would like to build a model in R to simulate the seed dispersal by one plant. The plant produced 5 seeds and the probability of falling inside the eight closest space was 0.8 and in the next space 0.2 and in the rest space 0: 0 0 0 0 0 0 0.2 0.2 0.2 0.2 0.2 0 0.2 0.8 0.8 0.8 0.2 0 0.2 0.8 1 0.8 0.2 0 0.2 0.8 0.8 0.8 0.2 0 0.2 0.2 0.2 0.2 0.2 0 0 0 0 0 0 0 So I do not know how to program in R to choose these 5 places (randomly) knowing the probability of event. I have built a matrix that each value has assigned the probability of mortality (0.7) except when there is a plant (1) so, the 5 seeds that felt in each place around that plant would have to be multiplied by 0.7 and if the result was >1 then we would have a new individual for the next run. 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 1 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 I am trying to do loops but I am a beginner and I am very lost. Could anyone say me what code should I use? Thank you very much, Judit Barroso [[alternative HTML version deleted]]
Michael Bedward
2010-Nov-11 06:52 UTC
[R] Difficult doubt about choose distances randomly in a matrix with a probability of event
Hello Judit,
The code below is a toy simulation function that takes as arguments a
matrix representing the initial population, a dispersal kernel, global
survival probability and max number of iterations to run. It doesn't
implement exactly what you described in your post but if you study the
code you should be able to use it as a starting point.
WARNING - it was slapped together quickly so watch out for bugs.
I spend a lot of time writing plant population simulation code in R
(mainly for woodland tree species). Once you get started you will find
that it's a very nice platform for working with both simple and
complex models and / or prototyping models before moving them into
another programming language.
You might also want to post on the r-sig-ecology list.
Have fun.
Michael
plantSim <- function(pop, dispKernel, dispOrigin=NULL, nseed=5,
survival=1.0, niter=100) {
# Arguments:
# pop - matrix for initial population map
# (1 is occupied cell; 0 is vacant cell)
#
# dispKernel - matrix of probabilities for dispersal kernel
#
# dispOrigin - vector of two integers identifying the source cell of
the dispersal kernel;
# if null, the centre cell will be used
#
# nseed - number of seeds to disperse from each individual (constant)
#
# survival - individual probability of survival per time step (constant)
#
# niter - maximum number of iterations to run
#
# Returns a matrix with cols for time, number of individuals, number
of successful dispersals,
# number of deaths
result <- matrix(0, nrow=niter+1, ncol=4)
colnames(result) <- c("time", "N",
"dispersals", "deaths")
result[1, ] <- c(0, sum(pop == 1), 0, 0)
dnr <- nrow(dispKernel)
dnc <- ncol(dispKernel)
if (is.null(dispOrigin))
dispOrigin <- c(ceiling(dnr/2), ceiling(dnc)/2)
MAPROWS <- nrow(pop)
MAPCOLS <- ncol(pop)
for (iter in 1:niter) {
ndeaths <- 0
ndisp <- 0
occupied <- which(pop == 1, arr.ind=TRUE)
nstart <- nrow(occupied)
# check for extinction
if (nrow(occupied) == 0) {
# trim results matrix and finish early
result <- result[1:iter, ]
break
}
for (iocc in 1:nrow(occupied)) {
cell <- occupied[iocc, ]
# test survival
if (runif(1) >= survival) {
# mortality
pop[cell[1], cell[2]] <- 0
ndeaths <- ndeaths + 1
} else {
# survived - do seed dispersal
for (iseed in 1:nseed) {
# choose a dispersal kernel cell
repeat {
row <- sample(dnr, 1)
col <- sample(dnc, 1)
if (runif(1) < dispKernel[row, col]) {
destOffset <- c(row - dispOrigin[1], col - dispOrigin[2])
break
}
}
destCell <- cell + destOffset
# check cell is in bounds - if not the seed is lost
# (modify this to whatever boundary conditions you want)
if (destCell[1] >= 1 & destCell[1] <= MAPROWS &
destCell[2] >= 1 & destCell[2] <= MAPCOLS) {
# cell becomes occupied if vacant
if (pop[destCell[1], destCell[2]] == 0) {
pop[destCell[1], destCell[2]] <- 1
ndisp <- ndisp + 1
}
}
}
}
}
# end of time step reporting
result[iter+1, ] <- c(iter, nstart-ndeaths+ndisp, ndisp, ndeaths)
}
# return results
result
}
On 11 November 2010 05:19, Barroso, Judit
<judit.barroso at exchange.montana.edu> wrote:> I would like to build a model in R to simulate the seed dispersal by one
plant.
> The plant produced 5 seeds and the probability of falling inside the eight
closest space was 0.8 and in the next space 0.2 and in the rest space 0:
> 0
>
> 0
>
> 0
>
> 0
>
> 0
>
> 0
>
> 0.2
>
> 0.2
>
> 0.2
>
> 0.2
>
> 0.2
>
> 0
>
> 0.2
>
> 0.8
>
> 0.8
>
> 0.8
>
> 0.2
>
> 0
>
> 0.2
>
> 0.8
>
> 1
>
> 0.8
>
> 0.2
>
> 0
>
> 0.2
>
> 0.8
>
> 0.8
>
> 0.8
>
> 0.2
>
> 0
>
> 0.2
>
> 0.2
>
> 0.2
>
> 0.2
>
> 0.2
>
> 0
>
> 0
>
> 0
>
> 0
>
> 0
>
> 0
>
> 0
>
> So I do not know how to program in R to choose these 5 places (randomly)
knowing the probability of event.
>
> I have built a matrix that each value has assigned the probability of
mortality (0.7) except when there is a plant (1) so, the 5 seeds that felt in
each place around that plant would have to be multiplied by 0.7 and if the
result was >1 then we would have a new individual for the next run.
> 0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7
> 0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7
> 0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7
> 0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?1 ? ? ? ? ? ? ?0.7 ? ? ? ? ?0.7
> 0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7
> 0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7 ? ? ? ? ?0.7
>
> I am trying to do loops but I am a beginner and I am very lost. Could
anyone say me what code should I use?
>
> Thank you very much,
> Judit Barroso
>
> ? ? ? ?[[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.
>