[R] Difficult doubt about choose distances randomly in a matrix with a probability of event

[Michael Bedward]

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.
>