{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "(ODE_Lotka_Volterra_fit_multiple_ways)= \n", "# ODE Lotka-Volterra With Bayesian Inference in Multiple Ways\n", "\n", ":::{post} January 16, 2023\n", ":tags: ODE, PyTensor, gradient-free inference\n", ":category: intermediate, how-to\n", ":author: Greg Brunkhorst\n", ":::" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running on PyMC v5.1.2+24.gf3ce16f26\n" ] } ], "source": [ "import arviz as az\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import pymc as pm\n", "import pytensor\n", "import pytensor.tensor as pt\n", "\n", "from numba import njit\n", "from pymc.ode import DifferentialEquation\n", "from pytensor.compile.ops import as_op\n", "from scipy.integrate import odeint\n", "from scipy.optimize import least_squares\n", "\n", "print(f\"Running on PyMC v{pm.__version__}\")" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%load_ext watermark\n", "az.style.use(\"arviz-darkgrid\")\n", "rng = np.random.default_rng(1234)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Purpose\n", "The purpose of this notebook is to demonstrate how to perform Bayesian inference on a system of ordinary differential equations (ODEs), both with and without gradients. The accuracy and efficiency of different samplers are compared.\n", "\n", "We will first present the Lotka-Volterra predator-prey ODE model and example data. Next, we will solve the ODE using scipy.odeint and (non-Bayesian) least squares optimization. Next, we perform Bayesian inference in PyMC using non-gradient-based samplers. Finally, we use gradient-based samplers and compare results. \n", "\n", "### Key Conclusions\n", "Based on the experiments in this notebook, the most simple and efficient method for performing Bayesian inference on the Lotka-Volterra equations was to specify the ODE system in Scipy, wrap the function as a Pytensor op, and use a Differential Evolution Metropolis (DEMetropolis) sampler in PyMC. " ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Background\n", "### Motivation\n", "Ordinary differential equation models (ODEs) are used in a variety of science and engineering domains to model the time evolution of physical variables. A natural choice to estimate the values and uncertainty of model parameters given experimental data is Bayesian inference. However, ODEs can be challenging to specify and solve in the Bayesian setting, therefore, this notebook steps through multiple methods for solving an ODE inference problem using PyMC. The Lotka-Volterra model used in this example has often been used for benchmarking Bayesian inference methods (e.g., in this Stan [case study](https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html), and in Chapter 16 of *Statistical Rethinking* {cite:p}mcelreath2018statistical." ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "### Lotka-Volterra Predator-Prey Model\n", "The Lotka-Volterra model describes the interaction between a predator and prey species. This ODE given by:\n", "\n", "\n", "\\begin{aligned}\n", "\\frac{d x}{dt} &=\\alpha x -\\beta xy \\\\ \n", "\\frac{d y}{dt} &=-\\gamma y + \\delta xy\n", "\\end{aligned}\n", "" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "The state vector $X(t)=[x(t),y(t)]$ comprises the densities of the prey and the predator species respectively. Parameters $\\boldsymbol{\\theta}=[\\alpha,\\beta,\\gamma,\\delta, x(0),y(0)]$ are the unknowns that we wish to infer from experimental observations. $x(0), y(0)$ are the initial values of the states needed to solve the ODE, and $\\alpha,\\beta,\\gamma$, and $\\delta$ are unknown model parameters which represent the following: \n", "* $\\alpha$ is the growing rate of prey when there's no predator.\n", "* $\\beta$ is the dying rate of prey due to predation.\n", "* $\\gamma$ is the dying rate of predator when there is no prey.\n", "* $\\delta$ is the growing rate of predator in the presence of prey. " ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "### The Hudson's Bay Company data\n", "The Lotka-Volterra predator prey model has been used to successfully explain the dynamics of natural populations of predators and prey, such as the lynx and snowshoe hare data of the Hudson's Bay Company. Since the dataset is small, we will hand-enter the values." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
yearlynxhare
01900.04.030.0
11901.06.147.2
21902.09.870.2
31903.035.277.4
41904.059.436.3
\n", "