library(greta.dynamics)
greta_sitrep()ℹ checking if python available✔ python (version 3.7) availableℹ checking if TensorFlow available✔ TensorFlow (version 1.14.0) availableℹ checking if TensorFlow Probability available✔ TensorFlow Probability (version 0.7.0) availableℹ checking if greta conda environment available✔ greta conda environment availableℹ greta is ready to use!This example is taken from the ode_solve.R helpfile, and
run here to demonstrate the outputs from this in documentation.
# replicate the Lotka-Volterra example from deSolve
library(deSolve)
LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion <- rIng * Prey * Predator
    GrowthPrey <- rGrow * Prey * (1 - Prey / K)
    MortPredator <- rMort * Predator
    
    dPrey <- GrowthPrey - Ingestion
    dPredator <- Ingestion * assEff - MortPredator
    
    return(list(c(dPrey, dPredator)))
  })
}
pars <- c(
  rIng = 0.2, # /day, rate of ingestion
  rGrow = 1.0, # /day, growth rate of prey
  rMort = 0.2, # /day, mortality rate of predator
  assEff = 0.5, # -, assimilation efficiency
  K = 10
) # mmol/m3, carrying capacity
yini <- c(Prey = 1, Predator = 2)
times <- seq(0, 30, by = 1)
out <- ode(yini, times, LVmod, pars)
# simulate observations
jitter <- rnorm(2 * length(times), 0, 0.1)
y_obs <- out[, -1] + matrix(jitter, ncol = 2)# fit a greta model to infer the parameters from this simulated data
# greta version of the function
lotka_volterra <- function(y, t, rIng, rGrow, rMort, assEff, K) {
  Prey <- y[1, 1]
  Predator <- y[1, 2]
  
  Ingestion <- rIng * Prey * Predator
  GrowthPrey <- rGrow * Prey * (1 - Prey / K)
  MortPredator <- rMort * Predator
  
  dPrey <- GrowthPrey - Ingestion
  dPredator <- Ingestion * assEff - MortPredator
  
  cbind(dPrey, dPredator)
}
# priors for the parameters
rIng <- uniform(0, 2) # /day, rate of ingestion
rGrow <- uniform(0, 3) # /day, growth rate of prey
rMort <- uniform(0, 1) # /day, mortality rate of predator
assEff <- uniform(0, 1) # -, assimilation efficiency
K <- uniform(0, 30) # mmol/m3, carrying capacity
# initial values and observation error
y0 <- uniform(0, 5, dim = c(1, 2))
obs_sd <- uniform(0, 1)
# solution to the ODE
y <- ode_solve(lotka_volterra, y0, times, rIng, rGrow, rMort, assEff, K)
# sampling statement/observation model
distribution(y_obs) <- normal(y, obs_sd)# we can use greta to solve directly, for a fixed set of parameters (the true
# ones in this case)
values <- c(
  list(y0 = t(1:2)),
  as.list(pars)
)
vals <- calculate(y, values = values)[[1]]
plot(vals[, 1] ~ times, type = "l", ylim = range(vals))
lines(vals[, 2] ~ times, lty = 2)
points(y_obs[, 1] ~ times)
points(y_obs[, 2] ~ times, pch = 2)# or we can do inference on the parameters:
# build the model (takes a few seconds to define the tensorflow graph)
m <- model(rIng, rGrow, rMort, assEff, K, obs_sd)
# compute MAP estimate
o <- opt(m)
o$par
$par$rIng
[1] 0.2053214
$par$rGrow
[1] 1.030247
$par$rMort
[1] 0.193122
$par$assEff
[1] 0.4651422
$par$K
[1] 9.913288
$par$obs_sd
[1] 0.1036765
$value
[1] -44.57787
$iterations
[1] 55
$convergence
[1] 0