Lotka-Volterra

analysis.py

# code adapted from https://github.com/stan-dev/example-models/blob/master/knitr/lotka-volterra/lotka-volterra-predator-prey.Rmd
# this is a close translation of the original R code

import micropip
# ggplot style graphics for python
await micropip.install("plotnine") 

import pandas as pd
import numpy as np
from plotnine import *

# Read the data
lynx_hare_df = pd.read_csv("https://raw.githubusercontent.com/stan-dev/example-models/refs/heads/master/knitr/lotka-volterra/hudson-bay-lynx-hare.csv", comment="#")

z_init_draws = draws.get("z_init") 
z_draws = draws.get("z")
y_init_rep_draws = draws.get("y_init_rep")
y_rep_draws = draws.get("y_rep")

predicted_pelts = np.zeros((21, 2))
min_pelts = np.zeros((21, 2))
max_pelts = np.zeros((21, 2))

for k in range(2):
    predicted_pelts[0, k] = np.mean(y_init_rep_draws[:, :, k])
    min_pelts[0, k] = np.quantile(y_init_rep_draws[:, :, k], 0.25)
    max_pelts[0, k] = np.quantile(y_init_rep_draws[:, :, k], 0.75)
    
    for n in range(1, 21):
        predicted_pelts[n, k] = np.mean(y_rep_draws[:, :, n-1, k])
        min_pelts[n, k] = np.quantile(y_rep_draws[:, :, n-1, k], 0.25)
        max_pelts[n, k] = np.quantile(y_rep_draws[:, :, n-1, k], 0.75)

# Melt the dataframe - similar to R's melt from reshape2
lynx_data = pd.DataFrame({'year': range(len(lynx_hare_df)), 'species': 'Lynx', 'pelts': lynx_hare_df.iloc[:, 2]})
hare_data = pd.DataFrame({'year': range(len(lynx_hare_df)), 'species': 'Hare', 'pelts': lynx_hare_df.iloc[:, 1]})
lynx_hare_melted_df = pd.concat([lynx_data, hare_data])

# Add 1899 to year
lynx_hare_melted_df['year'] = lynx_hare_melted_df['year'] + 1899

# Create observe dataframe
Nmelt = len(lynx_hare_melted_df)
lynx_hare_observe_df = lynx_hare_melted_df.copy()
lynx_hare_observe_df['source'] = 'measurement'

# Create prediction dataframe
years = list(range(1900, 1921)) * 2
species = ['Lynx'] * 21 + ['Hare'] * 21
pelts = np.concatenate([predicted_pelts[:, 1], predicted_pelts[:, 0]])
min_pelts_flat = np.concatenate([min_pelts[:, 1], min_pelts[:, 0]])
max_pelts_flat = np.concatenate([max_pelts[:, 1], max_pelts[:, 0]])

lynx_hare_predict_df = pd.DataFrame({
    'year': years,
    'species': species,
    'pelts': pelts,
    'min_pelts': min_pelts_flat,
    'max_pelts': max_pelts_flat,
    'source': 'prediction'
})

# Combine the dataframes
lynx_hare_observe_df['min_pelts'] = np.nan
lynx_hare_observe_df['max_pelts'] = np.nan
lynx_hare_observe_predict_df = pd.concat([lynx_hare_observe_df, lynx_hare_predict_df])

# First plot
p1 = (ggplot(lynx_hare_observe_predict_df, 
            aes(x='year', y='pelts', color='source'))
     + geom_vline(xintercept=1900, color='grey')
     + geom_hline(yintercept=0, color='grey')
     + facet_wrap('~ species', ncol=1)
     + geom_ribbon(aes(ymin='min_pelts', ymax='max_pelts'), alpha=0.1)
     + geom_line()
     + geom_point()
     + ylab("pelts (thousands)")
     + theme(legend_position="bottom")
     + labs(caption="Posterior predictive checks, including posterior means and 50% intervals along\nwith the measured data.\nIf the model is well calibrated, as this one appears to be, 50% of the points are expected\nto fall in their50% intervals."))

p1.show()

# Second plot - orbit from z_draws
data_list = []
for m in range(100):  # First 100 draws
    for i in range(12):  # First 12 time points
        data_list.append({
            'lynx': z_draws[0,m, i, 0],
            'hare': z_draws[0,m, i, 1],
            'draw': m + 1
        })
        
df = pd.DataFrame(data_list)

p2 = (ggplot(df, aes(x='lynx', y='hare', color='factor(draw)'))
     + geom_vline(xintercept=0, color='grey')
     + geom_hline(yintercept=0, color='grey')
     + geom_path(size=0.75, alpha=0.2)
     + xlab("lynx pelts (thousands)")
     + ylab("hare pelts (thousands)")
     + theme(legend_position="none")
     + labs(caption="Plot of expected population orbit for one hundred draws from the posterior.\nEach draw represents a different orbit determined by the differential equation system parameters.\nTogether they provide a sketch of posterior uncertainty for the *expected* population dynamics.\nIf the ODE solutions were extracted per month rather than per year,\nthe resulting plots would appear fairly smooth."))

p2.show()

# Third plot - orbit from y_rep_draws
data_list = []
for m in range(100):  # First 100 draws
    for i in range(12):  # First 12 time points
        data_list.append({
            'lynx': y_rep_draws[0, m, i, 0], 
            'hare': y_rep_draws[0, m, i, 1],
            'draw': m + 1
        })

df = pd.DataFrame(data_list)

p3 = (ggplot(df, aes(x='lynx', y='hare', color='factor(draw)'))
     + geom_vline(xintercept=0, color='grey')
     + geom_hline(yintercept=0, color='grey')
     + geom_path(size=0.75, alpha=0.2)
     + xlab("lynx pelts (thousands)")
     + ylab("hare pelts (thousands)")
     + theme(legend_position="none")
     + labs(caption="Plot of expected pelt collection orbits for one hundred draws of system parameters from the posterior.\nEven if plotted at more fine-grained time intervals, error would remove any apparent smoothness.\nExtreme draws as seen here are typical when large values have high error on the multiplicative scale."))

p3.show()

analysis.R

# code adapted from https://github.com/stan-dev/example-models/blob/master/knitr/lotka-volterra/lotka-volterra-predator-prey.Rmd

library(posterior)
library(reshape2)
library(ggplot2)

lynx_hare_df <-
  read.csv("https://raw.githubusercontent.com/stan-dev/example-models/refs/heads/master/knitr/lotka-volterra/hudson-bay-lynx-hare.csv",
           comment.char="#")

z_init_draws <- extract_variable_array(draws, "z_init")
z_draws <- extract_variable_array(draws, "z")
y_init_rep_draws <- extract_variable_array(draws, "y_init_rep")
y_rep_draws <- extract_variable_array(draws, "y_rep")
predicted_pelts <- matrix(NA, 21, 2)
min_pelts <- matrix(NA, 21, 2)
max_pelts <- matrix(NA, 21, 2)


for (k in 1:2) {
  predicted_pelts[1, k] <- mean(y_init_rep_draws[ ,, k])
  min_pelts[1, k] <- quantile(y_init_rep_draws[ ,, k], 0.25)
  max_pelts[1, k] <- quantile(y_init_rep_draws[ ,, k], 0.75)
  for (n in 2:21) {
    predicted_pelts[n, k] <- mean(y_rep_draws[ ,, n - 1, k])
    min_pelts[n, k] <- quantile(y_rep_draws[ ,, n - 1, k], 0.25)
    max_pelts[n, k] <- quantile(y_rep_draws[ ,, n - 1, k], 0.75)
  }
}

lynx_hare_melted_df <- melt(as.matrix(lynx_hare_df[, 2:3]))
colnames(lynx_hare_melted_df) <- c("year", "species", "pelts")
lynx_hare_melted_df$year <-
  lynx_hare_melted_df$year +
  rep(1899, length(lynx_hare_melted_df$year))

Nmelt <- dim(lynx_hare_melted_df)[1]
lynx_hare_observe_df <- lynx_hare_melted_df
lynx_hare_observe_df$source <- rep("measurement", Nmelt)

lynx_hare_predict_df <-
  data.frame(year = rep(1900:1920, 2),
             species = c(rep("Lynx", 21), rep("Hare", 21)),
             pelts = c(predicted_pelts[, 2],
                       predicted_pelts[, 1]),
             min_pelts = c(min_pelts[, 2], min_pelts[, 1]),
             max_pelts = c(max_pelts[, 2], max_pelts[, 1]),
             source = rep("prediction", 42))

lynx_hare_observe_df$min_pelts = lynx_hare_predict_df$min_pelts
lynx_hare_observe_df$max_pelts = lynx_hare_predict_df$max_pelts

lynx_hare_observe_predict_df <-
  rbind(lynx_hare_observe_df, lynx_hare_predict_df)

ggplot(data = lynx_hare_observe_predict_df,
         aes(x = year, y = pelts, color = source)) +
  geom_vline(xintercept = 1900, color = "grey") +
  geom_hline(yintercept = 0, color = "grey") +
  facet_wrap( ~ species, ncol = 1) +
  geom_ribbon(aes(ymin = min_pelts, ymax = max_pelts),
	      colour = NA, fill = "black", alpha = 0.1) +
  geom_line() +
  geom_point() +
  ylab("pelts (thousands)") +
  theme(legend.position="bottom") +
  labs(caption="Posterior predictive checks, including posterior means and 50% intervals along with the measured data.\nIf the model is well calibrated, as this one appears to be, 50% of the points are expected to fall in their 50% intervals.")

# ---------------------------

df <- data.frame(list(lynx = z_draws[1,1, 1:12 , 1], hare = z_draws[1,1, 1:12, 2], draw = 1))
for (m in 2:100) {
  df <- rbind(df, data.frame(list(lynx = z_draws[m,1, 1:12 , 1], hare = z_draws[m,1, 1:12, 2], draw = m)))
}
ggplot(df) +
  geom_vline(xintercept = 0, color = "grey") +
  geom_hline(yintercept = 0, color = "grey") +
  geom_path(aes(x = lynx, y = hare, colour = draw), linewidth = 0.75, alpha = 0.2) +
  xlab("lynx pelts (thousands)") +
  ylab("hare pelts (thousands)") +
  theme(legend.position="none") +
  labs(caption="Plot of expected population orbit for one hundred draws from the posterior.\nEach draw represents a different orbit determined by the differential equation system parameters.\nTogether they provide a sketch of posterior uncertainty for the *expected* population dynamics.\nIf the ODE solutions were extracted per month rather than per year, the resulting plots would appear fairly smooth.")

# ---------------------------

df <- data.frame(list(lynx = y_rep_draws[1,1, 1:12 , 1], hare = y_rep_draws[1,1, 1:12, 2], draw = 1))
for (m in 2:100) {
  df <- rbind(df, data.frame(list(lynx = y_rep_draws[m,1, 1:12 , 1], hare = y_rep_draws[m,1, 1:12, 2], draw = m)))
}

ggplot(df) +
  geom_vline(xintercept = 0, color = "grey") +
  geom_hline(yintercept = 0, color = "grey") +
  geom_path(aes(x = lynx, y = hare, colour = draw), linewidth = 0.75, alpha = 0.2) +
  xlab("lynx pelts (thousands)") +
  ylab("hare pelts (thousands)") +
  theme(legend.position="none") +
  labs(caption="Plot of expected pelt collection orbits for one hundred draws of system parameters from the posterior.\nEven if plotted at more fine-grained time intervals, error would remove any apparent smoothness.\nExtreme draws as seen here are typical when large values have high error on the multiplicative scale.")

data.json

{
  "N": 20,
  "ts": [
    1,
    2,
    3,
    4,
    5,
    6,
    7,
    8,
    9,
    10,
    11,
    12,
    13,
    14,
    15,
    16,
    17,
    18,
    19,
    20
  ],
  "y_init": [
    30,
    4
  ],
  "y": [
    [
      47.2,
      6.1
    ],
    [
      70.2,
      9.8
    ],
    [
      77.4,
      35.2
    ],
    [
      36.3,
      59.4
    ],
    [
      20.6,
      41.7
    ],
    [
      18.1,
      19
    ],
    [
      21.4,
      13
    ],
    [
      22,
      8.3
    ],
    [
      25.4,
      9.1
    ],
    [
      27.1,
      7.4
    ],
    [
      40.3,
      8
    ],
    [
      57,
      12.3
    ],
    [
      76.6,
      19.5
    ],
    [
      52.3,
      45.7
    ],
    [
      19.5,
      51.1
    ],
    [
      11.2,
      29.7
    ],
    [
      7.6,
      15.8
    ],
    [
      14.6,
      9.7
    ],
    [
      16.2,
      10.1
    ],
    [
      24.7,
      8.6
    ]
  ]
}

data.py

import pandas as pd
import numpy as np

# pandas can directly download, avoiding the need for `requests` etc
lynx_hare_df = pd.read_csv("https://raw.githubusercontent.com/stan-dev/example-models/refs/heads/master/knitr/lotka-volterra/hudson-bay-lynx-hare.csv", 
                           comment="#")

print(lynx_hare_df.head())

N = len(lynx_hare_df) - 1
ts = np.arange(1, N+1)
y_init = [lynx_hare_df.iloc[0, 1], lynx_hare_df.iloc[0, 2]]
y = lynx_hare_df.iloc[1:N+1, 1:3].values

# hare, lynx order (swapping columns)
y = np.column_stack((y[:, 1], y[:, 0]))

data = {
    'N': N,
    'ts': ts,
    'y_init': y_init,
    'y': y
}

data.R

lynx_hare_df <-
  read.csv("https://raw.githubusercontent.com/stan-dev/example-models/refs/heads/master/knitr/lotka-volterra/hudson-bay-lynx-hare.csv",
           comment.char="#")
head(lynx_hare_df)
N <- length(lynx_hare_df$Year) - 1
ts <- 1:N
y_init <- c(lynx_hare_df$Hare[1], lynx_hare_df$Lynx[1])
y <- as.matrix(lynx_hare_df[2:(N + 1), 2:3])
y <- cbind(y[ , 2], y[ , 1]); # hare, lynx order

data <- list(N = N, ts = ts, y_init = y_init, y = y)

main.stan

// Read the accompanying case study:
// https://mc-stan.org/learn-stan/case-studies/lotka-volterra-predator-prey.html
functions {
  vector dz_dt(real t, // time
               vector z, real alpha, real beta, real gamma, real delta) {
    real u = z[1];
    real v = z[2];
    
    real du_dt = (alpha - beta * v) * u;
    real dv_dt = (-gamma + delta * u) * v;
    
    return [du_dt, dv_dt]';
  }
}
data {
  int<lower=0> N; // number of measurement times
  array[N] real ts; // measurement times > 0
  array[2] real y_init; // initial measured populations
  array[N, 2] real<lower=0> y; // measured populations
}
transformed data {
  real rel_tol = 1e-5;
  real abs_tol = 1e-3;
  int max_num_steps = to_int(5e2);
}
parameters {
  real<lower=0> alpha, beta, gamma, delta;
  vector<lower=0>[2] z_init; // initial population
  array[2] real<lower=0> sigma; // measurement errors
}
transformed parameters {
  array[N] vector[2] z = ode_rk45_tol(dz_dt, z_init, 0.0, ts, rel_tol,
                                      abs_tol, max_num_steps, alpha, beta,
                                      gamma, delta);
}
model {
  alpha ~ normal(1, 0.5);
  gamma ~ normal(1, 0.5);
  beta ~ normal(0.05, 0.05);
  delta ~ normal(0.05, 0.05);
  sigma ~ lognormal(-1, 1);
  z_init ~ lognormal(log(10), 1);
  for (k in 1 : 2) {
    y_init[k] ~ lognormal(log(z_init[k]), sigma[k]);
    y[ : , k] ~ lognormal(log(z[ : , k]), sigma[k]);
  }
}
generated quantities {
  array[2] real y_init_rep;
  array[N, 2] real y_rep;
  for (k in 1 : 2) {
    y_init_rep[k] = lognormal_rng(log(z_init[k]), sigma[k]);
    for (n in 1 : N) {
      y_rep[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
    }
  }
}

meta.json

{"title":"Lotka-Volterra","dataSource":"generated_by_r"}

sampling_opts.json

{"num_chains":4,"num_warmup":1000,"num_samples":1000,"init_radius":2}