Examples

The Gompertz model

PartiallyObservedMarkovProcesses.Examples.parus_dataConstant
parus_data

Parus major (Great Tit) census (all individuals), Wytham Wood, Oxfordshire. Global Population Dynamics Database dataset #10163. (NERC Centre for Population Biology, Imperial College (2010) The Global Population Dynamics Database Version 2. http://www.sw.ic.ac.uk/cpb/cpb/gpdd.html).

Original source: McCleery, R. & Perrins, C. (1991) Effects of predation on the numbers of Great Tits, Parus major. In: Bird Population Studies, edited by Perrins, C.M., Lebreton, J.-D. & Hirons, G.J.M. Oxford. Univ. Press. pp. 129–147.

source
PartiallyObservedMarkovProcesses.Examples.gompertzFunction
gompertz()

gompertz is a PompObject containing Parus major data and a simple Gompertz population model. The population model has a single scalar state variable, $X_t$, which obeys

\[X_t = X_{t-1}^S\,K^{1-S}\,\varepsilon_t,\]

where $S = e^{-r\delta{t}}$ and $\varepsilon_t \sim \mathrm{LogNormal}(0,\sigma_p)$. The time-step is one unit: $\delta{t}=1$. The data are assumed to be drawn from a log-normal distribution. In particular,

\[\mathrm{pop}_t \sim \mathrm{LogNormal}(\log{X_t},\sigma_m).\]

Parameters

  • r: the growth rate
  • K: the equilibrium population density
  • X0: the initial population density
  • σₚ: process noise s.d.
  • σₘ: measurement noise s.d.
source

View the Parus data:

using PartiallyObservedMarkovProcesses, RCall
using PartiallyObservedMarkovProcesses.Examples
P = gompertz()
d = melt(P)
R"""
$d |>
  pivot_longer(-year) |>
  ggplot(aes(x=year,y=value))+
  geom_line()+
  facet_wrap(~name,scales="free_y",ncol=1)+
  labs(y="")+
  theme_bw() -> pl
print(pl)
"""

gompertz_data

View a few representative simulations:

using PartiallyObservedMarkovProcesses, RCall
using PartiallyObservedMarkovProcesses.Examples
P = gompertz()
Q = simulate(P;params=(r=4.5,K=210.0,σₚ=0.7,σₘ=0.1,X0=150.0),nsim=5)
d = melt(Q,:parset,:rep)
R"""
$d |>
  pivot_longer(-c(year,rep,parset)) |>
  ggplot(aes(x=year,y=value,group=rep,color=factor(rep)))+
  geom_line()+
  geom_point()+
  facet_wrap(~name,scales="free_y",ncol=1)+
  labs(y="",color="replicate")+
  theme_bw() -> pl
print(pl)
"""

gompertz_sims

A simple SIR model

PartiallyObservedMarkovProcesses.Examples.sirFunction
sir(
    β = 0.5, γ = 0.25, N = 10000,
    ρ = 0.3, k = 10,
    S0 = 0.9, I0 = 0.01, R0 = 0.1,
    δt = 0.1, t0 = 0.0,
    times = range(start=1.0,stop=90,step=1.0)
   )

sir returns a PompObject containing simulated SIR data.

Parameters

  • β: transmission rate
  • γ: recovery rate
  • N: population size
  • ρ: reporting rate
  • k: overdispersion coefficient (negative binomial size parameter)
  • S0, I0, R0: relative proportions of susceptible, infected, recovered (respectively) in the population at t=t0.
  • δt: Euler stepsize
  • t0: zero-time
  • times: vector of observation times
source

sir diagram

The Rosenzweig-MacArthur model

Deterministic version

PartiallyObservedMarkovProcesses.Examples.drmcaFunction
drmca(
     r = 1, K = 1e4, A = 1e3,
     b = 1e-3, c = 1, m = 0.8,
     σ = 0.01,
     N0 = 3000, P0 = 4, t0 = 0.0,
     times = range(start=0,stop=500,step=0.2),
     integrator = AutoTsit5(Rosenbrock23())
   )

Parameters

  • r: intrinsic growth rate of prey
  • K: carrying capacity for prey
  • A: half-saturation prey density
  • c: predator foraging rate
  • b: predator yield (predators born per prey item killed)
  • m: predator death rate
  • σ: measurement noise magnitude
  • N0, P0: initial densities
  • t0: zero-time
  • times: vector of observation times
  • integrator: integration algorithm.

Observables

  • n: prey density
  • p: predator density

State variables

  • $X = \log(N)$
  • $Y = \log(P)$

Details

drmca returns a PompObject containing simulated data from a deterministic Rosenzweig-MacArthur model implemented as a vectorfield. Specifically, the model is the classical Rosenzweig-MacArthur model

\[\begin{aligned} \frac{dN}{dt} &= r N \left(1-\frac{N}{K}\right) - \frac{c N P}{1+N/A} \\ \frac{dP}{dt} &= \frac{b c N P}{1+N/A} - m P \\ \end{aligned}\]

There is also a stochastic version of this model available as an example (see rmca).

In this system, the predator is inviable unless $R = \frac{bcA}{m} > 1$. Even if the predator is viable, the environment is too impoverished to support predators unless $R>1+\frac{A}{K}$. If the environment is rich enough, and if moreover $R>\frac{1+\frac{A}{K}}{1-\frac{A}{K}}$, then the nontrivial equilibrium of the system is unstable. For the default parameters, we have $R = 1.25$ and $\frac{A}{K} = 0.1$, so the latter condition holds.

source
P = drmca(σ=0.1,times=range(0,400.0,step=1.0))
d = melt(P)
R"""
$d |>
  pivot_longer(-time) |>
  ggplot(aes(x=time,y=value))+
  geom_path()+
  facet_wrap(~name,scales="free_y",ncol=1)+
  labs(y="")+
  theme_bw() -> pl
print(pl)
"""

drmca_dynamics

Stochastic version

PartiallyObservedMarkovProcesses.Examples.rmcaFunction
rmca(
     r = 1, K = 1e4, A = 1e3,
     b = 1e-3, c = 1, m = 0.8,
     V = 100, σ = 0.01,
     N0 = 3000, P0 = 4, t0 = 0.0,
     δt = 0.01,
     times=range(start=0,stop=500,step=0.2)
   )

Parameters

  • r: intrinsic growth rate of prey
  • K: carrying capacity for prey
  • A: half-saturation prey density
  • c: predator foraging rate
  • b: predator yield (predators born per prey item killed)
  • m: predator death rate
  • V: system size
  • σ: measurement noise magnitude
  • N0, P0: initial densities
  • t0: zero-time
  • δt: Euler stepsize
  • times: vector of observation times

Observables

  • n: prey density
  • p: predator density

State variables

  • $X = \log(N)$
  • $Y = \log(P)$

Details

rmca returns a PompObject containing simulated data from a Rosenzweig-MacArthur model implemented as an Itô diffusion. Specifically, if $N$ and $P$ are prey and predator densities, respectively, then $dN = dG - dC - dS$ and $dP = b dS - dM$, where

\[\begin{aligned} dG &= r N dt + \sqrt{\frac{1}{V} r N} dW_1 \\ dC &= \frac{r N^2}{K} dt + \sqrt{\frac{1}{V} \frac{r N^2}{K}} dW_2 \\ dS &= \frac{c N P}{1+N/A} dt + \sqrt{\frac{1}{V} \frac{c N P}{1+N/A}} dW_3 \\ dM &= m P dt + \sqrt{\frac{1}{V} m P} dW_4 \\ \end{aligned}\]

Here, the $dW_i$ are increments of independent standard Wiener processes. Thus, the process noise scales demographically. Specifically, the system size, $V$, converts the densities $N$, $P$ into numbers. It controls the relative magnitude (coefficient of variation) of the demographic process noise. Moreover, $V$ determines a lower threshold on the population sizes, such that if ever $N V < 1$ or $P V < 1$, the population is taken to be extinct. Otherwise, it plays no role in the dynamics. The measurement error is assumed to scale environmentally:

\[\begin{aligned} n &\sim \mathrm{LogNormal}(\log{N},\sigma) \\ p &\sim \mathrm{LogNormal}(\log{P},\sigma) \\ \end{aligned}\]

Note that, in the limit $V\to\infty$, the Itô diffusion becomes the classical Rosenzweig-MacArthur ordinary differential equation model (see drmca):

\[\begin{aligned} \frac{dN}{dt} &= r N \left(1-\frac{N}{K}\right) - \frac{c N P}{1+N/A} \\ \frac{dP}{dt} &= \frac{b c N P}{1+N/A} - m P \\ \end{aligned}\]

source

A sample simulation.

P = rmca(σ=0.1,times=range(0,400.0,step=1.0))
d = melt(P)
R"""
$d |>
  pivot_longer(-time) |>
  ggplot(aes(x=time,y=value))+
  geom_path()+
  facet_wrap(~name,scales="free_y",ncol=1)+
  labs(y="")+
  theme_bw() -> pl
print(pl)
"""

rmca_dynamics

Multivariate Brownian motion

PartiallyObservedMarkovProcesses.Examples.brownian_motionFunction
brownian_motion(;times, t0 = 0, x0, σ, τ)

returns a PompObject for a multivariate Brownian motion process with observations at times. The zero-time is t0 and the starting state at that time is x0. The intensity of the Brownian motion is given by the matrix σ. The measurement error is multivariate normal with variance τ^⊤ τ.

source