University of Southern California Dissertations and Theses

From least squares to Bayesian methods: refining parameter estimation in the Lotka-Volterra model

(USC Thesis Other)

From least squares to Bayesian methods: refining parameter estimation in the Lotka-Volterra model

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content FROM LEAST SQUARES TO BAYESIAN METHODS:
REFINING PARAMETER ESTIMATION IN THE LOTKA-VOLTERRA MODEL
by
Shengmiao Huang
A Thesis Presented to the
FACULTY OF THE USC DANA AND DAVID DORNSIFE COLLEGE OF
LETTERS ARTS AND SCIENCE
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
MASTER OF SCIENCE
(APPLIED MATHEMATICS)
May 2024
© 2024 Shengmiao Huang

TABLE OF CONTENTS
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
1 Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 The Lotka-Volterra Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Ordinary Differential Equations and Bayesian Inference . . . . . . . . . . . . 4
2 Chapter 2: Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1 Without Bayesian Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.1 Solving ODE Using Scipy.odeint Function . . . . . . . . . . . . . . 6
2.1.2 Least Square Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Gradient-free Bayesian Inference . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 The Bayesian Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.2 PyMC Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.3 Metropolis Sampler . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3 Chapter 3: Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.1 The Lynx and Hare Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Explore The Lynx and Hare Data Using The LV Model . . . . . . . . . . . . 13
3.3 Least Square Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.4 With Gradient-free Bayesian Inference . . . . . . . . . . . . . . . . . . . . . 16
4 Chapter 4: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
ii

LIST OF FIGURES
Figure 1.1: Exponential Growth Model vs. Logistic Growth Model . . . . . . . . 2
Figure 1.2: An Example of Classic Lotka-Volterra Model . . . . . . . . . . . . . . 3
Figure 3.1: Example Model Run . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Figure 3.2: Least Square Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Figure 3.3: Trace Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Figure 3.4: Trace Plot and Rank Plot . . . . . . . . . . . . . . . . . . . . . . . . 17
Figure 3.5: Trajectory Graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
iii

ABSTRACT
This study explores the application of Bayesian inference to the Lotka-Volterra model for
parameter estimation, comparing to traditional methods that often rely solely on point estimates. Initially, the Lotka-Volterra model, a fundamental representation in ecological modeling for predator-prey dynamics, is solved using classical techniques with the Scipy.odeint
function, followed by parameter refinement through least squares optimization. Subsequently, the research shifts to a Bayesian framework employing gradient-free methods using
the PyMC library, enabling a probabilistic approach to parameter estimation. This shift allows for the incorporation of prior knowledge and observational data, offering a more robust
understanding of parameter uncertainties. Results indicate that while traditional methods
provide quick and effective estimates, they may underestimate parameter variability and
exhibit bias under model misspecification. Conversely, Bayesian inference delivers a comprehensive probabilistic view, revealing not only the most likely estimates but also their feasible
ranges, thereby providing deeper insights into the ecological dynamics. The study underscores the ability of Bayesian methods in handling complex ecological models, particularly
in scenarios requiring a nuanced understanding of uncertainty for effective decision-making
in ecological management and conservation.
iv

Chapter 1: Introduction
1.1 The Lotka-Volterra Model
Mathematical models have long served as a beacon in the pursuit of understanding
the ebb and flow of life within an ecosystem, guiding ecologists through the complexities
of biological interactions. Among these models, the simplest yet profoundly influential are
the exponential growth model and logistic growth model. The exponential growth model
represents the growth of a population without any environmental constraints, that is, infinite
resources and spaces without the presence of natural enemies [1]. However, this condition is
ideal and barely met in reality. So the logistic growth model, which represents the growth of
a population within a constrained environment, was developed. The logistic growth model
offers an introductory glimpse into the mathematical domain of ecology [2].
The logistic growth equation
dP
dt = rP(1 −
P
K
)
elegantly captures the essence of resource-limited expansion, where
• P stands for the population size;
• r is the rate of maximum population growth;
• K is the carrying capacity of the environment
are all real numbers.
Over time, this model predicts a characteristic S-shaped curve, an initial exponential
growth that plateaus as the population nears the carrying capacity, intuitively mirroring
natural populations. Figure 1.1 below visually illustrates the difference between exponential
1

growth and logistic growth models:
Figure 1.1: Exponential Growth Model vs. Logistic Growth Model
From the seeds of understanding planted by logistic growth, we advance to the more
nuanced Lotka-Volterra model, a keystone in the archway to complex population dynamics.
The Lotka-Volterra model, often referred to as the predator-prey equations, represents a
fundamental concept within mathematical biology and ecology. Originally formulated independently by Alfred J. Lotka in 1925 and Vito Volterra in 1926, the model describes the
dynamics of biological systems in which two species interact, one as a predator and the
other as prey [3][4][5]. The equations provide a theoretical framework for understanding and
predicting the population cycles of the interacting species.
Figure 1.2 below illustrates a classic Lotka-Volterra model. We can see that, at first, the
population of the predators reaches peak and starts to decrease when the population of the
prey drops to a certain level (being predated); then when the population of predators drops
to a certain level (no prey to be predated), the population of the prey starts to increase again.
Then the cycle repeats thus the simple structured ecosystem involving these two species are
balanced [6].
2

Figure 1.2: An Example of Classic Lotka-Volterra Model
The model comprises a pair of first-order, nonlinear, differential equations:
dx
dt = αx − βxy
dy
dt = −γy + δxy
where
• x and y represent the prey and predator population sizes, respectively.
• α represents the growth rate of the prey population size without the presence of predators;
• β represents the decreasing rate of the prey population size in the presence of predators;
• γ represents the decreasing rate of the predator population size without the presence
of prey;
• δ represents the growth rate of the predator population size in the presence of prey.
3

are all real valued parameters.
One famous application in ecology is the study of the Canadian lynx and snowshoe hare
pelts data collected by the Hudson’s Bay Company over nearly a century. This dataset is
organized and analyzed in many different studies (by MacLulich, Elton&Nicholso, Odum, et
al.), which is also the dataset we are going to use in the simulation [7]. The cycles observed
in this data set were effectively modeled by the Lotka-Volterra equations, demonstrating the
model’s practical relevance in predicting natural phenomena.
Beyond ecological modeling, the Lotka-Volterra equations find applications in various
fields, showcasing their versatility in modeling competitive and predatory dynamics in different systems. Mahtani in his paper adapts the Lotka–Volterra model to simulate COVID–19
spread, treating the infected individuals as predators and susceptible individuals as prey,
demonstrating that the model can predict basic dynamic behaviors of the pandemic and
thus offering an application of the Lotka–Volterra model in epidemiology [8]. Moreover,
Comesa applied the Lotka-Volterra model to the banking sector, modeling interactions between Mother Banks, Subsidiary Banks, and Individuals/Companies. The model uses differential equations to analyze capital transfers and dynamic equilibrium, demonstrating the
versatility of the Lotka-Volterra framework in financial contexts [9].
1.2 Ordinary Differential Equations and Bayesian Inference
The pair of equations in the Lotka-Volterra Model above is a classical form of Ordinary
Differential Equations (ODEs). The ODE models are pivotal in science and engineering for
modeling the time evolution of physical variables across diverse domains.
Given their complexity and the intrinsic uncertainties in real-world data, Bayesian inference emerges as a natural choice for estimating ODE model parameters. This probabilistic
approach not only quantifies the uncertainty in parameter estimates but also incorporates
prior knowledge and observational data to refine these estimates.
In this thesis, I will employ the Lotka-Volterra model to illustrate how Bayesian inference can be applied to systems of ODEs to estimate the parameters of such models more
robustly, compared to traditional methods that often provide point estimates and require
4

assumptions about parameter distributions and error terms. This method incorporates prior
knowledge and observational data within a coherent statistical framework, yielding distributions of possible parameter values rather than single estimates. This approach enhances
the understanding of model uncertainty and reliability, which is crucial for making informed
decisions in ecology and resource management.
5

Chapter 2: Methods
2.1 Without Bayesian Inference
2.1.1 Solving ODE Using Scipy.odeint Function
Scipy.odeint is a function from the Scipy.integrate module of Scipy library, which
is used to solve systems ODEs.
Specifically, the function odeint solves ODEs in the form of dy
dt = f(y, t), where y is
a vector of dependent variables, and t is the time. Regarding our intention, the dependent
variable vector y represents the population Canadian Lynx and snowshoe hare, and t represents time over years. The function odeint will integrate these equations from one point in
time to another and return the y values for each requested point. For two species, we can
define our ODE model as follows:
def model(Y, t, params):
x, y = Y
alpha, beta, gamma, delta = params
dxdt = alpha * x - beta * x * y
dydt = -gamma * y + delta * x * y
return [dxdt, dydt]
For function odeint to work, we need to input initial parameters according to our system
of equations: θ = [α, β, γ, δ, x(0), y(0)], where x(0) and y(0) are the initial populations of the
two species we are interested in. We also need to input an array of time points at which the
solver should provide the values of y, and the first element of this array must be the point
at which the initial conditions are given.
6

There could be multiple ways of arriving at initial guesses. For example, some published literature may provide initial estimates for parameters based on findings from similar
ecosystems or species; we could also try out some different values of the parameters to see
which combination better fits the observed data points.
2.1.2 Least Square Solution
After successfully solving the Lotka-Volterra ODEs using the Scipy.odeint function,
our next step is to refine the model parameters to fit the observed data better. This transition
from ODE solving to parameter estimation is crucial for enhancing the model’s accuracy and
reliability.
To achieve this, we employ the least squares optimization method, integrated seamlessly
with the ODE solver to minimize the discrepancy between the model’s predictions and the
actual data. The Scipy.optimize.least squares function is particularly suited for this
task as it is designed to handle complex models and can effectively deal with the nonlinear
nature of the Lotka-Volterra equations.
In this setup, the Scipy.odeint function calculates the predicted populations of the
species over time based on initial guesses for the parameters. These predictions are then
compared to the actual observed data points. The least squares method adjusts the model
parameters iteratively, aiming to minimize the sum of the squares of the differences (residuals) between the observed values and the model predictions. Each iteration involves recalculating the model predictions with the updated parameters until the best-fit parameters,
which minimize the residuals, are found.
2.2 Gradient-free Bayesian Inference
The least square method in general provides a good-fitting model to the given datasets,
which can be used for prediction and analysis purpose in many applications. However, least
square method is a point estimation, which does not take the uncertainty of parameters into
consideration. On the other hand, Bayesian inference techniques provide a probabilistic view
of parameter estimates, usually in the form of posterior distributions. These distributions
7

provide a full picture of the uncertainty around parameter estimates, showing the most likely
value and how likely all possible values are [10].
In transitioning to a Bayesian framework using PyMC, we establish a probabilistic model
that incorporates our prior beliefs about the parameters and the likelihood of observing the
data given these parameters. This approach not only provides a rigorous statistical basis for
parameter estimation but also quantifies the uncertainty associated with these estimates.
2.2.1 The Bayesian Model
To those who might not be familiar with Bayesian models, it is important to understand
how they are set up and how they work.
• Setting Up the Bayesian Model
The process starts with defining prior distributions for each parameter in the LotkaVolterra equations. Priors represent our previous knowledge or assumptions about the
values these parameters might take before considering the current data. The choice of
priors is crucial as it can significantly influence the posterior outcomes, especially in
cases where data are sparse or noisy.
• Specifying Priors
For each parameter in our model — namely, the growth rates (α and γ) and interaction
coefficients (β and δ) — we assign a prior distribution. These distributions are typically
chosen based on historical data, expert knowledge, or previous studies. In our case, we
select normal distributions for these priors due to their mathematical convenience and
the central tendency they imply, which is a reasonable assumption in many biological
systems.
• Likelihood Function
After defining the priors, we integrate the observed data into our model through the
likelihood function. This function calculates the probability of observing the data
given the parameters specified by the model. For the Lotka-Volterra model, where
the data involve species counts over time, a normal distribution is typically used for
the likelihood, assuming that the data are normally distributed around the model’s
predictions with some measurement error.
8

• Model Compilation and Sampling
Finally, the model is compiled, and posterior sampling is conducted using PyMC’s
built-in sampling algorithms. This step involves generating samples from the posterior distributions of the parameters, which incorporate both our prior beliefs and the
information from the observed data.
By setting up the Bayesian model in this manner, we not only estimate the parameters
but also obtain a full probabilistic description of the parameter space, which is essential for
making informed decisions under uncertainty in reality applications.
2.2.2 PyMC Model
We utilize PyMC, a powerful open-source programming library for Bayesian statistical
modeling and probabilistic machine learning, for Lotka-Voltarre model parameter estimation
with Bayesian inference. PyMC was firstly developed by Chris Fonnesbeck starting from
2003, and then joint by Huard, Patil, and Salvatier for the later development of PyMC 2.0
and 3.0 [11][12][13]. It uses advanced sampling algorithms like Markov Chain Monte Carlo
and integrates with scientific computing libraries such as Numpy and Scipy. It is particularly
popular in application fields such as ecology, where complex data analysis is required [14].
However, the function scipy.integrate.odeint we used cannot be directly used in
the PyMC model because it requires the input and output types of the variable to compile.
Therefore, we need first to tell PyMC the variable types using the Pytensor wrapper for the
function to be used in the PyMC model.
Firstly, we define the type of data we are inputting and what we expect to get out. We
use the @as op decorator to specify that both inputs and outputs are double float tensors:
@as_op(itypes=[pt.dvector], otypes=[pt.dmatrix])
def pymc_ode_model(params):
return odeint(func=model, y0=params[-2:], t=data.year, args=(params[:-2],))
Next, we implement the pymc ode model function, which encapsulates the call to Scipy.odeint.
Then, by imputing theta, we can perform numerical integration to obtain the model predictions over time.
9

Then, we can select the PyMC model using the ODE solver. Based on the results of the
least square solution we obtained earlier, we assign the prior distribution to the parameters,
which help inform the model about the range within which the parameters are likely to fall.
with pm.Model() as pymc_model:
# Define priors for the ODE parameters and initial conditions
alpha = pm.TruncatedNormal("alpha", mu=params[0],
sigma=0.1, lower=0)
beta = pm.TruncatedNormal("beta", mu=params[1],
sigma=0.01, lower=0)
gamma = pm.TruncatedNormal("gamma", mu=params[2],
sigma=0.1, lower=0)
delta = pm.TruncatedNormal("delta", mu=params[3],
sigma=0.01, lower=0)
x0 = pm.Normal("x0", mu=params[4], sigma=1)
y0 = pm.Normal("y0", mu=params[5], sigma=1)
# Define the standard deviation of the observation error
sigma_error = pm.HalfNormal("sigma_error", sigma=10)
# Combine all parameters into a single vector
parameters = pm.math.stack([alpha, beta, gamma, delta, x0, y0])
# ODE model solution wrapped with Pytensor operator
ode_solution = pymc_ode_model(parameters)
We define the likelihood function using the normal distribution and give the observed
data to this likelihood function. PyMC is then used to draw samples from the posterior
distribution of the parameters.
Y_obs = pm.Normal("Y_obs", mu=ode_solution, sigma=sigma_error,
observed=data[[’hare’, ’lynx’]].values)
10

2.2.3 Metropolis Sampler
This is a traditional Markov Chain Monte Carlo (MCMC) algorithm that proposes small
random steps from the current parameter state and accepts or rejects these steps based on
how likely the resulting parameters are given the observed data. It allows for sampling from
a probability distribution by constructing a Markov Chain with the desired distribution as
its equilibrium distribution [15].
Even though there is ready-to-use MCMC algorithm available, it is essential to understand what is it and how it works. The concept of MCMC was originated with Andrey
Markov in the early 1900s, and the methods became prominent with the Metropolis algorithm in 1953 (Metropolis et al.) and the Metropolis–Hastings algorithm in 1970 (Hastings)
[16][17]. As not mentioned in details in Brunkhorst’s paper [22], I will walk through the
theory of MCMC to give a general idea for people who might not be familiar with it.
Specifically, we start by initializing the algorithm with the theta we obtained from the
least squared solution step. Then, for each iteration, we generate a new sample from a
normal distribution. We choose normal distribution because it is symmetric around the
current point, which is the mean of the distribution, and the standard deviation is used as
a tuning parameter.
In implementing the Metropolis sampler for our Bayesian inference, a critical step is the
selection of the proposal distribution’s variance. The variance of the proposal distribution
plays a crucial role in determining the efficiency and effectiveness of the sampling process.
A too small variance leads to a high acceptance rate but poor exploration of the parameter
space, resulting in correlated samples and slow convergence. Conversely, a too large variance
may result in a low acceptance rate, causing the sampler to reject many proposals and
stagnate, failing to explore significant regions of the parameter space.
To optimize this variance, we typically start with a pilot run of the sampler, using an
initial guess based on theoretical considerations or empirical data. During this pilot run, we
monitor the acceptance rate. An ideal target for the acceptance rate is around 25% for singleparameter problems or up to 50% when sampling from many parameters simultaneously.
Adjustments to the variance are then made based on this observed acceptance rate. If the
11

rate is too low, we decrease the variance; if it’s too high, we increase the variance and then
iterate this process until the acceptance rate falls within the desired range.
Then we calculate an acceptance ratio, which is the probability of moving to the new
sample compared to the probability of staying at the current sample. This is determined
by the likelihood of the observed data under the proposed parameters and the priors on
the parameters. The ratio is given by P(new sample)/P(current sample), where P is the
probability law of the parameter value under the targeted distribution.
If the ratio is greater or equal to 1, the new sample is more probable, and we accept
the move; if the ratio is less than 1, then the new sample is less probable, and we accept the
move with probability equal to the ratio. To do this, we draw a uniform ramble number u
between 0 and 1. If u is less than or equal to the ratio, we move to the new point; otherwise,
we stay at the current point.
We iterate as many times as we desire. Initially, samples are typically discarded to allow
the chain to converge to the target distribution (posterior distribution). The subsequent
samples are considered from the target distribution, and the collection of these valid samples
forms an empirical distribution that approximates the target distribution.
with pymc_model:
all_vars = pymc_model.free_RVs
trace_M = pm.sample(3000, step=pm.Metropolis(vars=all_vars))
12

Chapter 3: Results and Discussion
3.1 The Lynx and Hare Data
We use the methods described above to analyze the lynx and hare data. This data set,
obtained from MacLulich in 1937 [18], was primarily compiled from Hudson’s Bay Company’s
historical pelt-trading records of Canadian lynx and snowshoe hare. The hare data before
1903 was from fur records, while data after 1903 was from questionnaires. The lynx data post1903 was also incomplete as presented by Elton and Nicholson in 1942 [19]. The presentation
of this hare–lynx data varies across authors, with Odum re-editing MacLulich’s data in 1953,
leading to slight difference [20].
In this paper, I will use the ready-to-use data compiled and organized in the paper I’m
doing the literature review on by Greg Brunkhorst [22].
3.2 Explore The Lynx and Hare Data Using The LV Model
We use the two dimensional LV-model to analyze the lynx and hare data. First we
initialize the parameters as
params = [0.52, 0.026, 0.84, 0.026]
y0 = [40, 5]
t = np.arange(1900, 1921, 0.01)
where 40 and 5 indicate the “guessed” hare and lynx population size in year 1900
respectively. For the values of the interaction parameters α, β, γ, and δ, I used Brunkhorst’s
paper as reference for the initial run [22].
13

Then we can solve the Lotka-Volterra ODEs with the initial values using the Scipy.odeint
function and get an idea about the model as shown in Figure 3.1:
result = odeint(model, y0, t, args=(params,))
Figure 3.1: Example Model Run
The model somewhat captures the data points, but can be further improved.
From the model with our reasonable initial guess, the sum of squared residuals for Lynx
population is 656.62, and the sum of squared residuals for Hare population 1867.01.
3.3 Least Square Solution
We use the following command to find the least square estimate of the parameters on
nthe lynx and hare data.
def residual(theta):
return (
data[["hare", "lynx"]] -
odeint(func=model, y0=[40, 5], t=data.year, args=(params,))
).values.flatten()
result = least_squares(residual, params)
The optimal values provided by the least square solution for our θ are
• α ≈ 0.469;
14

• β ≈ 0.023;
• γ ≈ 0.893;
• δ ≈ 0.025,
which gives a better fitting graph as shown in Figure 3.2 below:
Figure 3.2: Least Square Solution
In the model provided by the least square solution method, the sum of squared residuals
for Lynx population is 198.51, and the sum of squared residuals for Hare is 639.32, which
are much less than the residuals provided by the model using reasonable initial guess at the
beginning. We observed that it provides a pretty good fit to the data and a clear visualization
for interpretation purpose.
The least squares method is widely used in parameter estimation due to its simplicity
and efficiency. However, it is important to understand the limitations and assumptions
associated with this method. One key assumption of the least squares method is that the
model errors are normally distributed. When these assumptions are violated, the parameter
estimates can be biased, and the confidence intervals can be inaccurate.
Furthermore, the least squares method is sensitive to outliers and can be influenced by
extreme values in the data. This is particularly relevant in ecological studies like the one
we explored in this paper, where outliers and variations are common. Robustness to such
variations is crucial for reliable parameter estimation.
In conclusion, while the least squares method offers a computationally efficient and
straightforward approach for parameter estimation, it is essential to consider its limitations
and the context of the data being analyzed.
15

3.4 With Gradient-free Bayesian Inference
In assessing the convergence of our MCMC simulations, we employ diagnostic tools to
ensure the reliability of our parameter estimates. Trace plots and Rank plots are used to
visually inspect the stability of the sampling paths, looking for signs of convergence to a
steady state without drift. The diagnostic is critical for confirming the robustness of our
Bayesian inference process, ensuring that the posterior distributions are well-characterized
and representative of the underlying data [21].
Figure 3.3 and Figure 3.4 provides a numerical and visual presentation of the distribution
of each parameter simulated by the MCMC process:
Figure 3.3: Trace Summary
• hdi 3% and hdi 97%: the 3% and 97% highest density intervals, representing the range within
which 94% of the parameter values fall, providing a measure of the parameter’s credible interval.
• mcse sd: the Monte Carlo standard error for the standard deviation, indicating the error in
the estimate of the standard deviation due to sampling variability.
• ess bulk: the effective sample size for the bulk of the distribution, reflecting te number of
independent sample with high auto-correlation.
• ess tail: the effective sample size for the tail of the distribution, indicating the number of
independent samples in the extremes of the distribution.
• r hat: the potential scale reduction factor, used to diagnose convergence. Values close to 1
indicate that the chains have converged.
16

Figure 3.4: Trace Plot and Rank Plot
The trace plots on the left are bell-shaped around the mean for each parameter, and the rank plots
on the right are pretty well leveled for each chain of each parameter, which indicates the sampled
values for each parameter of the MCMC simulation converges to the posterior distribution.
17

The graphs on the lefts are trace plots, which show the sequence of sampled values
for each parameter across the iterations of the MCMC simulation. The curves are a little
wiggly but do form a band around a central value without trending upward or downward, it
indicates that the chain has likely converged to the posterior distribution.
The graphs on the right are rank plots, where the colors represent different chains, and
the vertical position of each dot shows the rank of a sample when all chains are merged and
sorted. For example, if ten values, from one to ten, are sampled from the simulation, then
we place the lowest value one at rank one and the highest value ten at rank ten. The points
somewhat form a horizontal line, indicating that the ranks are uniformly distributed across
chains, thus convergence may have been reached.
Then we use the simulated parameters to draw a trajectory graph as shown in Figure
3.5:
Figure 3.5: Trajectory Graph
The “shaded” regions around the model predictions indicate the uncertainty associated with these
estimates. These regions are generated by sampling the trajectories from the posterior distribution
of the parameters.
The spread of these trajectories provides insight into the confidence we can place in the
model’s prediction. A narrow shaded region indicates high confidence, while a wider region
suggests greater uncertainty.
This graph demonstrates the power of Bayesian inference in ecological modeling. It
allows for a robust estimation of model parameters, providing insights into the underlying
biological processes and helping to understand the uncertainty in model predictions.
In conclusion, the Bayesian inference, though computationally more demanding due to
the iterative nature of MCMC methods, provided a deeper insight into the parameter dynam18

ics. The posterior distributions obtained from Bayesian inference revealed not only the most
likely estimates of the parameters but also their feasible ranges, offering a comprehensive
view of the model’s behavior under uncertainty.
19

Chapter 4: Conclusion
The findings from this case study underscore the difference between gradient-free Bayesian
inference and traditional least squares solutions in complex ecological models like the LotkaVolterra. The least squares solution method, though computationally more convenient and
easier to understand and interpret in many applications, it ignores the variability in the
parameters. The Bayesian approach, in contrast, is significantly more computationally expensive and could be complicated to understand for people who are not familiar with statistical models. But with its robustness to model assumptions and ability to incorporate prior
knowledge, the Bayesian approach offers a more nuanced and comprehensive understanding
of ecological dynamics. These insights are invaluable in ecological management and conservation strategies, where understanding the full spectrum of potential outcomes is crucial for
effective decision-making.
20

References
[1] Snider, S. B. & Brimlow, J. N. (2013) An Introduction to Population Growth. Nature
Education Knowledge 4(4):3
[2] Vandermeer, J. (2010) How Populations Grow: The Exponential and Logistic Equations.
Nature Education Knowledge 3(10):15
[3] Kingsland S. Alfred J. Lotka and The Origins of Theoretical Population Ecology. Proc
Natl Acad Sci U S A. 2015 Aug 4;112(31):9493-5. doi: 10.1073/pnas.1512317112. PMID:
26243872; PMCID: PMC4534218.
[4] Lotka AJ. Elements of Mathematical Biology. Dover; New York: 1956.
[5] Volterra V. Fluctuations in The Abundance of a Species Considered Mathematically.
Nature. 1926;118(2972):558–560.
[6] The Lotka-Volterra Predator-Prey Model. Hong Kong University of Science and Technology, 17 July 2022, https://math.libretexts.org/@go/page/93493.
[7] Zhang, Zhibin & Tao, yi & Li, Zhenqing. (2007). Factors Affecting Hare-Lynx Dynamics
in The Classic Time Series of The Hudson Bay Company, Canada. Climate Research.
34. 83-89. 10.3354/cr034083.
[8] Mahtani, Krish. “An Epidemiological Application of the Lotka-Volterra Model to Predict Population Dynamics of COVID-19.” Parabola, vol. 59, no.1, 2023
[9] Comes, C˘alin-Adrian. “Banking System: Three Level Lotka-Volterra Model.” Procedia
Economics and Finance, vol. 3, 2012, pp. 251-255. International Conference Emerg21

ing Markets Queries in Finance and Business. DOI: https://doi.org/10.1016/S2212-
5671(12)00148-7.
[10] Taboga, Marco (2021). “Bayesian Inference”, Lectures on Probability Theory and Mathematical Statistics. Kindle Direct Publishing. Online appendix.
https://www.statlect.com/fundamentals-of-statistics/Bayesian-inference.
[11] Salvatier, Wiecki & Fonnesbeck (2016) Salvatier J, Wiecki TV, Fonnesbeck C. Probabilistic Programming in Python Using PyMC3. PeerJ Computer Science. 2016;2(2):e55.
doi: 10.7717/peerj-cs.55.
[12] “History.” PyMC, The PyMC Development Team, 2024,
www.pymc.io/about/history.html.
[13] “PyMC: Past, Present, and Future.” PyMC, The PyMC Development Team, 2024,
www.pymc.io/blog/PyMC Past Present Future.html.
[14] Abril-Pla O, Andreani V, Carroll C, Dong L, Fonnesbeck CJ, Kochurov M, Kumar R,
Lao J, Luhmann CC, Martin OA, Osthege M, Vieira R, Wiecki T, Zinkov R. PyMC: a
Modern, and Comprehensive Probabilistic Programming Framework in Python. PeerJ
Comput Sci. 2023 Sep 1;9:e1516. doi: 10.7717/peerj-cs.1516. PMID: 37705656; PMCID:
PMC10495961.
[15] Taboga, Marco (2021). “Markov Chain Monte Carlo (MCMC) Methods”, Lectures
on Probability Theory and Mathematical Statistics. Kindle Direct Publishing. Online
appendix. https://www.statlect.com/fundamentals-of-statistics/Markov-Chain-MonteCarlo.
[16] Metropolis, N., et al. “Equation of State Calculations by Fast Computing Machines.”
The Journal of Chemical Physics, vol. 21, no. 6, 1953, pp. 1087-1092.
[17] Hastings, W. K. “Monte Carlo Sampling Methods Using Markov Chains and Their
Applications.” Biometrika, vol. 57, no. 1, 1970, pp. 97-109.
22

[18] MacLulich, D. A. Fluctuations in the Numbers of the Varying Hare (Lepus Americanus), University of Toronto Press, 1937, pp. 5–136. JSTOR,
http://www.jstor.org/stable/10.3138/j.ctvfrxkmj.3.
[19] Elton, Charles, and Mary Nicholson. “The Ten-Year Cycle in Numbers of the Lynx
in Canada.” Journal of Animal Ecology, vol. 11, no. 2, 1942, pp. 215–44. JSTOR,
https://doi.org/10.2307/1358.
[20] E. P. Odum (1953), Fundamentals of Ecology, Philadelphia, W. B. Saunders.
[21] Taboga, Marco (2021). “Markov Chain Monte Carlo (MCMC) Diagnostics”, Lectures
on Probability Theory and Mathematical Statistics. Kindle Direct Publishing. Online
appendix. https://www.statlect.com/fundamentals-of-statistics/Markov-Chain-MonteCarlo-diagnostics.
[22] Greg Brunkhorst . “ODE Lotka-Volterra With Bayesian Inference in Multiple Ways”.
In: PyMC Examples. Ed. by PyMC Team. DOI: 10.5281/zenodo.5654871
[23] Christopher H Remien, Mariah J Eckwright, and Benjamin J Ridenhour. Structural
Identifiability of The Generalized Lotka–Volterra Model for Microbiome Studies. Royal
Society Open Science, 8(7):201378, 2021.
[24] Heilman, Steven. MATH 541B, Graduate Mathematical Statistics II, Fall 2023. 2023.
PDF file.
[25] Jacob D Davis, Daniel V Olivenca, Sam P Brown, and Eberhard O Voit. Methods of
Quantifying Interactions Among Populations Using Lotka-Volterra Models. Frontiers in
Systems Biology, 2:1021897, 2022.
[26] Jerome Mounier, Christophe Monnet, Tatiana Vallaeys, Roger Arditi, Anne Sophie
Sarthou, Arnaud Helias, and Francoise Irlinger. Microbial Interactions within a Cheese
Microbial Community. Applied and Environmental Microbiology, 74(1):172–181, 2008.
[27] Vanni Bucci, Belinda Tzen, Ning Li, Matt Simmons, Takeshi Tanoue, Elijah Bogart,
Luxue Deng, Vladimir Yeliseyev, Mary L Delaney, Qing Liu, et al. Mdsine: Microbial
23

Dynamical Systems Inference Engine for Microbiome Time Series Analyses. Genome
biology, 17:1–17, 2016.
[28] S. P. Hastings and J. B. McLeod Classical Methods in Ordinary Differential Equations:
With Applications to Boundary Value Problems, Amer. Math. Soc., Rhode Island, 2011.
[29] Taboga, Marco (2021). “Bayesian Estimation of The Parameters of The Normal
Distribution”, Lectures on Probability Theory and Mathematical Statistics. Kindle Direct Publishing. Online appendix. https://www.statlect.com/fundamentals-ofstatistics/normal-distribution-Bayesian-estimation.
[30] Navarro, Danielle. 2023. “The Metropolis-Hastings Algorithm.” April 12, 2023.
https://blog.djnavarro.net/posts/2023-04-12 metropolis-hastings.
24

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Stochastic Variational Inference as a solution to the intractability problem of Bayesian inference

PDF

Data-driven learning for dynamical systems in biology

PDF

Increase colorectal cancer prediction accuracy with the influence (I)-score

PDF

Parameter estimation in second-order stochastic differential equations

PDF

Obtaining breath alcohol concentration from transdermal alcohol concentration using Bayesian approaches

PDF

Scalable latent factor models for inferring genetic regulatory networks

PDF

Robust estimation of high dimensional parameters

PDF

Population modeling and Bayesian estimation for the deconvolution of blood alcohol concentration from transdermal alcohol biosensor data

PDF

Prohorov Metric-Based Nonparametric Estimation of the Distribution of Random Parameters in Abstract Parabolic Systems with Application to the Transdermal Transport of Alcohol

PDF

An abstract hyperbolic population model for the transdermal transport of ethanol in humans: estimating the distribution of random parameters and the deconvolution of breath alcohol concentration

PDF

New methods for asymmetric error classification and robust Bayesian inference

PDF

Analysis of ergodic and mixing dynamical systems

PDF

The existence of absolutely continuous invariant measures for piecewise expanding operators and random maps

PDF

Nonparametric estimation of an unknown probability distribution using maximum likelihood and Bayesian approaches

PDF

Bayesian models for a respiratory biomarker with an underlying deterministic model in population research

PDF

A survey on the computational hardness of linear-structured Markov decision processes

PDF

Finite dimensional approximation and convergence in the estimation of the distribution of, and input to, random abstract parabolic systems with application to the deconvolution of blood/breath al...

PDF

Small area cancer incidence mapping using hierarchical Bayesian methods

PDF

A Bayesian approach for estimating breath from transdermal alcohol concentration

PDF

A nonlinear pharmacokinetic model used in calibrating a transdermal alcohol transport concentration biosensor data analysis software

Asset Metadata

Creator Huang, Shengmiao (author)

Core Title From least squares to Bayesian methods: refining parameter estimation in the Lotka-Volterra model

School College of Letters, Arts and Sciences

Degree Master of Science

Degree Program Applied Mathematics

Degree Conferral Date 2024-05

Publication Date 06/10/2024

Defense Date 06/05/2024

Publisher Los Angeles, California (original), University of Southern California (original), University of Southern California. Libraries (digital)

Tag Bayesian inference,least square solution,Lotka-Volterra model,MCMC,OAI-PMH Harvest,parameter estimation

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Sun, Fengzhu (committee chair), Haskell, Cymra (committee member), Heilman, Steven (committee member)

Creator Email shengmia@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC113992547

Unique identifier UC113992547

Identifier etd-HuangSheng-13067.pdf (filename)

Legacy Identifier etd-HuangSheng-13067

Document Type Thesis

Format theses (aat)

Rights Huang, Shengmiao

Internet Media Type application/pdf

Type texts

Source 20240610-usctheses-batch-1166 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu