0

Apologies in advance if this question has already been asked, but I can't seem to find anything online to help!

I'm trying to add stochasticity to a function I've generated, that contains a logistic growth model.

This is my function:

    ricker <- function(K, R, B0, t, E, Q)
    {
  # create empty matrix with only biomass of first year; this is where all         other values will be added
  B <- matrix(nrow = t, ncol = length(Q))
  B[1,1:length(Q)] <- B0

  for (i in 1:t) {
    for (j in 1:length(Q)) {
      B[i + 1,j]  <- 
{
  B[i,j] + R * B[i,j] * (1 - B[i,j]/K) - Q[j] * E * B[i,j]

  }
}
}
return(B)
}

The matrix created above has 34 rows and 34 columns. The following is the code I'm using to attempt to add stochasticity to my data.

# create variables/parameters
k <- 420000 
b0 <- as.matrix(3363) 
T <- 34 # number of years in data 
e <- fp 
r=1.2

pb <- txtProgressBar(style=3)

for (loop in 1:10) {

  B <- ricker(K=k, R=r, B0=b0, t=T, E=0.5, Q=fp)

  B[loop] <- B[loop] - rnorm(1,0,0.25)

  if (B[loop] <0) {

   B [loop] <- 0
  }
  setTxtProgressBar(pb,value=l/T)
}
#'
matplot(B,type="l", col=1)

plot of output from function

I want to end up with graph like this:

enter image description here

However, the output isn't any different to the initial output and I can't figure out what I'm missing. Any help on the matter would be greatly appreciated!!

Thanks!!

Brooke
  • 41
  • 1
  • 4
  • What exactly are you asking here? What does it mean to "add stochasticity"? You just mean subtract a random normal value? The random values you are generating are so small compared to the values of B you are plotting. Nearly all will be <1 so it may not look different on a plot. What are you using to verify your results? – MrFlick Feb 24 '16 at 17:44
  • @MrFlick, I've uploaded an example of the kind of graph I'm attempting to get out of my function. I've tried experimenting with simpler code, and by increasing the rnorm values, but I can't seem to get my populations to fluctuate as I want. – Brooke Feb 24 '16 at 21:09
  • What R norm values did you try? How about `rnorm(1, 3000, 2000)`? – MrFlick Feb 24 '16 at 21:10

0 Answers0