September 14, 2017
\[X_0 \overset{\theta}{\rightarrow} X_1 \overset{\theta}{\rightarrow} X_2 \overset{\theta}{\rightarrow} \cdots\]
\[\begin{bmatrix}A_0\\B_0\\C_0\end{bmatrix} \overset{\theta}{\rightarrow} \begin{bmatrix}A_1\\B_1\\C_1\end{bmatrix} \overset{\theta}{\rightarrow} \begin{bmatrix}A_2\\B_2\\C_2\end{bmatrix} \overset{\theta}{\rightarrow} \cdots\]
## X.bull. X.bear. X.stagnant. ## 1 1 0 0 ## 2 0 1 0 ## 3 1 0 0 ## 4 0 0 1 ## 5 1 0 0 ## 6 0 1 0
\[\begin{array} \text{B} & B & B & B & 1 & 2 & 3 & 4 & B & \text{Signal} \\ A & C & G & A & A & C & T & C & A & \text{Observations} \end{array}\]
A Ricker model is a population model that relates the population at time \(t\) to the population at time \(t + 1\). \[N_{t+1} =r N_t e^{\left(-\frac{N_t}{k}\right)}\]
\(r\) is the intrinsic growth rate
\(k\) is the carrying capacity
If we assume that the growth rate \(r\) is log-normally distributed, then \[N_{t+1} = r N_t \exp(-N_t/k + \epsilon_t),\] where \(\epsilon_t \sim \text{Normal}(0, \sigma)\).
In any real system, we need to model measurement error: \[y_t \sim \text{Poisson}(\phi N_t)\]
library(ggplot2) data <- read.csv("./data/rickerdata.csv") ggplot(data, aes(time, pop)) + geom_line()
library(pomp) rickerModel <- pomp(data = data, times = 'time', t0 = 0)
rickerModel <- pomp(rickerModel, statenames = c('N'), obsnames = c('pop'), paramnames = c('r','k','sigma'))
init <- function (params, t0,...){ return(c('N' = params[['N.0']] )) } rickerModel <- pomp(rickerModel, initializer = init)
rproc <- discrete.time.sim( step.fun = function(x, t, params, delta.t, ...) { sigma <- params[['sigma']] r <- params[['r']] k <- params[['k']] N <- x[['N']] e <- rnorm(1, mean=0, sd=sigma) result <- c('N' = r * N * exp(e - N/k)) return(result) }, delta.t = 1 ) rickerModel <- pomp(rickerModel, rprocess = rproc)
\(y_t \sim \text{Poisson}(\phi N_t)\)
rmeas <- function (x, t, params, ...) { N <- x[['N']] phi <- params[['phi']] result <- c('pop' = rpois(n=1, lambda=phi * N)) return(result) } rickerModel <- pomp(rickerModel, rmeasure = rmeas)
dmeasure
should be the PDF of the Poisson distribution we used for rmeasure
. This is available in R as the dpois
function.dmeas <- function(y, x, t, params, log, ...) { N <- x[['N']] pop <- y[['pop']] phi <- params[['phi']] return(dpois(pop, lambda = phi*N, log=log)) } rickerModel <- pomp(rickerModel, dmeasure = dmeas)
p <- c(N.0=1, r=20, k=1,sigma=0.1,phi=200) sim <- simulate(rickerModel, params = p, as.data.frame = T) ggplot(sim, aes(time, pop)) + geom_line()
sim <- simulate(rickerModel, params = p, nsim=8, as.data.frame=TRUE,include.data=TRUE) ggplot(data=sim,aes(x=time,y=pop,group=sim,color=(sim=="data")))+ geom_line()+guides(color=FALSE)+ facet_wrap(~sim,ncol=3)