search

cboettig / regimeshifts

:package: R package providing simulations and data from ecological regime shifts
Other
3 stars 1 forks source link

Add May's model #1

Open cboettig opened 5 years ago

cboettig commented 5 years ago

@amy-tong Let's start with May's model. You can find a statement of the model on page 5 of https://github.com/cboettig/noise-phenomena/blob/master/appendixA/appendixA.pdf, (the stochastic oscillator). The code there uses nimble to simulate the model, but let's just start with code that does this as a for loop for simplicity to avoid the extra software dependency.

amy-tong commented 5 years ago

@cboettig

Seems like there are quite a few more issues than I anticipated for this code, so I don't have a coherent first draft of it just yet. Here are some questions I have:

  1. Which paper did May's model come from? None of the references listed included May as an author so I'm a bit confused. Would be great to see the paper so that I can label each variable used appropriately and better understand how the function works.
  2. Would you prefer that the may function just take in x[t], return x[t+1], and have an external code segment deal with looping (like how dai works), or should the for loop be within may?

Also, I'm not sure if the paper with May's model can answer the following questions, but if it can't:

  1. When was the y array initialized and what does y[t+1] stand for? From what I can tell it currently represents the random variable of a distribution function, but then how does x[t+1] <- max(y[t+1],0) work (since x[t+1] is a non-random variable I think)?
  2. x[1] = 1.2, but what about y[1]?

Thanks for the help!

cboettig commented 5 years ago

@amy-tong no worries, and all good questions!

  1. May (1977): https://www.nature.com/articles/269471a0 . (Sorry, you'll find the full citation in the manuscript itself, eg. https://doi.org/10.1111/ele.13085 , or preprint is in the repo. I had just linked you directly to the appendix since it has code.

  2. Good idea! Yeah, let's do it like Dai case, where iteration is done externally. We can then write a generic helper function that will take the name of a model (dai, may, etc) and do the iteration for that model.

  3. Yeah, so that block looks like R code, but is really black magic, as you may have guessed from the stochastic term and the ~. That block comes from nimble, a really cool package which doesn't actually evaluate the R code directly, but just "inspects" it, translates it into C++ code for a fully Bayesian model, and compiles that to an executable it can run. This strategy (known as "metaprogramming") is very fast and very powerful, but it's overkill here. It comes in really handy once you want to estimate the model parameters from data, but isn't really necessary for simulating. We may later add a wrapper level to let users do nimble-based estimation of these models too.

In this case y is initialized as a vector of all zeros, and note that the code actually never calls y[1] (R is 1 indexed), and y[t+1] is first used only after it is initialized from the equation above. The whole use of the y vector is really just a temporary variable separating the step between calculating the deterministic part of x[t+1] and the truncated normal part of x[t+1]. So x[t+1] is a random variable, only mu[t] is non-random (the mean or so-called deterministic skeleton).

Sounds like you are off to a good start!