Tuesday, January 22, 2013

Fun with vectorizing

I found myself staying up way too late and having way too much fun helping a friend with her R homework. One of the things we did was to vectorize a function! Her assignment was to create plots using the Ricker model, with varying values for r:
So the first thing we did was to write a function that would calculate the population size over a certain time period for a single value of r:

Ricker = function(r, K = 100, n0 = 50, tf = 100) {
    n = numeric(length = tf)
    n[1] = n0
    for(t in 1:(tf-1)) {
          # We have tf - 1 because when t = 1, we are actually
          # calculating n[2], so if we want to end on n[tf], t 
          # should end on tf - 1
         n[t+1] = n[t]*exp(r*(1-n[t]/K))
    }
    n
}

Changing this to work with a longer vector of r values was fairly simple:

vectorizedRicker = function(r, K = 100, n0 = 50, tf = 100) {
     NumofRs = length(r)
      # First, we're changing our output to be
      # a matrix instead of a vector
     n = matrix(0, NumbofRs, tf)
      # We could vectorize over n0 as well, or any of the 
      # other parameters, if we set the number of rows for n
      # to be equal to the length of the longest vector
      # of parameter inputs
     n[,1] = n0
     for(t in 1:(tf-1)) {
           # Now, we just change the for loop to reflect that it
           # is working on entire columns of a matrix rather
           # than just a single element of a vector
          n[,t+1] = n[,t]*exp(r*(1 - n[,t]/K))
     }
     n
}

1 comment: