In cboettig/nonparametric-bayes: Nonparametric Bayes inference for ecological models

library(pomp)
#opts_knit$set(upload.fun = socialR::flickr.url)
#options(flickr_tags = "nonparametric-bayes")

process.sim <- function(x, t, params, delta.t, ...){
# r = 0.75, K = 10, Q = 3, H = 1, a = 1.3 
  r <- params["r"]
  K <- params["K"]
  Q <- params["Q"]
  H <- params["H"]
  a <- params["a"]
  sigma <- params["sigma"]
  X <- x["X"]
  epsilon <- exp(rlnorm(n = 1, mean = 0, sd = sigma))  

  c(X = unname(
    X * exp(r * delta.t * (1 - X / K) - a * X ^ (Q - 1) / (X ^ Q + H ^ Q)) * epsilon)) 
}

measure.sim <- function(x, t, params, ...){
  tau <- params["tau"]
  X <- x["X"]
  y <- c(Y = unname(rlnorm(n = 1, meanlog = log(X), sd=tau)))
}

Create a container of class pomp

may <- pomp(data = data.frame(time = 1:100, Y = NA), times = "time", rprocess = discrete.time.sim(step.fun = process.sim, delta.t = 1), rmeasure = measure.sim, t0 = 0)

Parameters and an initial condition

theta <- c(r = 0.75, K = 10, Q = 3, H = 1, a = 1.3, 
           sigma = 0.2, tau = 0.1, X.0 = 5) 

Simulate

dat <- simulate(may, params=theta, obs=TRUE, states=TRUE)
library(reshape2)
library(ggplot2)
dt <- melt(dat)[3:5]
names(dt) = c("time", "value", "variable")
ggplot(dt, aes(time, value, lty=variable)) + geom_line()

Note that simulate returns an instance of our object, with the data column filled in. This is not true if we give it arguments such as state or obs = TRUE as above. However, when we query for observed and true states in this way, the simulated data is not stored in our object. We must assign it with may <- pomp(data = ... where data is a data frame with columns "time" and "Y". This is much simpler if we use the default return:

may <- simulate(may, params=theta)

Measurement density

measure.density <- function(y, x, t, params, log, ...) {
    tau <- params["tau"]
    ## state at time t:
    X <- x["X"]
    ## observation at time t:
    Y <- y["Y"]
    ## compute the likelihood of Y|X,tau
    f <- dlnorm(x = Y, meanlog = log(X), sdlog = tau, log = log)

}

Add to our container

may <- pomp(may, dmeasure = measure.density)

Optionally, work in the transformed variable space instead by providing parameter and inverse parameter transforms.

may  <- pomp(may, parameter.transform = function(params, 
    ...) {
    exp(params)
}, parameter.inv.transform = function(params, ...) {
    log(params)
})

Run the particle filter. Note that this evaluates the likelihood of our simulated data (stored in the object) at the given (true) parameters.

pf <- pfilter(may, params = theta, Np = 1000)
logLik(pf)

Here we go, estimate the model by particle filter. This'll be slow...

theta.guess <- theta + c(abs(rnorm(length(theta)-1)),0)
estpars <- names(theta[1:7]) # which pars should be estimated. lets try and get them all
rw.sd <- theta[1:7]/(50*theta[1:7]) # a named vector of random walk values. This sets it to 0.02 for all pars
system.time(mf <- 
  mif(may, Nmif = 10, # Number of iterations
      start = theta.guess, # Starting value of parameters
      transform = TRUE, # work in transformed parameter space (e.g. logs avoid neg par values)
      pars = estpars,  # which parameters to estimate (others fixed at starting value)
      rw.sd = rw.sd, # random walk width
      Np = 200, # Number of particles
      var.factor = 4, # scaling factor for random walk sd
      ic.lag = 10, # fixed lag
      cooling.factor = 0.999, # exponential cooling
      max.fail = 10) # internal tolerance
)

compare.mif(mf)

coef(mf)
logLik(mf)
pf_est <- pfilter(may, params = theta.guess, Np = 1000)
logLik(pf_est) # initial loglik

These functions apply to replicated mifs above.

theta.true <- coef(may)
theta.mif <- apply(sapply(mf,coef),1,mean)
loglik.mif <- replicate(n=10,logLik(pfilter(mf[[1]],params=theta.mif,Np=10000)))
bl <- mean(loglik.mif)
loglik.mif.est <- bl+log(mean(exp(loglik.mif-bl)))
loglik.mif.se <- sd(exp(loglik.mif-bl))/exp(loglik.mif.est-bl)

cboettig/nonparametric-bayes documentation built on May 13, 2019, 2:09 p.m.