In a previous post I used ellipses to put hard bounds on where a moose could have been between two GPS fixes, assuming only a maximum possible speed. That approach answers “where could it not have been” and never says anything about likelihood. This post is about the complementary probabilistic tool, and about the single most widely used stochastic model of animal movement: the Ornstein-Uhlenbeck (OU) process.
The OU process is a diffusion with a linear restoring force. Interpreted for moose, it says that the animal is pulled back towards a home-range centre and buffeted by random motion, and it is the special case of the general Fokker-Planck framework in which the density stays Gaussian and every quantity of interest is available in closed form. That tractability is exactly why it is a natural first model, and also why it is too restrictive to be the last one. In the next post I lift the restrictions one at a time and build a general Fokker-Planck model of moose movement; this post establishes the baseline that the general model has to beat.
From paths to densities
Start with the individual. An It^o stochastic differential equation (SDE) says that a small change in position over a small change in time has a predictable part and a random part:
\[dX_t = \mu(X_t, t)\,dt + \sigma(X_t, t)\,dW_t,\]
where \(\mu\) is the drift vector, \(\sigma\) the diffusion matrix, and \(W_t\) a Wiener process (Brownian motion). Each sample path of this SDE is a plausible trajectory for one animal: drift is the average or directed motion, diffusion is the unpredictable jitter.
The density \(p(x, t)\) of where the animal is likely to be then obeys the Fokker-Planck equation (also called the Kolmogorov forward equation):
\[\frac{\partial p}{\partial t} = -\nabla \cdot \bigl(\mu\, p\bigr) + \frac{1}{2}\,\partial_i \partial_j \bigl(D_{ij}\, p\bigr), \qquad D = \sigma\sigma^{\top}.\]
The SDE and the Fokker-Planck equation are two descriptions of the same stochastic process: the SDE is the language of paths, the PDE the language of densities. The OU process is what you get when the drift is linear and the diffusion is constant and isotropic.
The Ornstein-Uhlenbeck model
Moose are not diffusing gas molecules, but they do have a strong tendency to stay within a home range. The OU model encodes that with a drift that pulls the animal back towards a centre \(c\),
\[\mu(x) = -\kappa\,(x - c),\]
so that far from the centre the drift is large and inward, and near the centre it is negligible. The two parameters have clear interpretations:
- \(\kappa > 0\) is the fidelity to the home range, in units of \(1/\text{time}\). Its reciprocal \(1/\kappa\) is the characteristic time over which the animal forgets its current position and relaxes back towards the centre. It is the relaxation or memory time.
- \(\sigma\) is the restlessness, in units of \(\text{distance}/\sqrt{\text{time}}\). It is not a speed: it is the strength of the random component, and its units are those of a diffusion coefficient’s square root.
Because the drift is linear and the noise is constant, the model is exactly solvable. The stationary distribution of positions is Gaussian,
\[p_\infty(x) = \frac{1}{\sqrt{2\pi s^2}}\exp\!\left(-\frac{(x-c)^2}{2s^2}\right), \qquad s^2 = \frac{\sigma^2}{2\kappa},\]
so the home-range size \(s\) is set by the ratio of restlessness to fidelity. In two dimensions the same model gives an isotropic Gaussian around \(c\) with per-axis standard deviation \(s\), so the radius from the centre follows a Rayleigh distribution with mean \(s\sqrt{\pi/2}\). The full time-dependent solution is also known exactly: starting from a point at the centre, the position at time \(t\) is Gaussian with mean \(c\) and variance
\[s^2(t) = \frac{\sigma^2}{2\kappa}\left(1 - e^{-2\kappa t}\right),\]
which approaches \(s^2\) at rate \(\kappa\). This is the property that makes the OU model the reference point: from two numbers, \(\kappa\) and \(\sigma\), and a centre, you can write down the entire distribution at any time.
Simulating the paths
The Rust program below simulates a virtual population of moose as two-dimensional OU processes with an Euler-Maruyama scheme, which is the standard first-order discretization of an SDE. It depends only on the rand and rand_distr crates. Parameters are in kilometres and hours: a fidelity of \(\kappa = 0.5\ \text{h}^{-1}\) (a two-hour memory), a restlessness of \(\sigma = 1.0\ \text{km}/\sqrt{\text{h}}\), and a time step of \(36\) seconds over twenty hours. Every moose starts at the home-range centre.
//! Virtual moose movement from an Ornstein-Uhlenbeck process.
//!
//! Units are kilometres and hours. Each virtual moose has a home-range
//! centre and is pulled back towards it with strength `kappa`, while an
//! exploratory random walk of strength `sigma` pushes it away:
//!
//! ```text
//! dx = -kappa * (x - c) dt + sigma * dW
//! ```
//!
//! The collection of moose is simulated as independent It^o processes with
//! an Euler-Maruyama scheme. The density of the positions then evolves
//! according to the Fokker-Planck equation
//!
//! ```text
//! dp/dt = kappa * d/dx[(x - c) p] + (sigma^2 / 2) * d^2p/dx^2
//! ```
//!
//! which for this drift and diffusion has the closed-form stationary
//! solution `Normal(c, sigma / sqrt(2 * kappa))` per axis. The printed
//! summary compares the simulated final positions against that analytic
//! answer.
use std::fs::File;
use std::io::Write;
use rand::SeedableRng;
use rand::rngs::StdRng;
use rand_distr::{Distribution, Normal};
/// Parameters for the movement model.
#[derive(Clone, Copy, Debug)]
pub struct MooseConfig {
/// Home-range fidelity, in 1/h.
pub kappa: f64,
/// Exploratory jitter, in km/sqrt(h).
pub sigma: f64,
/// Time step, in h.
pub dt: f64,
/// Number of time steps.
pub steps: usize,
/// Number of virtual moose.
pub moose: usize,
/// Home-range centre, in km.
pub centre: f64,
/// Seed for reproducibility.
pub seed: u64,
}
impl Default for MooseConfig {
fn default() -> Self {
Self {
kappa: 0.5, // Strong fidelity: relaxation time 1/kappa = 2 h
sigma: 1.0, // Exploratory motion
dt: 0.01, // 36 s fixes, finer than a typical GPS schedule
steps: 2_000, // 20 h of simulated movement
moose: 10_000, // Virtual population size
centre: 0.0, // Home-range centre at the origin
seed: 20_261_016, // Arbitrary fixed seed
}
}
}
impl MooseConfig {
/// Total simulated time, in hours.
#[must_use]
pub fn t_end(&self) -> f64 {
self.dt * self.steps as f64
}
/// Analytic stationary standard deviation per axis, in km.
#[must_use]
pub fn stationary_sd(&self) -> f64 {
self.sigma / (2.0 * self.kappa).sqrt()
}
/// Analytic mean distance from the centre at stationarity, in km.
#[must_use]
pub fn mean_radius(&self) -> f64 {
self.stationary_sd() * std::f64::consts::FRAC_PI_2.sqrt()
}
}
/// Simulate the final positions of a population of virtual moose.
///
/// Every moose starts at the home-range centre and is stepped forward with
/// the Euler-Maruyama scheme.
#[must_use]
pub fn simulate(config: &MooseConfig) -> Vec<(f64, f64)> {
let mut positions = vec![(config.centre, config.centre); config.moose];
let mut rng = StdRng::seed_from_u64(config.seed);
// Standard normal draws for the increments of the Wiener process.
let normal = Normal::new(0.0_f64, 1.0).unwrap();
let sqrt_dt = config.dt.sqrt();
// Euler-Maruyama time stepping of the It^o process.
for _ in 0..config.steps {
for position in positions.iter_mut() {
let dw_x = normal.sample(&mut rng) * sqrt_dt;
let dw_y = normal.sample(&mut rng) * sqrt_dt;
position.0 += -config.kappa * (position.0 - config.centre) * config.dt
+ config.sigma * dw_x;
position.1 += -config.kappa * (position.1 - config.centre) * config.dt
+ config.sigma * dw_y;
}
}
positions
}
/// Summary statistics of a simulated population, with analytic references.
#[derive(Clone, Copy, Debug)]
pub struct Summary {
/// Sample mean of the x coordinate, in km.
pub mean_x: f64,
/// Sample mean of the y coordinate, in km.
pub mean_y: f64,
/// Sample standard deviation of x, in km.
pub sd_x: f64,
/// Sample standard deviation of y, in km.
pub sd_y: f64,
/// Sample mean distance from the centre, in km.
pub mean_radius: f64,
/// Home-range centre, in km.
pub centre: f64,
/// Analytic stationary standard deviation per axis, in km.
pub stationary_sd: f64,
/// Analytic mean distance from the centre at stationarity, in km.
pub mean_radius_analytic: f64,
}
/// Compute summary statistics of a simulated population.
#[must_use]
pub fn summarise(config: &MooseConfig, positions: &[(f64, f64)]) -> Summary {
let n = positions.len() as f64;
let mean_x = positions.iter().map(|p| p.0).sum::<f64>() / n;
let mean_y = positions.iter().map(|p| p.1).sum::<f64>() / n;
let sd_x = (positions.iter().map(|p| (p.0 - mean_x).powi(2)).sum::<f64>() / n).sqrt();
let sd_y = (positions.iter().map(|p| (p.1 - mean_y).powi(2)).sum::<f64>() / n).sqrt();
let mean_radius = positions
.iter()
.map(|p| p.0.hypot(p.1))
.sum::<f64>()
/ n;
Summary {
mean_x,
mean_y,
sd_x,
sd_y,
mean_radius,
centre: config.centre,
stationary_sd: config.stationary_sd(),
mean_radius_analytic: config.mean_radius(),
}
}
/// Count how many positions fall in each of `bins` equal bins on `[min, max)`.
#[must_use]
pub fn histogram(positions: &[(f64, f64)], min: f64, max: f64, bins: usize) -> Vec<usize> {
let bin_size = (max - min) / bins as f64;
let mut counts = vec![0_usize; bins];
for position in positions {
if position.0 >= min && position.0 < max {
let bin = ((position.0 - min) / bin_size) as usize;
counts[bin] += 1;
}
}
counts
}
fn print_config(config: &MooseConfig) {
println!("Movement model: dx = -kappa (x - c) dt + sigma dW (Ornstein-Uhlenbeck)");
println!(
"kappa = {:.3} 1/h, sigma = {:.3} km/sqrt(h), dt = {:.4} h, steps = {} (T = {:.1} h)",
config.kappa,
config.sigma,
config.dt,
config.steps,
config.t_end()
);
println!(
"{} virtual moose, seed {}, home-range centre at {:.1} km",
config.moose, config.seed, config.centre
);
println!(
"Analytic stationary sd per axis: sigma / sqrt(2 kappa) = {:.3} km",
config.stationary_sd()
);
}
fn print_summary(summary: &Summary) {
println!("\nFinal positions at t = T:");
println!("{:<28}{:>10}{:>10}", "quantity", "simulated", "analytic");
println!(
"{:<28}{:>10.3}{:>10.3}",
"mean x (km)", summary.mean_x, summary.centre
);
println!(
"{:<28}{:>10.3}{:>10.3}",
"mean y (km)", summary.mean_y, summary.centre
);
println!(
"{:<28}{:>10.3}{:>10.3}",
"sd x (km)", summary.sd_x, summary.stationary_sd
);
println!(
"{:<28}{:>10.3}{:>10.3}",
"sd y (km)", summary.sd_y, summary.stationary_sd
);
println!(
"{:<28}{:>10.3}{:>10.3}",
"mean radius from centre (km)", summary.mean_radius, summary.mean_radius_analytic
);
}
fn print_histogram(positions: &[(f64, f64)]) {
let bins = 12;
let min = -3.0;
let max = 3.0;
let bin_size = (max - min) / bins as f64;
let counts = histogram(positions, min, max, bins);
println!("\nMarginal distribution of x (each * is about 40 moose):");
for (i, &count) in counts.iter().enumerate() {
let bin_start = min + i as f64 * bin_size;
let bar = "*".repeat(count / 40);
println!(
"{:>5.2} to {:>5.2} | {} ({})",
bin_start,
bin_start + bin_size,
bar,
count
);
}
}
/// Final positions of every virtual moose, for downstream plotting.
fn write_positions(positions: &[(f64, f64)]) -> std::io::Result<()> {
let mut file = File::create("data/moose_positions.csv")?;
writeln!(file, "x,y")?;
for (x, y) in positions {
writeln!(file, "{x:.6},{y:.6}")?;
}
Ok(())
}
/// One representative trajectory, subsampled for downstream plotting.
fn write_track(config: &MooseConfig) -> std::io::Result<()> {
let mut rng = StdRng::seed_from_u64(config.seed);
let normal = Normal::new(0.0_f64, 1.0).unwrap();
let sqrt_dt = config.dt.sqrt();
let mut x = config.centre;
let mut y = config.centre;
let mut file = File::create("data/moose_track.csv")?;
writeln!(file, "t,x,y")?;
let record_every = 10;
for step in 0..=config.steps {
if step % record_every == 0 {
writeln!(file, "{:.3},{x:.6},{y:.6}", step as f64 * config.dt)?;
}
if step == config.steps {
break;
}
let dw_x = normal.sample(&mut rng) * sqrt_dt;
let dw_y = normal.sample(&mut rng) * sqrt_dt;
x += -config.kappa * (x - config.centre) * config.dt + config.sigma * dw_x;
y += -config.kappa * (y - config.centre) * config.dt + config.sigma * dw_y;
}
Ok(())
}
/// Simulate, print the summary, and write the plotting CSVs for `config`.
///
/// # Errors
///
/// Returns an error if the CSV files under `data/` cannot be written.
pub fn run_with_config(config: &MooseConfig) -> std::io::Result<()> {
let positions = simulate(config);
let summary = summarise(config, &positions);
print_config(config);
print_summary(&summary);
print_histogram(&positions);
write_positions(&positions)?;
write_track(config)?;
Ok(())
}
/// Run the model with the default parameters used in the blog post.
///
/// # Errors
///
/// Returns an error if the CSV files under `data/` cannot be written.
pub fn run() -> std::io::Result<()> {
run_with_config(&MooseConfig::default())
}Movement model: dx = -kappa (x - c) dt + sigma dW (Ornstein-Uhlenbeck)
kappa = 0.500 1/h, sigma = 1.000 km/sqrt(h), dt = 0.0100 h, steps = 2000 (T = 20.0 h)
10000 virtual moose, seed 20261016, home-range centre at 0.0 km
Analytic stationary sd per axis: sigma / sqrt(2 kappa) = 1.000 km
Final positions at t = T:
quantity simulated analytic
mean x (km) -0.003 0.000
mean y (km) -0.002 0.000
sd x (km) 1.003 1.000
sd y (km) 1.008 1.000
mean radius from centre (km) 1.261 1.253
Marginal distribution of x (each * is about 40 moose):
-3.00 to -2.50 | * (40)
-2.50 to -2.00 | **** (173)
-2.00 to -1.50 | ************ (480)
-1.50 to -1.00 | *********************** (930)
-1.00 to -0.50 | ************************************ (1446)
-0.50 to 0.00 | *********************************************** (1898)
0.00 to 0.50 | ************************************************* (1979)
0.50 to 1.00 | ************************************ (1449)
1.00 to 1.50 | ********************** (908)
1.50 to 2.00 | *********** (450)
2.00 to 2.50 | **** (169)
2.50 to 3.00 | * (45)Even at this coarse resolution the simulated moments match the analytic stationary values to within about a percent. This is a numerical check of the theory: the Wiener process was the only randomness in the program, and out of it came the Gaussian equilibrium predicted by the analytics. Notice, though, that the program integrates the SDE, not the density equation. The comment in the original scratch version of this code called the Euler-Maruyama update a “Fokker-Planck update step”, which is a common slip: the particles carry the density, but the equation they satisfy collectively is a different one.
The density equation for the OU model
To get the density equation itself onto the page, the Python script below solves the one-dimensional Fokker-Planck equation for the OU model on a grid with an explicit finite-volume scheme. The probability current is
\[J = -\kappa\,(x - c)\,p - \frac{\sigma^2}{2}\frac{\partial p}{\partial x},\]
and it is discretized in flux form with the advective part upwinded and the diffusive part centred. The flux is set to zero at the walls, so total probability is conserved exactly, and the initial condition is a point mass at the centre. The script then compares three things: the simulated moose from the Rust program, the numerical PDE solution, and the analytic Gaussian.
"""Fokker-Planck and analytic solutions of the Ornstein-Uhlenbeck moose model.
The companion Rust crate ``moose_ou`` simulates virtual moose whose positions
follow an Ornstein-Uhlenbeck (OU) process, that is, the It^o SDE
dx = -kappa (x - c) dt + sigma dW
The density of those positions, p(x, t), obeys the Fokker-Planck equation
dp/dt = kappa d/dx[(x - c) p] + (sigma^2 / 2) d^2p/dx^2.
This module solves that PDE directly on a grid with an explicit finite-volume
scheme (upwinded drift, centred diffusion, no-flux walls), and plots the PDE
solution against the simulated moose and against the analytic OU transition
density,
X_t ~ Normal(c, (sigma^2 / 2 kappa) (1 - exp(-2 kappa t))).
Run the Rust crate first so that ``data/moose_positions.csv`` and
``data/moose_track.csv`` exist.
"""
from __future__ import annotations
import math
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
KAPPA = 0.5
"""Home-range fidelity, in 1/h."""
SIGMA = 1.0
"""Exploratory jitter, in km/sqrt(h)."""
T_END = 20.0
"""Simulation length, in h."""
CENTRE = 0.0
"""Home-range centre, in km."""
def load_positions(path: str | Path) -> tuple[np.ndarray, np.ndarray]:
"""Load the final ``(x, y)`` positions of every moose from ``path``."""
data = np.loadtxt(path, delimiter=",", skiprows=1)
return data[:, 0], data[:, 1]
def load_track(path: str | Path) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Load a representative ``(t, x, y)`` trajectory from ``path``."""
data = np.loadtxt(path, delimiter=",", skiprows=1)
return data[:, 0], data[:, 1], data[:, 2]
def analytic_variance(t: float, *, kappa: float = KAPPA, sigma: float = SIGMA) -> float:
"""Variance of the OU process started at the centre.
The variance relaxes from zero towards ``sigma**2 / (2 kappa)`` like
``1 - exp(-2 kappa t)``.
"""
return sigma**2 / (2 * kappa) * (1 - np.exp(-2 * kappa * t))
def stationary_sd(*, kappa: float = KAPPA, sigma: float = SIGMA) -> float:
"""Stationary standard deviation per axis, ``sigma / sqrt(2 kappa)``."""
return sigma / math.sqrt(2 * kappa)
def analytic_pdf(
x: np.ndarray,
t: float,
*,
kappa: float = KAPPA,
sigma: float = SIGMA,
centre: float = CENTRE,
) -> np.ndarray:
"""Transition density of the OU process at time ``t``."""
variance = analytic_variance(t, kappa=kappa, sigma=sigma)
return np.exp(-0.5 * (x - centre) ** 2 / variance) / math.sqrt(
2 * math.pi * variance
)
def solve_fokker_planck(
extent: float = 5.0,
dx: float = 0.01,
dt: float = 5e-5,
*,
kappa: float = KAPPA,
sigma: float = SIGMA,
t_end: float = T_END,
centre: float = CENTRE,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Finite-volume solution of the 1D Fokker-Planck equation.
The probability current is written in flux form,
J = -kappa (x - c) p - (sigma^2 / 2) dp/dx,
with the advective part upwinded and the diffusive part centred, and the
time step is explicit. The flux is set to zero at the two walls, so mass
is conserved exactly. Upwinding adds a small numerical diffusion of size
``|v| dx / 2``, which is the main discretization error.
Returns ``(x, times, snapshots, masses, means, variances)``.
"""
x = np.arange(-extent, extent + dx / 2, dx)
n_steps = round(t_end / dt)
record_every = max(1, n_steps // 400)
p = np.zeros_like(x)
p[np.argmin(np.abs(x - centre))] = 1.0 / dx # Start as a point mass
diffusivity = sigma**2 / 2
def step(density: np.ndarray) -> np.ndarray:
# Probability current at the interior cell faces.
velocity_face = -kappa * (0.5 * (x[:-1] + x[1:]) - centre)
upwind = np.where(velocity_face > 0, density[:-1], density[1:])
advective = velocity_face * upwind
diffusive = -diffusivity * (density[1:] - density[:-1]) / dx
# No-flux walls, so the interior update telescopes to zero.
current = np.concatenate(([0.0], advective + diffusive, [0.0]))
return density - dt / dx * (current[1:] - current[:-1])
times = [0.0]
snapshots = [p.copy()]
masses = [float(np.sum(p) * dx)]
means = [float(np.sum(x * p) * dx)]
variances = [float(np.sum((x - centre) ** 2 * p) * dx)]
for step_index in range(n_steps):
p = step(p)
if (step_index + 1) % record_every == 0:
times.append((step_index + 1) * dt)
snapshots.append(p.copy())
masses.append(float(np.sum(p) * dx))
means.append(float(np.sum(x * p) * dx))
variances.append(float(np.sum((x - centre) ** 2 * p) * dx))
return (
x,
np.array(times),
np.array(snapshots),
np.array(masses),
np.array(means),
np.array(variances),
)
def diagnostics(
masses: np.ndarray,
means: np.ndarray,
variances: np.ndarray,
*,
kappa: float = KAPPA,
sigma: float = SIGMA,
t_end: float = T_END,
centre: float = CENTRE,
) -> dict[str, float]:
"""Compare the PDE solution at ``t_end`` with the analytic values."""
return {
"mass": float(masses[-1]),
"pde_variance": float(variances[-1]),
"analytic_variance": float(analytic_variance(t_end, kappa=kappa, sigma=sigma)),
"pde_mean": float(means[-1]),
"analytic_mean": float(centre),
}
def print_diagnostics(diag: dict[str, float]) -> None:
"""Print the PDE-versus-analytic comparison quoted in the blog post."""
print("Fokker-Planck PDE versus analytic solution:")
print(f" {'mass at T:':<20}{diag['mass']:.6f} (should be 1.000000)")
print(f" {'PDE variance at T:':<20}{diag['pde_variance']:.6f}")
print(f" {'analytic variance:':<20}{diag['analytic_variance']:.6f}")
print(f" {'PDE mean at T:':<20}{diag['pde_mean']:.6f}")
print(f" {'analytic mean:':<20}{diag['analytic_mean']:.6f}")
def plot_comparison(
x_positions: np.ndarray,
y_positions: np.ndarray,
t_track: np.ndarray,
x_track: np.ndarray,
y_track: np.ndarray,
solve: tuple[
np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray
],
output: str | Path,
*,
kappa: float = KAPPA,
sigma: float = SIGMA,
t_end: float = T_END,
centre: float = CENTRE,
) -> Path:
"""Render the four-panel comparison figure and return the written path."""
x, times, snapshots, _masses, _means, variances = solve
fig, axes = plt.subplots(2, 2, figsize=(13, 10))
# (a) One simulated moose trajectory.
ax = axes[0, 0]
points = ax.scatter(x_track, y_track, c=t_track, cmap="viridis", s=6)
ax.scatter([centre], [centre], marker="+", color="black", s=150, linewidths=2)
ax.annotate(
"home-range\ncentre",
(centre, centre),
xytext=(10, -25),
textcoords="offset points",
fontsize=9,
)
ax.set_aspect("equal", adjustable="box")
ax.set_xlabel("x (km)")
ax.set_ylabel("y (km)")
ax.set_title("One simulated moose, coloured by time (h)")
fig.colorbar(points, ax=ax, label="t (h)")
# (b) Final positions of all virtual moose against the analytic density.
ax = axes[0, 1]
ax.scatter(x_positions, y_positions, s=3, alpha=0.25, color="steelblue")
radius = np.linspace(0, 2 * math.pi, 200)
sd = stationary_sd(kappa=kappa, sigma=sigma)
for level in (1, 2):
ax.plot(
level * sd * np.cos(radius),
level * sd * np.sin(radius),
color="darkred",
linewidth=2,
)
ax.set_aspect("equal", adjustable="box")
ax.set_xlabel("x (km)")
ax.set_ylabel("y (km)")
ax.set_title("Final positions with 1 and 2 sd analytic contours")
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
# (c) 1D marginal: particles, Fokker-Planck PDE, and the analytic density.
ax = axes[1, 0]
bins = np.linspace(-4, 4, 41)
ax.hist(
x_positions,
bins=bins,
density=True,
color="steelblue",
alpha=0.6,
label="simulated moose",
)
x_fine = np.linspace(-4, 4, 400)
ax.plot(
x_fine,
analytic_pdf(x_fine, t_end, kappa=kappa, sigma=sigma, centre=centre),
color="black",
linewidth=2,
label="analytic density at t = T",
)
ax.plot(
x,
snapshots[-1],
color="darkred",
linestyle="--",
linewidth=2,
label="Fokker-Planck PDE",
)
ax.set_xlabel("x (km)")
ax.set_ylabel("density (1/km)")
ax.set_title(f"Marginal distribution at t = {t_end:.0f} h")
ax.legend(fontsize=9)
# (d) Relaxation of the variance towards its equilibrium value.
ax = axes[1, 1]
ax.plot(
times,
variances,
color="darkred",
linestyle="--",
label="Fokker-Planck PDE",
)
t_fine = np.linspace(0, t_end, 300)
ax.plot(
t_fine,
analytic_variance(t_fine, kappa=kappa, sigma=sigma),
color="black",
linewidth=2,
label="analytic",
)
ax.axhline(sigma**2 / (2 * kappa), color="grey", linewidth=1, label="equilibrium")
ax.set_xlabel("t (h)")
ax.set_ylabel("variance (km$^2$)")
ax.set_title("Relaxation of the density to equilibrium")
ax.legend(fontsize=9)
fig.suptitle(
"Ornstein-Uhlenbeck moose movement: particles versus the Fokker-Planck equation",
)
fig.tight_layout()
output_path = Path(output)
output_path.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(output_path, dpi=300, bbox_inches="tight", facecolor="white")
plt.close(fig)
return output_path
def run(
positions_path: str | Path = "data/moose_positions.csv",
track_path: str | Path = "data/moose_track.csv",
output: str | Path = "figures/moose_fokker_planck.png",
) -> Path:
"""Solve, compare, plot, and print diagnostics for the OU moose model."""
x_positions, y_positions = load_positions(positions_path)
t_track, x_track, y_track = load_track(track_path)
solve = solve_fokker_planck()
_x, _times, _snapshots, masses, means, variances = solve
output_path = plot_comparison(
x_positions, y_positions, t_track, x_track, y_track, solve, output
)
print_diagnostics(diagnostics(masses, means, variances))
print(f"\nSaved plot to: {output_path}")
return output_path
def main() -> None:
"""Console-script entry point."""
run()
if __name__ == "__main__":
main()
The three routes to the answer agree. The printed diagnostics are:
Fokker-Planck PDE versus analytic solution:
mass at T: 1.000000 (should be 1.000000)
PDE variance at T: 1.007964
analytic variance: 1.000000
PDE mean at T: -0.000000
analytic mean: 0.000000The mass is conserved to machine precision, as it must be for a no-flux boundary. The variance is high by about \(0.8\%\), and the reason is instructive: first-order upwinding of the drift term adds an artificial diffusion of size roughly \(|v|\,\Delta x/2\). Halving the grid spacing roughly halves that bias. This is a reminder that “solving the Fokker-Planck equation numerically” always means solving a discretization of it, and the discretization has its own physics.
Properties of the OU model
The OU model’s special structure is what makes it both useful and, in the end, restrictive.
- Exact solvability. Linear drift and constant diffusion keep the density Gaussian forever. The stationary and transition densities are normal, so predictions are closed-form.
- A Gaussian home range. The stationary distribution is a single isotropic bell centred on \(c\) with standard deviation \(s = \sigma/\sqrt{2\kappa}\). Its area is a home-range size, and its contours are a utilization distribution. It has no notion of shape beyond a circle.
- Exponential relaxation with a spectral gap. The Fokker-Planck operator for the OU model has Gaussian eigenfunctions with eigenvalues \(-n\kappa\), so every deviation from equilibrium decays at a known exponential rate and the slowest mode is \(1/\kappa\). Equilibrium is approached, and how fast is computable.
- Closed moment equations. The mean relaxes to \(c\) like \(e^{-\kappa t}\) and the variance approaches \(s^2\) like \(e^{-2\kappa t}\), so the first two moments can be tracked without solving the PDE at all.
- Reversibility. The OU process is in detailed balance: in equilibrium, every microscopic transition is balanced by its reverse, and there is no net circulation of probability. This is a strong restriction, and general models lose it.
- Positivity and conservation. As for any Fokker-Planck model, \(p \ge 0\) is preserved and total probability is conserved with no-flux boundaries. Adding a sink term turns it into a model of mortality or of animals leaving the study area.
- Infinite propagation speed. For \(\sigma > 0\) the equation is parabolic, so a point mass immediately has nonzero density everywhere. The moose could be anywhere, with minuscule but positive probability immediately, no matter how fast it would have to travel. This is in direct tension with the hard maximum-speed bound of the ellipse method.
- A diffusion limit. Without drift the model is Brownian motion and the variance grows linearly, \(2Dt\). The OU process is the mean-reverting version of that random walk, which is why it sits underneath both random-walk models and macroscopic diffusion models.
Strengths of the OU model
- Interpretable parameters. Fidelity and restlessness are statements about the animal, not statistical conveniences, and they carry units a field biologist can reason about.
- Closed-form answers for useful quantities. Home-range size, the time to settle, the overlap of two home ranges, and the expected time to encounter a fixed feature are all available without simulation.
- A utilization distribution for free. The stationary density is what home-range and habitat analyses actually want, produced from two individual-level parameters.
- It fits data. Because the transition density is Gaussian, the OU model is a linear Gaussian state-space model: it can be fitted to irregularly-spaced collar fixes by maximum likelihood or Kalman filtering, with measurement error handled explicitly (Hooten et al. 2017).
- It is the right null model. In movement ecology the mean-reverting diffusion is the baseline against which more elaborate models are compared, much as Brownian motion is the baseline in phylogenetics.
- It is a bridge. The OU process links individual movement rules to redistribution patterns through the continuum limit in the movement ecology literature (Turchin 1998; Nathan et al. 2008).
Defects of the OU model
Every strength above is also a defect when the animal does not cooperate, and moose mostly do not.
- The home range is a single Gaussian blob. Real ranges are often multi-modal, elongated along valleys and ridges, and asymmetrical. An isotropic OU process cannot represent any shape.
- The drift is linear. There is no habitat in the model. Food, cover, and terrain should bend the drift, not a spring to a point.
- The diffusion is constant and isotropic. Movement is equally easy in all directions and everywhere, which is false across snow, muskeg, forest, and open ground, and false along corridors.
- There are no boundaries. Rivers, lakes, highways, fences, and cliff lines are invisible to the model, even though they dominate where a moose can go.
- There are no interactions. The model describes one animal in a fixed environment. It cannot represent avoidance of conspecifics, a cow and calf staying together, or predator risk.
- There is no memory or behaviour. The process is Markovian with white noise and Gaussian steps. Real moose have behavioural states, autocorrelated moves, and heavy-tailed excursions that a single \(\sigma\) cannot produce.
- The parameters are constant. Real fidelity and restlessness vary with season, snow, rut, and time of day. A stationary model over months is an average that may describe no moose in particular.
- Only one animal at a time. Turning the density into a population model requires either independent copies (ignoring interaction) or coupling.
- It grants infinite speed, as noted above, so it is not a physical bound and can assign probability to impossible jumps.
The next post repairs these one at a time by writing the Fokker-Planck equation in its general form: a habitat potential instead of a spring, a spatial diffusion field instead of a constant, barriers as boundary conditions, interactions as nonlocal couplings, behaviour as multiple states, and jumps as fractional diffusion.
Conclusion
The Ornstein-Uhlenbeck process is the natural first model of a moose with a home range: two interpretable parameters, an exact Gaussian answer, a known relaxation rate, and a clean path to fitting collar data. It is also a caricature. The rest of this series generalizes the density equation until the model can carry habitat, barriers, interactions, behaviour, and jumps, using the OU model as the baseline whose defects define the agenda.