# Mathematical modelling and uncertainty quantification for analysis of
biphasic coral reef recovery patterns

David J. Warne111To whom correspondence should be addressed. E-mail:
<EMAIL_ADDRESS>School of Mathematical Sciences, Queensland University
of Technology, Brisbane, Queensland 4001, Australia  Centre for Data Science,
Queensland University of Technology, Brisbane, Queensland 4001, Australia
Kerryn Crossman  Australian Institute of Marine Science, Townsville,
Queensland Australia  Grace E. M. Heron  School of Mathematical Sciences,
Queensland University of Technology, Brisbane, Queensland 4001, Australia
Centre for Data Science, Queensland University of Technology, Brisbane,
Queensland 4001, Australia  Jesse A. Sharp  School of Mathematical Sciences,
Queensland University of Technology, Brisbane, Queensland 4001, Australia
Centre for Data Science, Queensland University of Technology, Brisbane,
Queensland 4001, Australia  Wang Jin  The Kirby Institute, University of New
South Wales, Sydney, New South Wales, Australia  Paul Pao-Yen Wu  School of
Mathematical Sciences, Queensland University of Technology, Brisbane,
Queensland 4001, Australia  Centre for Data Science, Queensland University of
Technology, Brisbane, Queensland 4001, Australia  Matthew J. Simpson  School
of Mathematical Sciences, Queensland University of Technology, Brisbane,
Queensland 4001, Australia  Centre for Data Science, Queensland University of
Technology, Brisbane, Queensland 4001, Australia  Kerrie Mengersen  School of
Mathematical Sciences, Queensland University of Technology, Brisbane,
Queensland 4001, Australia  Centre for Data Science, Queensland University of
Technology, Brisbane, Queensland 4001, Australia  Juan-C. Ortiz  Australian
Institute of Marine Science, Townsville, Queensland Australia

###### Abstract

Coral reefs are increasingly subjected to major disturbances threatening the
health of marine ecosystems. Substantial research is underway to develop
intervention strategies that assist reefs in recovery from, and resistance to,
inevitable future climate and weather extremes. To assess potential benefits
of interventions, mechanistic understanding of coral reef recovery and
resistance patterns is essential. Recent evidence suggests that more than half
of the reefs surveyed across the Great Barrier Reef (GBR) exhibit deviations
from standard recovery modelling assumptions when the initial coral cover is
low ($\leq 10$%). New modelling is necessary to account for these observed
patterns to better inform management strategies. We consider a new model for
reef recovery at the coral cover scale that accounts for biphasic recovery
patterns. The model is based on a multispecies Richards’ growth model that
includes a change point in the recovery patterns. Bayesian inference is
applied for uncertainty quantification of key parameters for assessing reef
health and recovery patterns. This analysis is applied to benthic survey data
from the Australian Institute of Marine Science (AIMS). We demonstrate
agreement between model predictions and data across every recorded recovery
trajectory with at least two years of observations following disturbance
events occurring between 1992–2020. This new approach will enable new insights
into the biological, ecological and environmental factors that contribute to
the duration and severity of biphasic coral recovery patterns across the GBR.
These new insights will help to inform managements and monitoring practice to
mitigate the impacts of climate change on coral reefs.

Keywords: coral reef recovery; biphasic recovery patterns; reef recovery
modelling; climate change; marine monitoring.

##  1 Introduction

Environmental stressors, including climate change, continue to drive increased
frequency and severity of disturbance events affecting coral reefs worldwide.
This is leading to long-term declines in coral cover and changes in community
composition (De’ath et al.,, 2009, 2012; Duran et al.,, 2017; Johns et al.,,
2014; Lapointe et al.,, 2019; Tebbett et al.,, 2023). To ensure the future
health of coral reefs and connected marine ecosystems we need to continue
developing capabilities to effectively prioritise targeted intervention and
restoration efforts (McLeod et al.,, 2022). Forecasts of coral reef health
under different potential environmental scenarios are an essential part of
this decision-making process to ensure limited resources are efficiently
allocated (Thompson and Dolman,, 2010; Vercelloni et al.,, 2020; van Woesik et
al.,, 2018). In addition, predictions of expected versus observed coral cover
is an important health indicator of recovery (Castro-Sanguino et al.,, 2021;
Logan et al.,, 2020; Thompson et al.,, 2020). To produce reliable forecasts
and predictions for decision-makers, we need models that are informed by our
understanding of coral recovery patterns of reefs following perturbations by
disturbance events, such as storm swells (Cheal et al.,, 2017; De’ath et al.,,
2012; Fabricius et al.,, 2008), predation by crown-of-thorns starfish,
Acanthaster spp. (COTS) (De’ath et al.,, 2012; Pratchett et al.,, 2021;
Seymour and Bradbury,, 1999; Vercelloni et al.,, 2017), and acute periods of
thermal stress resulting in coral bleaching (Hughes et al.,, 2017; Mumby et
al.,, 2021).

Recently, Warne et al., (2022) demonstrated evidence for biphasic recovery
patterns amongst almost $50\%$ of reefs recovering from very low cover ($\leq
10\%$), across Australia’s Great Barrier Reef (GBR). Such recovery patterns
were characterised by a transient period of reduced recovery that persists for
approximately four years on average (Warne et al.,, 2022). Such recovery
patterns represent deviations from standard recovery modelling approaches that
are based on variations of classical population growth models, such as
exponential, Gompertz, or Logistic growth models (MacNeil et al.,, 2019; Ortiz
et al.,, 2018; Simpson et al.,, 2022; Thompson and Dolman,, 2010; van Woesik
et al.,, 2018). Warne et al., (2022) also demonstrate coral recovery forecasts
that ignore biphasic patterns could be overestimating future median cover
estimates by up to 22% in absolute coral cover. Since the statistical
methodology applied by Warne et al., (2022) is based upon multi-segment
regression analysis (Aminikhanghahi and Cook,, 2016; Hinkley,, 1970; Jin et
al.,, 2017), the study was restricted to a small subset of severely disturbed
reefs with at least 5 observations of recovery. This motivates the development
of a recovery modelling framework that can account for the effects of biphasic
recovery in a broader range of scenarios.

In the biological and ecological modelling literature, biphasic models have
been widely studied and applied in a variety of settings (Bodnar et al.,,
2013; Brouwer et al.,, 2017; Honsey et al.,, 2016; de A. Møller et al.,, 2021;
Phaiboun et al.,, 2015), the most general of which has been demonstrated
recently by Murphy et al., (2022). Here, we present a multi-species biphasic
coral recovery model building on the modelling frameworks of Murphy et al.,
(2022) and Simpson et al., (2023). Using Bayesian analysis for model
calibration, prediction, and uncertainty quantification, we demonstrate the
approach is directly applicable to inform health assessment of the reefs
within the GBR based upon coral cover monitoring data (Jonker et al.,, 2008).
Through embedding our recovery model within a Bayesian analysis framework, we
provide a powerful tool for uncertainty quantification of key recovery
parameters and coral cover predictions that account for this parameter
uncertainty.

##  2 Materials and methods

In this section we briefly describe the benthic survey data used in this work.
A description of our mathematical model that captures biphasic recovery
patterns is then presented along with the Bayesian methods for model
calibration and assessment of model fitness.

###  2.1 Benthic survey data

Coral cover time-series data are obtained from benthic surveys undertaken by
the Australian Institute of Marine Sciences (AIMS) Long Term Monitoring
Program (LTMP) and Marine Monitoring Program (MMP) between 1992 to 2020.
Combined, these collection programs include taxonomically classified cover
data obtained from over 373 sites across 135 reefs. These monitored reefs are
characteristic of the inshore, middle and outer shelf areas of the GBR (Figure
1(a)). The dataset includes a mixture of annual and biennial surveys at depths
between 2m-9m. At each site, five permanent transects (dimension 5m$\times$50m
for LTMP and 5m$\times$20m for MMP) are photographed at regular intervals (1m
intervals for LTMP and 0.5m intervals for MMP). Coral cover is estimated and
classified using a 5-point stencil overlaid in each image (See Jonker et al.,
(2008) for details on the data collection protocol). Example coral cover time
series are shown for Agincourt 1 (Figure 1(b)), Thetford (Figure 1(c)) and
Lady Musgrave reefs (Figure 1(d)).

Figure 1: (a) Location of reef sites that are monitored as part of the Long-
Term Monitoring Program (LTMP) and the Marine Monitoring Program (MMP). The
MMP monitors the inshore reefs (red), and the LTMP monitors the Middle (green)
and Outer (blue) reef shelves. Example coral cover time-series are shown for
(b) Agincourt 1, (c) Thetford, and (d) Lady Musgrave reefs with their
locations indicated with arrows (black) in (A). Here, total hard coral cover
(black solid), the fast growing Acropridae family (black dotted), and other
hard corals (black dashed) are shown over time along with official disturbance
event records (red dashed). Error bars indicate the standard error of the
coral cover record based on the variation at the transect level.

Following Ortiz et al., (2018), and Warne et al., (2022), disturbances were
identified by statistically significant decline ($p\leq 0.05$) in cover
according to a one-sided paired t-test at the transect level ($n=5$) for
consecutive visits (Figure 1B–D). A recovery trajectory is defined as the
sequence of monitoring surveys between consecutive disturbances inclusive of
the disturbance event defining the start of the recovery. Across the GBR, the
patterns of coral recovery between disturbance events varies widely, as shown
in the examples in Figure 1. This motivates a modelling approach that can
account for deviations from standard growth patterns typically defined using a
Gompertz or logistic model (Thompson et al.,, 2016, 2020; Simpson et al.,,
2023). One limitation of the methodology developed in Warne et al., (2022) was
the restriction to recovery trajectories with low initial cover ($\leq 10\%$)
and have $n\geq 5$ benthic survey visits before the next disturbance. Our new
modelling framework reduces these restrictions, requiring only that the number
of visits between subsequent disturbances is $n\geq 3$ with no restrictions on
the initial coral cover are necessary. From the LTMP and MMP benthic surveys,
we obtain $N=737$ recovery trajectories among $336$ reef sites and $124$ reefs
across the GBR.

Throughout this manuscript, we denote $\mathbf{X}^{(i)}(t_{j})$ as the vector
of taxonomically classified coral cover data at the reef site level for the
$i$-th recovery trajectory, $i=1,2,\ldots,N$, as surveyed at time $t_{j}$, and
$j=0,1,\ldots,n^{(i)}$. Here, $j=0$ relates to the initial visit in which the
reef site was perturbed by a disturbance event and $n^{(i)}$ is the number of
subsequent visits in the $i$-th recovery trajectory. Thus the $i$-th observed
recovery trajectory is
$\mathcal{D}^{(i)}=[\mathbf{X}^{(i)}(t_{0}),\mathbf{X}^{(i)}(t_{1}),\ldots,\mathbf{X}^{(i)}(t_{n^{(i)}})]$
containing the time-series of benthic survey data for that recovery
trajectory.

###  2.2 Model development

Our mathematical model of coral reef recovery accounts for biphasic recovery
patterns as identified in previous work (Warne et al.,, 2022), and allows for
more general single-phase recovery patterns beyond standard approaches such as
logistic or Gompertz growth models (Simpson et al.,, 2022; Thompson and
Dolman,, 2010; Thompson et al.,, 2020). We first describe a single species
version of the model, before extending this to a multi-species model that also
accounts for interactions and competition between coral types.

####  2.2.1 Generalised logistic growth

We begin with the Richards’ growth model (Richards,, 1959; Tjørve and Tjørve,,
2010; Simpson et al.,, 2022), also known as the generalised logistic growth
model, given by

$\underbrace{\frac{\text{d}C(t)}{\text{d}t}}_{\text{Rate of change in
cover}}=\underbrace{\alpha C(t)}_{\text{Coral
growth}}\times\underbrace{\frac{1}{\gamma}\left[1-\left(\frac{C(t)}{K}\right)^{\gamma}\right]}_{\text{Competition
for available space}},\quad t>t_{0}$ (1)

where $K\in(0,100]$ (% of area) is the maximum % of area the coral can
feasibly cover, i.e., the available space, $C(t)\in(0,100]$ (% of area) is the
hard coral cover at time $t$ (years), $\alpha>0$ (1/years) is the intrinsic
rate of cover increase and $\gamma>0$ is a non-dimensional generalisation
parameter that controls the shape of the sigmoid growth curve. For example,
logistic growth is recovered for $\gamma=1$ and Gompertz growth occurs in the
limit as $\gamma\to 0$. The effect of $\gamma$ on the recovery curve is shown
in Figure 2. Note that Richards’ growth model is usually presented without the
$1/\gamma$ factor as it can be absorbed into the definition of the rate
parameter, $\alpha$. However, Simpson et al., (2022) show that this leads to
non-identifiability in the rate parameter and furthermore renders comparison
of $\alpha$ across different values of $\gamma$ to be meaningless.

Figure 2: The effect of the generalisation parameter, $\gamma$, in Richards’
growth model. The logistic growth model corresponds to $\gamma=1$ (blue
solid). For $\gamma\to 0$ we obtain the Gompertz growth model with rapid
initial recovery and stronger effect of spatial competition. Other values for
$\gamma$ result is alternative recovery models. Here, $K=90\%$,
$\alpha=0.5\,(\text{Years}^{-1})$, and $c_{0}=5\%$.

####  2.2.2 Incorporating biphasic recovery

To incorporate biphasic growth into the model, we build on recent developments
by Murphy et al., (2022) who consider a piecewise defined population dynamic
model,

$\frac{\text{d}C(t)}{\text{d}t}=\begin{cases}f_{1}(C(t))&\text{if }t\leq T\\\
f_{2}(C(t))&\text{if }t>T\end{cases},$ (2)

where the rate of change of a population, $C(t)$, switches from $f_{1}(C(t))$
to $f_{2}(C(t))$ at the change point $t=T$. It is important to note that the
solution to Equation (2) is still continuous at this change point despite the
derivative potentially being discontinuous.

To mathematically model the biphasic recovery patterns identified by Warne et
al., (2022), we require $f_{1}(C(t))<f_{2}(C(t))$ that is for the time
$t\in(t_{0},T]$ the coral recovery rate parameter is lower than when
$t\in(T,\infty)$. We achieve this through defining
$f_{1}(C(t))=\alpha_{d}f_{2}(C(t))$ where $\alpha_{d}\in(0,1)$ is the recovery
scale factor, and $f_{2}(C(t))$ corresponds to the Richards’ growth model
(Equation (1)). We then define a duration, $T_{d}>0$, after which the recovery
reverts back to Equation (1). Throughout, we refer to this time $T_{d}$ as the
relative change point, since it is duration of the first phase relative to the
initial time $t_{0}$, that is $T=T_{d}+t_{0}$. Substituting Equation (1) into
the general biphasic model of Murphy et al., (2022) (Equation (2)) leads to
the following biphasic Richards’ growth model,

$\frac{\text{d}C(t)}{\text{d}t}=\begin{cases}\dfrac{\alpha_{d}\alpha}{\gamma}C(t)\left[1-\left(\dfrac{C(t)}{K}\right)^{\gamma}\right]&\text{if
}t_{0}<t\leq t_{0}+T_{d}\\\\[10.00002pt]
\dfrac{\alpha}{\gamma}C(t)\left[1-\left(\dfrac{C(t)}{K}\right)^{\gamma}\right]&\text{if
}t>t_{0}+T_{d}\end{cases}.$ (3)

Given a known initial condition $C(t_{0})=c_{0}$ and values for the parameters
$\boldsymbol{\theta}=[\alpha,\gamma,T_{d},\alpha_{d},K]$, Equation (3) admits
the analytical solution,

$C(t)=\begin{cases}K\left\\{1+[(K/C_{0})^{\gamma}-1]\exp{\left(-\alpha_{d}\alpha(t-t_{0})\right)}\right\\}^{-1/\gamma}&\text{if
}t_{0}<t\leq t_{0}+T_{d}\\\
K\left\\{1+[(K/C_{d})^{\gamma}-1]\exp{\left(-\alpha(t-T_{d})\right)}\right\\}^{-1/\gamma}&\text{if
}t>t_{0}+T_{d}\end{cases},$ (4)

where
$C_{d}=C(t_{0}+T_{d})=K\left\\{1+[(K/C_{0})^{\gamma}-1]\exp{\left(-\alpha_{d}\alpha
T_{d}\right)}\right\\}^{-1/\gamma}$.

####  2.2.3 Accounting for multiple coral family groups

Due to the substantial differences in the recovery rate of the Acroporidae
family of corals, it is common to model their recovery distinctly from other
hard coral families (Thompson and Dolman,, 2010; Thompson et al.,, 2020; Warne
et al.,, 2022). We consider the joint recovery of $M$ coral family groupings
with $C_{m}(t)$ denoting the coral cover of family group $m\in[1,2,\ldots,M]$.
This leads to a coupled system of $M$ biphasic Richards’ growth model
(Equation (3)) with an shared maximum coral cover term to account for
competition for space between species,

$\frac{\text{d}C_{m}(t)}{\text{d}t}=\begin{cases}\dfrac{\alpha_{d,m}\alpha_{m}}{\gamma_{m}}C_{m}(t)\left[1-\left(\dfrac{\sum_{k=1}^{M}C_{k}(t)}{K}\right)^{\gamma_{m}}\right]&\text{if
}t_{0}<t\leq t_{0}+T_{d,m}\\\\[10.00002pt]
\dfrac{\alpha_{m}}{\gamma_{m}}C_{m}(t)\left[1-\left(\dfrac{\sum_{k=1}^{M}C_{k}(t)}{K}\right)^{\gamma_{m}}\right]&\text{if
}t>t_{0}+T_{d,m}\end{cases},$ (5)

where $m=1,2,\ldots,M$ and
$\boldsymbol{\theta}=[\alpha_{1},\gamma_{1},T_{d,1},\alpha_{d,1},\alpha_{2},\gamma_{2},T_{d,2},\alpha_{d,2},\ldots,\alpha_{M},\gamma_{M},T_{d,M},\alpha_{d,M},K]$.
For this work, we only consider the case of $M=1$ or $M=2$. When $M=1$, we
consider only the total hard coral cover, and use the analytical solution
given in Equation (4). When $M=2$, we use $C_{1}(t)$ to denote the coral cover
of the Acroporidae family and $C_{2}(t)$ to denote the coral cover of all
other hard corals. See Section 4, Thompson and Dolman, (2010), and Thompson et
al., (2020) for discussion around this choice of splitting.

The analytic solution (Equation (4)) for the single species model (Equation
(3)) does not extend to the multispecies case (Equation (5)). Therefore, we
rely on numerical schemes to solve the coupled system. We apply the Runge-
Kutta-Fehlberg method with adaptive time-stepping (Fehlberg,, 1969) as
implemented in the R package, deSolve (Soetaert et al.,, 2010).

####  2.2.4 The observation process model

Let $\mathbf{C}(t)=[C_{1},C_{2},\ldots,C_{M}(t)]^{\text{T}}$ be the vector of
coral cover for each of the $M$ modelled family groups. We then model
observation noise using a Gaussian distribution, that is the observed cover,
$\mathbf{X}(t)$, at time $t$ is given by

$\mathbf{X}(t)\sim\mathcal{N}(\mathbf{C}(t),\boldsymbol{\Sigma}(t)),$ (6)

where
$\boldsymbol{\Sigma}(t)=\operatorname*{diag}{\\{\hat{s}_{m}(t)^{2}/5\\}}$ is a
diagonal matrix of standard error estimates where $\hat{s}_{m}(t)^{2}$ is the
empirical variance in coral cover for family group $m$ at time $t$ taken
across the five transects (See Section 2.1). Alternative observation processes
could be considered, such as binomial noise, however, the observed
heteroscedasticity in the cover data tends to be different than expected for
multinomial noise models, that is, the data demonstrates both overdispersion
and underdispersion. Using a Gaussian noise approximation allows for a more
flexible treatment of the variability that is also computationally efficient.

###  2.3 Bayesian uncertainty quantification

Inference is performed within a Bayesian framework, that is, for each recovery
trajectory $i\in[1,2,\ldots,N]$ we require samples from the joint posterior
distribution (Gelman et al.,, 2013) with density,

$p(\boldsymbol{\theta}^{(i)}\mid\mathcal{D}^{(i)})=\frac{\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})p(\boldsymbol{\theta}^{(i)})}{p(\mathcal{D}^{(i)})},$

where $\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})$ is the
likelihood of the observed benthic survey data, $\mathcal{D}^{(i)}$, under the
model given parameters, $\boldsymbol{\theta}^{(i)}$,
$p(\boldsymbol{\theta}^{(i)})$ is the prior probability density and
$p(\mathcal{D}^{(i)})$ is the evidence, that is, the marginal likelihood given
by

$p(\mathcal{D}^{(i)})=\int\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})p(\boldsymbol{\theta}^{(i)})\textrm{d}\boldsymbol{\theta}^{(i)}.$

Conditional on the model solution
$\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})$ (which will be numerical for
$M>1$) we obtain the log-likelihood as a multivariate Gaussian around the
model,

$\log\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})=-\frac{1}{2}\left\\{n^{(i)}\log(2\pi)+\sum_{j=0}^{n^{(i)}}\left[\log{|\boldsymbol{\Sigma}^{(i)}(t_{j})|}+\Delta\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})^{\text{T}}\boldsymbol{\Sigma}^{(i)}(t_{j})^{-1}\Delta\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})\right]\right\\},$
(7)

where
$\Delta\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})=\mathbf{X}^{(i)}(t_{j})-\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})$
is the difference between the observations and modelled mean recovery at time
$t_{j}$ with $j=0,1,\ldots,n^{(i)}$.

We assume independent priors, that is,

$p(\boldsymbol{\theta}^{(i)})=\prod_{m=1}^{M}p(\alpha_{m}^{(i)})p(\gamma_{m}^{(i)})p(T_{d,m}^{(i)})p(\alpha_{d,m}^{(i)}).$
(8)

Note that the maximum coral cover, $K^{(i)}$, is not inferred since it is
standard practice to estimate this directly from the LMTP and MMP Data
(Thompson et al.,, 2020) (See Section 4 for discussion). For the remaining
parameters we have, $\alpha_{m}^{(i)}\sim\mathcal{U}(0,1)$,
$\alpha_{d,m}^{(i)}\sim\mathcal{U}(0,0.9)$,
$\gamma_{m}^{(i)}\sim\mathcal{U}(0,2)$,
$T_{d,m}^{(i)}\sim\mathcal{U}(0,D^{(i)})$ where
$D^{(i)}=t_{n}^{(i)}-t_{0}^{(i)}$ is the duration of the $i$th recovery
trajectory and $\mathcal{U}(a,b)$ denotes a uniform distribution over the
interval $[a,b]$. By using uniform priors, we aim for the priors to be vague
or uninformative within the physically feasible ranges, however, we note that
general guarantees of uninformativity are not practically available (Efron,,
1986; Warne et al.,, 2019).

For our model, the evidence term is intractable. As a result, we can only
evaluate the posterior density up to a normalising constant using Equation (7)
and Equation (8). To deal with this, we apply Markov chain Monte Carlo (MCMC).
In our case, the commonly applied sampler available in the Stan environment
(Carpenter et al.,, 2017) or other samplers based on Hamiltonian Monte Carlo
(Duane et al.,, 1987) are not appropriate due to a discontinuity that arises
in $\nabla\log\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})$ from
the change-points $T_{d,m}$. While numerical smoothing of the discontinuities
is feasible, we avoid this by selecting a more general MCMC sampler that does
not require $\log\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})$ to
be differentiable. Specifically, we use the robust adaptive random-walk
Metropolis-Hastings MCMC scheme of Vihola, (2012) that is implemented within
the R package, adaptMCMC (Scheidegger,, 2021).

For each trajectory,
$\mathcal{D}^{(1)},\mathcal{D}^{(2)},\ldots,\mathcal{D}^{(N)}$, we use four
independent chains to assess standard convergence diagnostics of effective
sample size (ESS) and the Gelman $\hat{R}$ statistic (Gelman and Rubin,, 1992;
Vehtari et al.,, 2021) using the coda package (Plummer et al.,, 2006).
Sampling is terminated when $\hat{R}\leq 1.1$ and $\text{ESS}\geq 200$ for
each parameter. This results in $N=737$ sets of posterior samples. Thinning is
applied for chains with very large ESS to save memory requirements while still
ensuring the convergence criteria are satisfied.

To assess the model fit, we sample from the posterior predictive distribution
with density,

$p(\mathcal{D}_{s}^{(i)}\mid\mathcal{D}^{(i)})=\int
p(\mathcal{D}_{s}^{(i)}\mid\boldsymbol{\theta}^{(i)})p(\boldsymbol{\theta}^{(i)}\mid\mathcal{D}^{(i)})\,\text{d}\boldsymbol{\theta}^{(i)},$
(9)

where $\mathcal{D}_{s}^{(i)}$ is simulated trajectory data using the biphasic
recovery model (Equation (5)) and observation process (Equation (6)) given a
random draw from the posterior samples as obtained via MCMC. Using these
simulations, we obtain a probability density estimate for the model
uncertainty accounting for all of the underlying parameter uncertainty. This
density is then used to estimate 50%, 90%, 95%, and 99% credible intervals
(CrI). This enables the model fit to be assessed by inspecting the location of
the observed data within these reported intervals.

##  3 Results

We calibrate the biphasic generalised recovery model using Bayesian analysis
for the $N=737$ recovery trajectories as described Section 2. The resulting
output is posterior distributions and posterior predictive samples for each of
the $N=737$ recovery trajectories. In this section, we provide summaries of
the overall model fitness in terms of convergence diagnostics and posterior
predictive probabilities. We also present example outputs using four
representative recovery trajectories (Figure 3): site 2 of Gannett Cay reef
between 2001–2009 (Figure 3(a)), site 1 of Lady Musgrave reef between
1992–2003 (Figure 3(b)), and site 1 of Thetford reef between 1994–1999 (Figure
3(c)) and between 2001–2009 (Figure 3(d)). For ease of description we refer to
these recovery trajectories by their reef name and start year, that is Gannett
Cay 2001 (Figure 3(a)), Lady Musgrave 1992 (Figure 3(b)), Thetford 1994
(Figure 3(c)), and Thetford 2001 (Figure 3(d)).

Figure 3: Example recovery trajectory data used to highlight the efficacy of
our model in capturing different recovery patterns observed for hard coral
(solid black lines) broken down into the Acroporidae (dotted black lines) and
all other hard corals (dashed black lines). Error bars indicate standard
errors for the coral cover estimates.

These example recovery trajectories are selected to demonstrate the range of
bisphasic recovery patterns that are observed across the GBR while also
providing interesting interplay between the single species and two species
models.

###  3.1 MCMC convergence diagnostics

For each recovery trajectory we obtain samples from the joint posterior
distribution under both the single-species (Equation (3)) and the two-species
(Equation (5)) versions of the biphasic Richards’ growth model. For the single
species model we infer the intrinsic growth rate $\alpha$, the shape
parameter, $\gamma$, the relative change-point, $T_{d}$, and the scale factor,
$\alpha_{d}$. For the two-species model we infer equivalent parameters for the
Acroporidae family, $\alpha_{A},\gamma_{A},T_{d,A},\alpha_{d,A}$, for all
other hard corals, $\alpha_{C},\gamma_{C},T_{d,C},\alpha_{d,C}$. In both
cases, we used standard methods to estimate maximum coral cover for reefs
based on abiotic cover and transient silt records (Jonker et al.,, 2008;
Thompson and Dolman,, 2010; Thompson et al.,, 2020).

As described in Section 2.3 we terminate sampling when four independent chains
achieve $\widehat{R}\leq 1.1$ and ESS $\geq 200$. Of the $N=737$ recovery
trajectories, $38$ did not meet this convergence criteria after 200,000 MCMC
iterations. Most of these failed convergences where due to ESS $<200$, with
only 11 failing due to $\hat{R}>1.1$. The lower ESS indicates slower mixing in
the MCMC chains and a likely cause is the parameter non-identifiability issues
that arises when there a very weak biphasic recovery pattern leading to a
trade-off between $T_{d}$, $\alpha_{d}$ and $\gamma$, that is the data does
not have enough information to distinguish between a single-phase model with
$T_{d}=0$ and a biphasic model $T_{d}>0$. In the very few cases where
$\hat{R}>1.1$, it seems that the proposed model is a poor fit, and we discuss
potential improvements to account for these cases in the discussion (Section
4). That is, we are able to reliably characterise model parameters for $95\%$
of the available recovery trajectories ($N=699$) with only $1.5\%$ ($N=11$)
failing due to serious convergence issues arising from poor model fit.

Figure 4: Scatter plots summarising the MCMC convergence diagnostics for each recovery trajectory after a maximum of 200,000 iterations per chain and four independent chains. Convergence thresholds are indicated: $\hat{R}\leq 1.1$ (vertical red dashed line), $\text{ESS}\geq 200$ (horizontal red dashed line). (a) For the majority of recovery trajectories, the MCMC diagnostics comfortably meet the thresholds. (b) Of the small number of recovery trajectories that fail in their MCMC convergence diagnostics, most are close to the threshold.  Figure 5: Trace plots showing good mixing and agreement between four independent MCMC chains (green, orange, purple, and magenta solid lines) in all parameters. (a) Gannett Cay 2001, (b) Lady Musgrave 1992, (c) Thetford 1994, and (d) Thetford 2001. Table 1: Convergence diagnostic results for the four example recovery trajectories. Recovery Trajectory | Parameter |  $\hat{R}$ [Upper $95\%$CI] | ESS
---|---|---|---
| $\alpha$ |  $1.02$ $[1.06]$ | 429
Gannett Cay site 2 2001 | $\alpha_{d}$ |  $1.01$ $[1.02]$ | 1092
| $\gamma$ |  $1.01$ $[1.04]$ | 990
| $T_{d}$ |  $1.02$ $[1.05]$ | 227
| $\alpha$ |  $1.00$ $[1.00]$ | 1074
Lady Musgrave site 1 1992 | $\alpha_{d}$ |  $1.00$ $[1.01]$ | 707
| $\gamma$ |  $1.00$ $[1.01]$ | 446
| $T_{d}$ |  $1.00$ $[1.01]$ | 454
| $\alpha$ |  $1.00$ $[1.01]$ | 969
Thetford site 1 1994 | $\alpha_{d}$ |  $1.08$ $[1.15]$ | 287
| $\gamma$ |  $1.00$ $[1.02]$ | 1549
| $T_{d}$ |  $1.01$ $[1.03]$ | 589
| $\alpha$ |  $1.01$ $[1.01]$ | 576
Thetford site 1 2001 | $\alpha_{d}$ |  $1.01$ $[1.03]$ | 395
| $\gamma$ |  $1.01$ $[1.04]$ | 365
| $T_{d}$ |  $1.01$ $[1.03]$ | 402

In the context of the four examples trajectories (Figure 3), the convergence
diagnostics in the form of MCMC trace plots, parameter-wise ESS and
$\widehat{R}$ statistics are shown in Figure 5 and Table 1 for the single
species model. The trace plots show good mixing and agreement across all four
independent MCMC chains (Figure 5). In addition the estimated $\widehat{R}$
statistics including the upper $95\%$CI are well below the standard threshold
of $1.1$ (with exception of the $\widehat{R}$ upper $95\%$CI for $\alpha_{d}$
at Thetford 1994), and the ESS for all parameters exceed 200. Given these
results, it is reasonable to conclude our convergence criteria are appropriate
and that the Markov chains have reached stationarity.

###  3.2 Assessment of model fitness

For each of the $N=699$ recovery trajectories in which the MCMC sampling
satisfied convergence criteria (Section 3.1), we perform posterior predictive
checks. That is, we evaluate the location of the observed coral cover within
the posterior predictive distribution. For the model to be representing the
data accurately we would expect, for any $0<\beta<100$, that at least
$\beta\%$ of the observations are within the $\beta\%$CrI of the relevant
posterior predictive distribution.

For each observation, we calculate the smallest $\beta$ such that the
observation is within $\beta\%$CrI. To do this, samples are generated for the
posterior predictive distribution (Equation (9)) and the observed quantiles
for each coral cover observation is calculated using,

$Q_{\textrm{obs}}^{i,j}=\mathbb{P}(\mathbf{X}_{s}^{(i)}(t_{j})\leq\mathbf{X}^{(i)}(t_{j})\mid\mathcal{D}^{(i)})=\int_{0}^{\mathbf{X}^{(i)}(t_{j})}p(\mathbf{X}_{s}^{(i)}(t_{j})\mid\mathcal{D}^{(i)})\,\textrm{d}\mathbf{X}_{s}^{(i)}(t_{j}),$
(10)

where $\mathbf{X}_{s}^{(i)}(t_{j})$ and $\mathbf{X}^{(i)}(t_{j})$ are,
respectively, the simulated and observed coral cover for the $i$th recovery
trajectory at the $j$th time point. Given the observed quantile, then the
smallest $\beta$ for the CrI containing the observation is given by
$\beta^{i,j}=100(1-2Q_{\textrm{obs}}^{i,j})$ when
$Q_{\textrm{obs}}^{i,j}<0.5$, and $\beta^{i,j}=100(2Q_{\textrm{obs}}^{i,j}-1)$
otherwise. Hence, we can estimate the proportion of observations that within
the posterior predictive $\beta\%$CrI using

$\hat{p}_{\beta}=\frac{\sum_{i=1}^{N}\sum_{j=1}^{n^{(i)}}\mathbb{I}\left(\beta^{i,j}<\beta\right)}{N\sum_{i=1}^{N}n^{(i)}}.$
(11)

Here, $\mathbb{I}\left(\beta^{i,j}<\beta\right)=1$ if $\beta^{i,j}<\beta$ and
$\mathbb{I}\left(\beta^{i,j}<\beta\right)=0$ otherwise. The standard error for
the estimate of is given by

$\hat{s}_{\beta}=\sqrt{\frac{\hat{p}_{\beta}(1-\hat{p}_{\beta})}{N\sum_{i=1}^{N}n^{(i)}}}.$
(12)

Comparing estimates of $\hat{p}_{\beta}$ with $\beta=1,2,\ldots 99$ provides a
summary of how well posterior predictive distributions fit the data. As
demonstrated in Figure 6(a), we observe that for $\beta\in\\{1,91\\}$, we have
$\beta<\hat{p}_{\beta}\pm 1.96\hat{s}_{\beta}$, indicating the data is well
represented by the bulk of the posterior predictive distributions. When
$\beta\in\left\\{92,93,\ldots,97\right\\}$ we observe
$\hat{p}_{\beta}-1.96\hat{s}_{\beta}<\beta<\hat{p}_{\beta}+1.96\hat{s}_{\beta}$,
indicating the posterior predictive distributions are also tracking the more
extreme data points at the expected frequencies. In the extreme tails,
$\beta\in\left\\{98,99\right\\}$ we have $\beta>\hat{p}_{\beta}\pm
1.96\hat{s}_{\beta}$ and conclude the very extreme variations in the data are
under-represented, through the margin of inaccuracy here is still only a few
percent. We conclude that overall the performance of the biphasic Richards’
model in capturing the variation and dynamics of coral reef recovery is
extremely high.

Figure 6: Posterior predictive analysis over the entire GBR data set including
$N=699$ recovery trajectories for which the MCMC sampling converged according
to our criteria. (a) The estimated proportion of data that falls within the
$\beta\%$CrI, $\hat{p}_{\beta}$, (solid black lines, with error bars for
$95\%$ confidence intervals) for each level of $\beta$ (red line for
$\beta=\hat{p}_{\beta}$). (b) Bivariate density plot of the smallest
containing posterior predictive CrI for each coral cover observation (yellow
indicates high density and dark blue low density).

In Figure 6(b) a density plot shows the distribution of smallest containing
CrIs for each coral cover observation. For coral covers below $25\%$ cover the
pattern, the bulk of smallest containing CrIs are well below the $25\%$CrI and
the distribution skewed to the right, indicating high posterior probability
density for most of the observed data. As coral cover increases, especially
above $40\%$ cover, the distribution of smallest containing CrIs is more
diffuse. From this we conclude that the majority of the data overdispersion is
occurring in regions of higher cover. Importantly, the greatest accuracy in
the model is for cover below $30\%$, which is where the majority of the GBR
sits. It is also important to emphasise that at higher cover we still observe
good fitness, but the data is more variable in the extremes than the model
predicts.

Figure 7(a)–(d) demonstrates example posterior predictive distributions for
both the single species model (Equation (3)), and the two species model
(Equation (5)) that describes the competition for space between fast and slow
growing corals. The posterior predictive plots are compared with the observed
data in both the single-species and two-species representations for each of
the example recovery trajectories: Gannett Cay 2001 (Figure 7(a)), Lady
Musgrave 1992 (Figure 7(b)), Thetford 1994 (Figure 7(c)), and Thetford 2001
(Figure 7(d)). In each case, there is excellent agreement between the model
predictions and the observed data with majority of data points lying within
the 95%CrI of the posterior predictive distribution.

Figure 7: Posterior predictive distributions for example recovery
trajectories. (a)–(d) The observed data is given for the total cover for hard
coral (solid black lines) broken down into Acroporidae (dotted black lines),
and other hard corals (dashed black lines); and model predictions with
credible intervals (50%, 90%, 95% and 99% CrI) as shaded regions for total
hard coral (green), Acroporidae (red), and other hard corals (blue).

###  3.3 Parameter inference

For a given recovery trajectory, marginal posterior densities for the model
parameters are estimated directly from the MCMC samples obtained by targeting
the full posterior distribution (Sections 2.3 and 3.1). These marginal
densities can be used to obtain insight into the underlying mechanisms leading
to a particular recovery pattern for a recovery trajectory of interest. In
Figure 8, univariate posterior marginal probability densities are shown by
coral family group for each of the examples recovery trajectories: Gannett Cay
2001 (Figure 8(a),(e),(i),(m)), Lady Musgrave 1992 (Figure 8(b),(f),(j),(n)),
Thetford 1994 (Figure 8(c),(g),(k),(o)), and Thetford 1994 (Figure
8(d),(h),(l),(p)). Among these examples, our model identifies a range of
different recovery patterns that could be of interest to reef managers.

Across each of the example trajectories there are some common features in the
marginal posterior distributions for the recovery rates (Figure 8(a)-(d)).
Firstly, the two species model tends to indicate the recovery rate for the
Acroporidae family is greater than the other hard corals, and the recovery
rate for the single species model tends to sit between the two sub-groups with
the relative weighting between the two subgroups indicative of their relative
coral cover contribution (Figure 7 and Figure 8(a)-(d)). This is expected
since members of the Acroporidae family are known to be fast growing by
comparison to most other hard corals (Thompson and Dolman,, 2010; Ortiz et
al.,, 2021). An interesting deviation from this pattern is in the estimates
for the recovery rate for Thetford 2001 (Figure 8(d)) where the total hard
coral estimates are more dispersed than either of the sub-groups. An
explanation for this could be the more complex dynamics in the other hard
corals with an initial decline before an increase in cover coinciding with the
Acroporidae corals exceeding the cover of the other hard corals (Figure 7(d)).
This results in a more complex recovery pattern for the total hard coral
cover.

The shape parameter in the Richards’ growth model shows a variety of patterns
across the example recovery trajectories (Figure 8(e)–(h)) with some examples
approximating Gompertz growth ($\gamma\to 0$), others approximating logistic
growth ($\gamma=1$), and others tending toward other forms from Richards’
growth ($\gamma>1$). This is not surprising given the diversity in coral
species in the GBR. It is worth noting that the posterior distributions for
the shape parameter in both Thetford 1994 (Figure 8(g)) and Thetford 2001
(Figure 8(h)) are almost identical for the total and other hard corals, which
would be expected for the same reef. The difference observed in the posterior
for the shape parameter for the Acroporidae is likely due to the data
containing little information on the sigmoid shape since the recovery for
Acroporidae is near linear (Figure 7(d)). An alternate explanation could be a
change in the specific community composition in the Acroporidae group at
Thetford, we do not explore this further here. We highlight, however, that
identifying similarities and differences in the recovery patterns of a
specific reef over time is a key use case for coral recovery models in marine
ecology and monitoring (Thompson et al.,, 2020; Ortiz et al.,, 2018).

Figure 8: Univariate marginal posterior densities for (a)–(d) the recovery
rate $\alpha$, (e)–(h) shape parameter $\gamma$, (i)–(l) change-time $T_{d}$,
and (m)–(p) scale factor $\alpha_{d}$. Densities are given for total hard
coral (green) broken down into Acroporidae (red), and other hard corals
(blue). Results are shown for the four examples recovery trajectories:
(a),(e),(i),(m) Gannett Cay 2001, (b),(f),(j),(n) Lady Musgrave 1992,
(c),(g),(k),(o) Thetford 1994, and (d),(h),(l),(p) Thetford 2001.

For the parameters that relate to the biphasic recovery, the example recovery
trajectories demonstrate the range of different behaviours between different
locations and times, and between the single species and two species
representations. In the following paragraphs we highlight key features in each
example. One important property of our biphasic model (Equation (3) and
Equation (5)) is that as the scale factor $\alpha_{d}\to 1$, the change-time
becomes structurally non-identifiable. Thus whenever the posterior density for
the change-time is very flat, and the posterior density for the scale factor
skewed to the left with bulk above $0.6$, then we conclude there is no
evidence for biphasic recovery.

In the analysis of Gannett Cay 2001, the single species model indicates little
evidence of a biphasic recovery pattern for total hard coral cover. This
conclusion is drawn from the relatively flat posterior for the change time in
the (Figure 8(i)), and the bulk of the scale factor is tending toward larger
values (Figure 8(m)). However, when the data is interpreted with the two
species model, there is a strong signal for biphasic recovery patterns are
present in the Acroporidae group with strong posterior support for change-time
between $1$ to $4$ years (Figure 8(i)) and a small value for the scale factor
(Figure 8(m)). For the other hard corals, the posteriors for the change-time
and scale factor are very similar to the single species model with little
evidence for biphasic recovery.

At Lady Musgrave 1992, there is a strong indication for biphasic recovery
patterns in total hard coral cover with a change-time of between $1$ to $3$
years (Figure 8(j)) and scale factor less than $0.6$ (Figure 8(n)). When using
the two species model, it is clear that a biphasic recovery pattern in the
Acroporidae family, with a change-time of around $1.5$ to $2.5$ years and
scale factor less that $0.5$, provides the main driver for the overall
biphasic recovery for the total hard coral cover. This is expected at Lady
Musgrave where the Acroporidae family accounts for almost 100% of corals
(Figure 7(b)) (Simpson et al.,, 2022; Murphy et al.,, 2022; Simpson et al.,,
2023). However, the uncertainty in the biphasic parameters for the small
population of other hard corals in the two species model (Figure 7(b)) is
reflected in the larger variance in the change-time and scale factor posterior
densities for the total hard corals (Figure 8(j),(n)).

Analysis for Thetford 1994 and 2001, are consistent with each other in
relation to biphasic recovery. In both recovery trajectories, there is
evidence for biphasic recovery patterns in the total hard coral cover (Figure
8(k)–(l),(o)–(p)) (though there is substantially more uncertainty in the
results following 1994). Just as in the analysis of Lady Musgrave 1992, one of
the species in the two species model seems to be the main driver for biphasic
recovery pattern in the total hard coral cover, however, for Thetford it is
the other hard corals rather than the Acroporidae. We arrive at this
conclusion by noting, the change-time and scale factor posteriors densities
for the other hard corals follow qualitatively a similar pattern to the
equivalent posteriors for total hard corals (Figure 8(k)–(l),(o)–(p)). The
Thetford 2001 recovery trajectory is more informative in relation to the
biphasic recovery pattern (Figure 8(k),(o)), however, both Thetford 1994 and
2001 are qualitatively similar with no evidence for biphasic recovery patterns
in the Acroporidae family and agreement between other hard corals and total
hard corals in support of biphasic recovery patterns (Figure
8(k)–(l),(o)–(p)). Due to the consistency between distinct trajectories, this
would support the theory that this biphasic recovery pattern is a feature of
the environment or ecosystem at Thetford reef.

##  4 Discussion

In this work, we perform Bayesian parameter inference and posterior predictive
checks for $N=699$ recovery trajectories across the GBR. We find both the
single species and two species biphasic Richards’ growth model is flexible
enough to obtain good model posterior predictive performance. This is
impressive in of itself given the range of possible recovery patterns that
could arise from the complex interaction of different benthic species. Using
four recovery trajectories as exemplars, we highlight the utility of the
method for uncertainty quantification in key model parameters related to the
Richards’ growth shape, recovery rate, and biphasic recovery patterns.
Interpretation of these parameters both from a single species and two species
viewpoint could support reef management and interventions following
significant disturbance events.

Accurately predicting and interpreting recovery patterns for diverse locations
across the GBR is critical for management to support reef recovery, resistance
and adaptation into the future (Ortiz et al.,, 2021; Thompson et al.,, 2020).
Motivated by prior work demonstrating potential widespread deviations from
classical Gompertz or logistic recovery models (Warne et al.,, 2022; Murphy et
al.,, 2022; Ortiz et al.,, 2018), we develop a biphasic extension to the
multispecies Richards’ growth model. This framework is flexible enough to
capture a range of possible sigmoid recovery patterns and biphasic recovery
patterns where a period of inhibited recovery occurs in the year immediately
following the disturbance event (Warne et al.,, 2022). In this paper, we
formulate this new model and demonstrate its efficacy in accurately capturing
the recovery patterns of nearly every recovery trajectory available within
LTMP reef monitoring data, which is spatially and temporally the largest reef
monitoring dataset in the world. Therefore, our modelling approach is
describing well the recovery patterns observed in the largest and most
interconnected coral reef worldwide. Through Bayesian uncertainty
quantification we can reliably characterise the combination of parameters that
lead to particular recovery behaviours for a reef and provide forecasts of
coral recovery following a disturbance event.

In this work, we considered the LTMP data from a single species perspective,
based on total hard coral cover, and from a two species perspective that
separates the fast-growing Acroporidae coral families from the other hard
corals (Thompson and Dolman,, 2010; Thompson et al.,, 2020). We believe that
both models are important for management. Firstly, the single species model
enables a macroscopic view of the reef recovery patterns that is useful as a
reef health indicator (Thompson et al.,, 2020) and the notion of a biphasic
recovery has a direct interpretation in terms of expected future coral cover
that, in turn, can support intervention triage (Warne et al.,, 2022). In
addition, the single species model is extremely efficient for the purposes of
model calibration due to the lower dimensionality of the parameter space and
the analytic solution to the biphasic Richards’ growth model (Tjørve and
Tjørve,, 2010; Simpson et al.,, 2022). However, the two species model provides
more insight around the potential mechanisms driving any biphasic recoveries
that are observed using the single species model. This provides additional
information to management authorities regarding the types of interventions
that may be applicable. Unfortunately, this additional granularity of
interpretation comes with a computational cost as numerical schemes are needed
to solve the governing equations and there are more parameters that can lead
to slower convergence in the MCMC sampling. Therefore, we recommend a two-
stage process for monitoring in which the single species model is used to
triage reefs affected by biphasic recovery, then the two species model is
introduce to assist in planning interventions. Future work could also explore
the use of frequentist methods based on likelihood profiles (Simpson et al.,,
2023) as a preconditioning step to accelerate the MCMC sampling and reduce the
computational burden as proposed by Warne et al., (2024). Such improvements
would be important if further refinements are considered, such as the
inclusion of additional hard coral family divisions or the inclusion of other
benthic life forms that compete for space with coral such as macro-algae and
sponges (Thompson and Dolman,, 2010; Mumby,, 2009), or the inclusion of
spatial hierarchical patterns (Johns et al.,, 2014).

One limitation to our process is the pre-estimation of maximum coral cover.
While this is a standard procedure in the coral monitoring literature
(Thompson et al.,, 2016, 2020), there is evidence to suggest that these
estimates could be overestimates to the effective maximum coral cover
(Cresswell et al.,, 2024). In addition, it is also mathematically problematic
to pre-estimate long-term phenomena, in general. We also note that this pre-
estimation of the maximum coral cover is one of the main reasons for the
relatively few situations in which we obtain poor model fits or convergence
issues with MCMC sampling. Therefore, future work should also consider the
inclusion maximum coral cover parameter in the model calibration process. This
will be especially important if there is latent competition with other
species, like macro-algae, that are not currently explicitly modelled.

We have only considered deterministic models designed to capture the average
recovery patterns. As a result, the noise is assumed to be entirely due to the
monitoring process (Jonker et al.,, 2008; Thompson et al.,, 2016). However, at
a colony level, reef growth and competition with other species is a stochastic
process (Álvarez–Noriega et al.,, 2023; Bozec et al.,, 2021; McDonald et al.,,
2023). The net effect of this intrinsic noise could be captured at the reef
scale using stochastic differential equation (SDE) models. While the use of
SDEs generally requires more complex computational approaches to perform
parameter inference (Cranmer et al.,, 2020; Drovandi and Pettitt,, 2010;
Sisson et al.,, 2018; Wang et al.,, 2024; Warne et al.,, 2020, 2024), there
are known benefits in using SDE models when intrinsic noise is relevant, as
they lead to more accurate parameter inferences and can resolve issue related
to parameter non-identifiability (Browning et al.,, 2020; Munsky et al.,,
2009). This could be of considerable benefit, as practical identifiability
issues arise in the uncertainty quantification analysis for some parameters
(Figure 8) and is it not feasible to collect more data.

Prior work by Warne et al., (2022) demonstrated the importance of taking
biphasic recovery patterns in coral reefs into account. Through the
application of our Bayesian framework that accounts for this phenomena, the
construction of coral reef health indicators that are robust to deviations in
standard assumptions could be constructed (Simpson et al.,, 2022; Thompson and
Dolman,, 2010; Thompson et al.,, 2020). This new approach enables a more
nuanced view of expected reef recovery and provides a modelling basis to
investigate the biological, ecological and environmental factors that
contribute to biphasic recovery patterns (Warne et al.,, 2022). New measures
of reef decline can also be established with this approach. For example,
tracking the distribution of the change time and scale factor for a given reef
over multiple recoveries can identify steady increases in the duration of the
slower initial recovery phase that could indicate degradation in reef health.
We anticipate our modelling and analysis approach to have direct benefit for
the management of the GBR. Finally, through the use of decision-theoretic
approaches to sampling design (Akinlotan et al.,, 2024), understanding
biphasic recovery patterns could enable new monitoring protocols that maximise
information for decision makers.

##### Acknowledgements

This work was supported through the Australian Institute of Marine Science
(AIMS) and the Centre for Data Science at the Queensland University of
Technology. Computational resources were provided by the eResearch Office,
Queensland University of Technology.

## References

  * Akinlotan et al., (2024) Akinlotan, M. D., Warne, D. J., Helmstedt, K. J., Vollert, S. A., Chadès, I., Heneghan, R. F., Xiao, H., and Adams, M. P. (2024).  Beyond expected values: Making environmental decisions using value of information analysis when measurement outcome matters.  Ecological Indicators, 160:111828.
  * Aminikhanghahi and Cook, (2016) Aminikhanghahi, S. and Cook, D. J. (2016).  A survey of methods for time series change point detection.  Knowledge and Information Systems, 51(2):339–367.
  * Bodnar et al., (2013) Bodnar, M., Piotrowska, M. J., and Foryś, U. (2013).  Gompertz model with delays and treatment: Mathematical analysis.  Mathematical Biosciences and Engineering, 10:551.
  * Bozec et al., (2021) Bozec, Y., Hock, K., Mason, R. A. B., Baird, M. E., Castro‐Sanguino, C., Condie, S. A., Puotinen, M., Thompson, A., and Mumby, P. J. (2021).  Cumulative impacts across Australia’s great barrier reef: a mechanistic evaluation.  Ecological Monographs, 92(1):e01494.
  * Brouwer et al., (2017) Brouwer, A. F., Eisenberg, M. C., Remais, J. V., Collender, P. A., Meza, R., and Eisenberg, J. N. S. (2017).  Modeling biphasic environmental decay of pathogens and implications for risk analysis.  Environmental Science & Technology, 51(4):2186–2196.
  * Browning et al., (2020) Browning, A. P., Warne, D. J., Burrage, K., Baker, R. E., and Simpson, M. J. (2020).  Identifiability analysis for stochastic differential equation models in systems biology.  Journal of The Royal Society Interface, 17(173):20200652.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A Probabilistic Programming Language.  Journal of Statistical Software, 76(1):1–32.
  * Castro-Sanguino et al., (2021) Castro-Sanguino, C., Ortiz, J. C., Thompson, A., Wolff, N. H., Ferrari, R., Robson, B., Magno-Canto, M. M., Puotinen, M., Fabricius, K. E., and Uthicke, S. (2021).  Reef state and performance as indicators of cumulative impacts on coral reefs.  Ecological Indicators, 123:107335.
  * Cheal et al., (2017) Cheal, A. J., MacNeil, M. A., Emslie, M. J., and Sweatman, H. (2017).  The threat to coral reefs from more intense cyclones under climate change.  Global Change Biology, 23(4):1511–1524.
  * Cranmer et al., (2020) Cranmer, K., Brehmer, J., and Louppe, G. (2020).  The frontier of simulation-based inference.  Proceedings of the National Academy of Sciences, 117(48):30055–30062.
  * Cresswell et al., (2024) Cresswell, A. K., Haller-Bull, V., Gonzalez-Rivero, M., Gilmour, J. P., Bozec, Y.-M., Barneche, D. R., Robson, B., Anthony, K., Doropoulos, C., Roelfsema, C., Lyons, M., Mumby, P. J., Condie, S., Lago, V., and Ortiz, J.-C. (2024).  Reproducing within-reef variability in coral dynamics with a metacommunity modelling framework.  bioRxiv.org.
  * de A. Møller et al., (2021) de A. Møller, C., Christensen, B., and Rattray, F. (2021).  Modelling the biphasic growth of non-starter lactic acid bacteria on starter-lysate as a substrate.  International Journal of Food Microbiology, 337:108937.
  * De’ath et al., (2012) De’ath, G., Fabricius, K. E., Sweatman, H., and Puotinen, M. (2012).  The 27-year decline of coral cover on the Great Barrier Reef and its causes.  Proceedings of the National Academy of Sciences, 109(44):17995–17999.
  * De’ath et al., (2009) De’ath, G., Lough, J. M., and Fabricius, K. E. (2009).  Declining Coral Calcification on the Great Barrier Reef.  Science, 323(5910):116–119.
  * Drovandi and Pettitt, (2010) Drovandi, C. C. and Pettitt, A. N. (2010).  Estimation of parameters for macroparasite population evolution using approximate bayesian computation.  Biometrics, 67(1):225–233.
  * Duane et al., (1987) Duane, S., Kennedy, A., Pendleton, B. J., and Roweth, D. (1987).  Hybrid Monte Carlo.  Physics Letters B, 195(2):216–222.
  * Duran et al., (2017) Duran, A., Shantz, A. A., Burkepile, D. E., Collado-Vides, L., Ferrer, V. M., Palma, L., Ramos, A., and Gonzalez-Díaz, S. P. (2017).  Fishing, pollution, climate change, and the long-term decline of coral reefs off Havana, Cuba.  Bulletin of Marine Science, 94(2):213–228.
  * Efron, (1986) Efron, B. (1986).  Why isn’t everyone a Bayesian?  The American Statistician, 40(1):1–5.
  * Fabricius et al., (2008) Fabricius, K. E., De’ath, G., Puotinen, M. L., Done, T., Cooper, T. F., and Burgess, S. C. (2008).  Disturbance gradients on inshore and offshore coral reefs caused by a severe tropical cyclone.  Limnology and Oceanography, 53(2):690–704.
  * Fehlberg, (1969) Fehlberg, E. (1969).  Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems.  Technical Report 315, National aeronautics and space administration.
  * Gelman et al., (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, Aki, and Rubin, D. B. (2013).  Bayesian Data Analysis.  Chapman and Hall/CRC, 3rd edition.
  * Gelman and Rubin, (1992) Gelman, A. and Rubin, D. B. (1992).  Inference from Iterative Simulation Using Multiple Sequences.  Statistical Science, 7(4):457–472.
  * Hinkley, (1970) Hinkley, D. V. (1970).  Inference about the change-point in a sequence of random variables.  Biometrika, 57(1):1–17.
  * Honsey et al., (2016) Honsey, A. E., Staples, D. F., and Venturelli, P. A. (2016).  Accurate estimates of age at maturity from the growth trajectories of fishes and other ectotherms.  Ecological Applications, 27(1):182–192.
  * Hughes et al., (2017) Hughes, T. P., Kerry, J. T., Álvarez-Noriega, M., Álvarez-Romero, J. G., Anderson, K. D., Baird, A. H., Babcock, R. C., Beger, M., Bellwood, D. R., Berkelmans, R., Bridge, T. C., Butler, I. R., Byrne, M., Cantin, N. E., Comeau, S., Connolly, S. R., Cumming, G. S., Dalton, S. J., Diaz-Pulido, G., Eakin, C. M., Figueira, W. F., Gilmour, J. P., Harrison, H. B., Heron, S. F., Hoey, A. S., Hobbs, J.-P. A., Hoogenboom, M. O., Kennedy, E. V., Kuo, C.-y., Lough, J. M., Lowe, R. J., Liu, G., McCulloch, M. T., Malcolm, H. A., McWilliam, M. J., Pandolfi, J. M., Pears, R. J., Pratchett, M. S., Schoepf, V., Simpson, T., Skirving, W. J., Sommer, B., Torda, G., Wachenfeld, D. R., Willis, B. L., and Wilson, S. K. (2017).  Global warming and recurrent mass bleaching of corals.  Nature, 543(7645):373–377.
  * Jin et al., (2017) Jin, W., Shah, E. T., Penington, C. J., McCue, S. W., Maini, P. K., and Simpson, M. J. (2017).  Logistic Proliferation of Cells in Scratch Assays is Delayed.  Bulletin of Mathematical Biology, 79(5):1028–1050.
  * Johns et al., (2014) Johns, K. A., Osborne, K. O., and Logan, M. (2014).  Contrasting rates of coral recovery and reassembly in coral communities on the Great Barrier Reef.  Coral Reefs, 33(3):553–563.
  * Jonker et al., (2008) Jonker, M., Johns, K., and Osbourne, K. (2008).  Surveys of benthic reef communities using underwater digital photography and counts of juvenile corals.  Standard Operating Procedure. Australian Institute of Marine Science.
  * Lapointe et al., (2019) Lapointe, B. E., Brewton, R. A., Herren, L. W., Porter, J. W., and Hu, C. (2019).  Nitrogen enrichment, altered stoichiometry, and coral reef decline at Looe Key, Florida Keys, USA: a 3-decade study.  Marine Biology, 166(8):108.
  * Logan et al., (2020) Logan, M., Hu, Z., Brinkman, R., Sun, S., Sun, X., and Schaffelke, B. (2020).  Ecosystem health report cards: An overview of frameworks and analytical methodologies.  Ecological Indicators, 113:105834.
  * MacNeil et al., (2019) MacNeil, M. A., Mellin, C., Matthews, S., Wolff, N. H., McClanahan, T. R., Devlin, M., Drovandi, C., Mengersen, K., and Graham, N. A. J. (2019).  Water quality mediates resilience on the Great Barrier Reef.  Nature Ecology & Evolution, 3(4):620–627.
  * McDonald et al., (2023) McDonald, R. A., Neuhausler, R., Robinson, M., Larsen, L. G., Harrington, H. A., and Bruna, M. (2023).  Zigzag persistence for coral reef resilience using a stochastic spatial model.  Journal of The Royal Society Interface, 20(205):20230280.
  * McLeod et al., (2022) McLeod, I. M., Hein, M. Y., Babcock, R., Bay, L., Bourne, D. G., Cook, N., Doropoulos, C., Gibbs, M., Harrison, P., Lockie, S., van Oppen, M. J. H., Mattocks, N., Page, C. A., Randall, C. J., Smith, A., Smith, H. A., Suggett, D. J., Taylor, B., Vella, K. J., Wachenfeld, D., and Boström-Einarsson, L. (2022).  Coral restoration and adaptation in australia: The first five years.  PLOS ONE, 17(11):e0273325.
  * Mumby, (2009) Mumby, P. J. (2009).  Phase shifts and the stability of macroalgal communities on Caribbean coral reefs.  Coral Reefs, 28(3):761–773.
  * Mumby et al., (2021) Mumby, P. J., Mason, R. A. B., and Hock, K. (2021).  Reconnecting reef recovery in a world of coral bleaching.  Limnology and Oceanography: Methods, 19(10):702–713.
  * Munsky et al., (2009) Munsky, B., Trinh, B., and Khammash, M. (2009).  Listening to the noise: random fluctuations reveal gene network parameters.  Molecular Systems Biology, 5(1):318.
  * Murphy et al., (2022) Murphy, R. J., Maclaren, O. J., Calabrese, A. R., Thomas, P. B., Warne, D. J., Williams, E. D., and Simpson, M. J. (2022).  Computationally efficient framework for diagnosing, understanding and predicting biphasic population growth.  Journal of The Royal Society Interface, 19(197):20220560.
  * Ortiz et al., (2021) Ortiz, J. C., Pears, R. J., Beeden, R., Dryden, J., Wolff, N. H., del C. Gomez Cabrera, M., and Mumby, P. J. (2021).  Important ecosystem function, low redundancy and high vulnerability: The trifecta argument for protecting the great barrier reef's tabular _acropora_.  Conservation Letters, 14(5):e12817.
  * Ortiz et al., (2018) Ortiz, J.-C., Wolff, N. H., Anthony, K. R. N., Devlin, M., Lewis, S., and Mumby, P. J. (2018).  Impaired recovery of the Great Barrier Reef under cumulative stress.  Science Advances, 4(7):eaar6127.
  * Phaiboun et al., (2015) Phaiboun, A., Zhang, Y., Park, B., and Kim, M. (2015).  Survival kinetics of starving bacteria is biphasic and density-dependent.  PLOS Computational Biology, 11(4):e1004198.
  * Plummer et al., (2006) Plummer, M., Best, N., Cowles, K., and Vines, K. (2006).  CODA: Convergence Diagnosis and Output Analysis for MCMC.  R News, 6(1):7–11.
  * Pratchett et al., (2021) Pratchett, M. S., Caballes, C. F., Cvitanovic, C., Raymundo, M. L., Babcock, R. C., Bonin, M. C., Bozec, Y.-M., Burn, D., Byrne, M., Castro-Sanguino, C., Chen, C. C. M., Condie, S. A., Cowan, Z.-L., Deaker, D. J., Desbiens, A., Devantier, L. M., Doherty, P. J., Doll, P. C., Doyle, J. R., Dworjanyn, S. A., Fabricius, K. E., Haywood, M. D. E., Hock, K., Hoggett, A. K., Høj, L., Keesing, J. K., Kenchington, R. A., Lang, B. J., Ling, S. D., Matthews, S. A., McCallum, H. I., Mellin, C., Mos, B., Motti, C. A., Mumby, P. J., Stump, R. J. W., Uthicke, S., Vail, L., Wolfe, K., and Wilson, S. K. (2021).  Knowledge gaps in the biology, ecology, and management of the Pacific Crown-of-Thorns sea star Acanthaster sp. on Australia’s Great Barrier Reef.  The Biological Bulletin, 241(3):330–346.
  * Richards, (1959) Richards, F. J. (1959).  A flexible growth function for empirical use.  Journal of Experimental Botany, 10(2):290–301.
  * Scheidegger, (2021) Scheidegger, A. (2021).  adaptMCMC: Implementation of a generic adaptive Markov chain Monte Carlo sampler.  https://CRAN.R-project.org/package=adaptMCMC.
  * Seymour and Bradbury, (1999) Seymour, R. M. and Bradbury, R. H. (1999).  Lengthening reef recovery times from crown-of-thorns outbreaks signal systemic degradation of the Great Barrier Reef.  Marine Ecology Progress Series, 176:1–10.
  * Simpson et al., (2022) Simpson, M. J., Browning, A. P., Warne, D. J., Maclaren, O. J., and Baker, R. E. (2022).  Parameter identifiability and model selection for sigmoid population growth models.  Journal of Theoretical Biology, 535:110998.
  * Simpson et al., (2023) Simpson, M. J., Walker, S. A., Studerus, E. N., McCue, S. W., Murphy, R. J., and Maclaren, O. J. (2023).  Profile likelihood-based parameter and predictive interval analysis guides model choice for ecological population dynamics.  Mathematical Biosciences, 355:108950.
  * Sisson et al., (2018) Sisson, S. A., Fan, Y., and Beaumont, M., editors (2018).  Handbook of Approximate Bayesian Computation.  Chapman and Hall/CRC.
  * Soetaert et al., (2010) Soetaert, K., Petzoldt, T., and Setzer, R. W. (2010).  Solving differential equations in r: Package desolve.  Journal of Statistical Software, 33(9):1–25.
  * Tebbett et al., (2023) Tebbett, S. B., Connolly, S. R., and Bellwood, D. R. (2023).  Benthic composition changes on coral reefs at global scales.  Nature Ecology & Evolution, 7(1):71–81.
  * Thompson et al., (2016) Thompson, A. A., Costello, P., Davidson, J., Logan, M., Coleman, G., Gunn, K., and Schaffelke, B. (2016).  Annual report for inshore coral reef monitoring.  Technical report, Australian Institute for Marine Science.
  * Thompson and Dolman, (2010) Thompson, A. A. and Dolman, A. M. (2010).  Coral bleaching: one disturbance too many for near-shore reefs of the Great Barrier Reef.  Coral Reefs, 29(3):637–648.
  * Thompson et al., (2020) Thompson, A. A., Martin, K., and Logan, M. (2020).  Development of the coral index, a summary of coral reef resilience as a guide for management.  Journal of Environmental Management, 271:111038.
  * Tjørve and Tjørve, (2010) Tjørve, E. and Tjørve, K. M. (2010).  A unified approach to the Richards-model family for use in growth analyses: Why we need only two model forms.  Journal of Theoretical Biology, 267(3):417–425.
  * van Woesik et al., (2018) van Woesik, R., Köksal, S., Ünal, A., Cacciapaglia, C. W., and Randall, C. J. (2018).  Predicting coral dynamics through climate change.  Scientific Reports, 8(1):17997.
  * Vehtari et al., (2021) Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2021).  Rank-Normalization, Folding, and Localization: An Improved $\widehat{R}$ for Assessing Convergence of MCMC.  Bayesian Analysis, 16(2):667–718.
  * Vercelloni et al., (2017) Vercelloni, J., Caley, M. J., and Mengersen, K. (2017).  Crown-of-thorns starfish undermine the resilience of coral populations on the Great Barrier Reef.  Global Ecology and Biogeography, 26(7):846–853.
  * Vercelloni et al., (2020) Vercelloni, J., Liquet, B., Kennedy, E. V., González-Rivero, M., Caley, M. J., Peterson, E. E., Puotinen, M., Hoegh-Guldberg, O., and Mengersen, K. (2020).  Forecasting intensifying disturbance effects on coral reefs.  Global Change Biology, 26(5):2785–2797.
  * Vihola, (2012) Vihola, M. (2012).  Robust adaptive Metropolis algorithm with coerced acceptance rate.  Statistics and Computing, 22(5):997–1008.
  * Wang et al., (2024) Wang, X., Jenner, A. L., Salomone, R., Warne, D. J., and Drovandi, C. (2024).  Calibration of agent based models for monophasic and biphasic tumour growth using approximate Bayesian computation.  Journal of Mathematical Biology, 88(3):28.
  * Warne et al., (2019) Warne, D. J., Baker, R. E., and Simpson, M. J. (2019).  Using experimental data and information criteria to guide model selection for reaction–diffusion problems in mathematical biology.  Bulletin of Mathematical Biology, 81(6):1760–1804.
  * Warne et al., (2020) Warne, D. J., Baker, R. E., and Simpson, M. J. (2020).  A practical guide to pseudo-marginal methods for computational inference in systems biology.  Journal of Theoretical Biology, 496:110255.
  * Warne et al., (2022) Warne, D. J., Crossman, K. A., Jin, W., Mengersen, K., Osborne, K., Simpson, M. J., Thompson, A. A., Wu, P., and Ortiz, J.-C. (2022).  Identification of two-phase recovery for interpretation of coral reef monitoring data.  Journal of Applied Ecology, 59(1):153–164.
  * Warne et al., (2024) Warne, D. J., Maclaren, O. J., Carr, E. J., Simpson, M. J., and Drovandi, C. (2024).  Generalised likelihood profiles for models with intractable likelihoods.  Statistics and Computing, 34(1):50.
  * Álvarez–Noriega et al., (2023) Álvarez–Noriega, M., Madin, J. S., Baird, A. H., Dornelas, M., and Connolly, S. R. (2023).  Disturbance-induced changes in population size structure promote coral biodiversity.  The American Naturalist, 202(5):604–615.


Discrete-time stochastic optimal control remains a challenging problem for general, nonlinear systems under significant uncertainty, with practical solvers typically relying on the certainty equivalence assumption, replanning and/or extensive regularization.
Control as inference is an approach that frames stochastic control as an equivalent inference problem, and has demonstrated desirable qualities over existing methods, namely in exploration and regularization.
We look specifically at the input inference for control () algorithm,
and derive three key characteristics that enable advanced trajectory optimization: An `expert' linear Gaussian controller that combines the benefits of open-loop optima and closed-loop variance reduction when optimizing for nonlinear systems, inherent adaptive risk sensitivity from the inference formulation, and covariance control functionality with only a minor algorithmic adjustment.

§ INTRODUCTION

The control of stochastic environments is a ubiquitous problem across many domains, but remains challenging computationally for the general case.
Stochastic optimal control (SOC) solvers trade off computational complexity, exploitation of domain knowledge, use of simplifying assumptions and/or numerical sensitivity.
In this work, we discuss these trade-offs from the perspective of control as inference <cit.>, which seeks to frame stochastic control as a probabilistic inference problem.
This translation from optimization to inference can be seen as a subset of probabilistic numerics <cit.>, which utilize statistical methods to solve numerical problems, providing uncertainty quantification, regularization and faster convergence.
These viewpoints are made considering input inference for control () <cit.>, a fully probabilistic inference-based solver that frames SOC as input estimation.
Unlike classical control theory, an inference-based approach naturally lends itself to the manipulation of uncertainties, while also benefiting from mature approximate inference methods for complex dynamics and uncertainties where exact inference is intractable.
Considering nonlinear trajectory optimization, without access to closed-form solutions, a key quality is the ability to iteratively explore in a stable manner.
Many SOC algorithms require regularization heuristics such as line search or trust regions to achieve this, as well as initializing with a sufficiently random solution to encourage progress.
We show inference methods naturally achieve this, using belief akin to a trust region and leveraging an adaptive risk-seeking strategy for exploration.
Another numerical issue in nonlinear SOC is optimizing open- or closed-loop.
While open-loop strategies are brittle, yet simpler to compute, closed-loop controllers offer additional stability and therefore reduce the variance of the state distribution.
However optimizing with feedback tends to yield control-heavy solutions due to the interplay between exploration and local feedback during optimization.
We show how using the belief in the controller during optimization allows for this interplay to be managed, producing superior results on two nonlinear tasks.
Covariance control is also implemented using a minor adjustment to , due to its similarity to inference, enabling nonlinear distributional control that is simpler compared to alternative approaches.

The paper is structured as follows. Section II discuses the SOC literature and related theory, Section III described control as inference and the  algorithm, and Section III details the extensions to  than enable advanced trajectory optimization.

§ STOCHASTIC OPTIMAL CONTROL

In this section we outline SOC, baseline methods, the risk-based variant and covariance control.

§.§ The Finite-Horizon, Discrete-Time Setting

We specifically consider a stochastic discrete-time fully-observed nonlinear time-varying dynamical system, $\vf_t$, with state $\vx \in \mathbb{R}^{d_x}$ and input $\vu \in \mathbb{R}^{d_u}$, desiring the optimal controls over a time horizon $T$ that minimizes the cost functions $C_{0:T}$ in expectation

\begin{align}
&\min_{\vu_{0:T\text{-}1}}& &\mathbb{E}[C_T(\vx_T) + \sum\textstyle_{t=0}^{T\text{-}1}\;C_t(\vx_t, \vu_t)]\hspace{1cm}\\
&\text{s.t.}& &\vx_{t+1} = \vf_t(\vx_t, \vu_t).
\end{align}

From Bellman's dynamic programming perspective, the optimal controls
can be derived from the value function $V_t$ for the expected cost.
For $t = 0,\dots,T\text{-}1$,

\begin{align}
V_t(\vx) &= \min_{\vu}\mathbb{E}[C_t(\vx, \vu) + V_{t+1}(\vf_t(\vx, \vu))]\label{eq:value},
\end{align}

defining $V_T$ using the terminal cost.
For time-varying linear dynamics, quadratic costs and Gaussian disturbances, the classic LQG solution can be derived in closed-form from Equation <ref>.
However, in practice this value function is intractable globally given arbitrary dynamics, objectives and uncertainties.
As a result, methods involve either devising a means of approximating the value function, or tackling the constrained optimization problem directly instead, ignoring the temporal structure.

Differential dynamic programming (DDP) <cit.> exploits the temporal structure of the SOC task, iteratively constructing local Taylor series approximations of the dynamics and cost to compute a local approximate value function for closed-loop optimization.
The stochastic DDP extension (SDDP) to consider the impact of Gaussian disturbances on the expected value functions (SDDP) <cit.>.
To mitigate the computational burden of computing the Hessian, a Gauss-Newton approximation of DDP is used (iLQR <cit.>, iLQG <cit.>).
These methods all require regularization and line-search routines for stable convergence.
Guided policy search (GPS) is a method that uses iLQG to solve a variational form of trajectory optimization <cit.>.
This results in a maximum entropy linear Gaussian control law, and a KL divergence constraint on the trajectory distribution is used to stabilize the optimization.

For nonlinear, deterministic systems, sequential quadratic programming (SQP) has been used to compute the optimal action- or state-action sequence when framed as a constrained optimization problem, with quadratic convergence under mild assumptions <cit.>.
Linearizing along this trajectory, the LQG solution has then been applied to compute an approximate solution to the SOC problem (T-LQG <cit.>, T-PFC <cit.>), justified by a small noise performance bound <cit.>.
While mature and highly optimized SQP solvers can be used off the shelf, they can suffer computationally in the SOC domain.
By failing to exploit the temporal structure of the problem (multiple-shooting) or by considering numerically sensitive objectives (single-shooting), SQPs do not scale gracefully under long planning horizons.

§.§ Covariance Control

While SOC solvers are typically concerned with optimizing the expected cost over the mean trajectory, methods have also been devised to control higher order statistics.
Covariance control <cit.> specifically looks at constraining the mean and covariance of the terminal state to a target distribution $\vx_T^*$.
The linear Gaussian setting as been extensively studied, for both discrete <cit.> and continuous time <cit.>, where it can be shown that a solution exists should the system be controllable and $\sigm{\vx}{*}{T}{-}\sigEta{T}{>}0$, given process noise covariance $\mSigma_{\veta_t}$.
The hard constraint can be tackled by decomposing the problem into feedforward control for the mean, and linear feedback control for the covariance <cit.>.
The discrete-time case corresponds to minimizing the relative entropy between two MDPs and minimum-energy LQG <cit.>, where the terminal cost corresponds to the Lagrange multiplier of the constraint.
The nonlinear Gaussian case has been tackled using stochastic DDP <cit.> and through the combination of sequential convex programs and statistical linearization <cit.>.
The problem can also be viewed as a form of optimal transport and the Schrödinger bridge, which seeks to find the mapping (i.e. dynamical system) that transforms one distribution to another <cit.>.
The solution here is iterative, using forward and backward Riccati equations until both constraints are satisfied.

§.§ Risk-Sensitive Control

Introduced by Jacobson, risk-sensitive linear exponential quadratic Gaussian (LEQG) control <cit.> derives a control law that, unlike LQG, is dependent on the severity of uncertainty in the system dynamics. This sensitivity is determined by the scaling parameter $\sigma$ in the now exponentiated objective,

\begin{align}
\min_{\vu_{0:T\text{-}1}}
\text{-}
\frac{2}{\sigma}
\left[
\exp
\left(
\text{-}
\frac{\sigma}{2}
\left[
\right)
\right],
\label{eq:risk}
\end{align}

for quadratic costs $C_{0:T}$ and $\sigma\in\mathbb{R}$.
The consequence of this objective transformation is that, while they still form a LQG-like discrete algebraic Riccati equation (DARE), the weights of the value function now also depend on the covariance of the Gaussian disturbance $\veta_t$ and risk sensitivity $\sigma$, resulting in the transformation,

\begin{align}
\mP_{\text{LEQG},t} &= (\mP_{\text{LQG},t}\inv+\sigma\sigEta{t})\inv
\text{ given } V_t(\vx_t) &=
\vx_t\tran\mP_t\vx_t + 2\vp_t\tran\vx_t + p_t.
\end{align}

While the relationship between the linear feedback law and the value function is unchanged from the LQG case, the adjustment of the value function results in `risk-preferring' ($\sigma{>}0$) or `risk-averse' ($\sigma{<}0$) strategies under uncertainty.
On a high level, this can be viewed as a mean / variance trade-off in the evaluated cost.
Moreover, for $\sigma{=}0$ (risk-neutral) the control strategy reduces to the LQG result.
This behavior has lead to interest in risk-sensitive methods from domains such as quantitative finance and behavioral sciences.
The relationship between dynamics uncertainty and risk is what makes this formulation interesting from the inference perspective, which is discussed in Section <ref>.

§ CONTROL AS INFERENCE

\begin{align}
p(\mX, \mU, \mZ, \vtheta)=
\textstyle\prod_{t=0}^{T{\text{-}}1}
\textstyle\prod_{t=0}^{T{\text{-}}1}
\textstyle\prod_{t=0}^{T{\text{-}}1}
\label{eq:i2cmodel}
\end{align}

\begin{align}
&\text{Dynamics} &p(\vx_{t+1}|\vx_{t},\vu_t) &:
&\vx_{t+1} &= \vf_t(\vx_{t},\vu_t) + \veta_t,
&\veta_t &\sim \mathcal{N}(\bm{0}, \sigEta{t}),
\label{eq:dynamics}\\
&\text{Cost} &p(\vz_t|\vx_t,\vu_t,\vtheta)&:
&\vz_t &= \vg_t(\vx_{t},\vu_t) + \vxi_t,
&\vxi_t &\sim \mathcal{N}(\bm{0}, \sigXi(\vtheta)),
\label{eq:observation}\end{align}

The duality between optimal control and statistical inference techniques dates back to the work of Kalman <cit.> from his work on the LQG problem.
Later, relative entropy (KL divergence) was shown as a means of framing the SOC problem <cit.>, also known as `probabilistic control design' <cit.>.
While path integral control <cit.> strengthened the connection in continuous-time, adopting methods from probabilistic graphical models demonstrated the relation between message passing for inference and the discrete-time Riccati equations in SOC <cit.>.
AICO <cit.> demonstrated that, for open-loop nonlinear trajectory optimization, the linearization-based approximate message passing computation resembled Gauss-Newton SOC, but converged faster due to its exploratory forward pass and use of priors over control inputs.
However, AICO's exploration prior was fixed due to its dual role as the control regularization, and the translation of the cost to the distributions was hand tuned.
Input inference for control <cit.> built on AICO, extending the probabilistic perspective by framing optimal control as input estimation.
The closed-loop controller was defined as the conditional distribution of the state-action trajectory, as in posterior policy iteration (PPI) <cit.>, and crucially the graphical model was defined such that the controls had independent priors, allowing for iterative, variable exploration that corresponds to maximum entropy control.
Moreover, the translation of the cost function for the graphical model, through a scaling term, was jointly optimized during inference.
The combination of principled exploration and probabilistic regularization enable competitive performance against comparable algorithms for nonlinear optimal control.
To understand how  works, there are four critical components: The latent variable model, expectation maximization, message passing and the transformation of the control cost functions into likelihoods.

§.§ Latent Variable Models

Like in state estimation,  uses a sequential latent variable model for the state-action trajectory over time, however in the control as inference setting there is no data.
Instead, the `known' quantity is the desired trajectory $\mZ{=}\vz_{0:T}$, which is some transformation $\vg_t(.)$ of the state-action space.
The transformation is tasks specific, e.g. $\vz{=}[\vx, \vu]\tran$ for LQR, $\vz{=}\vu$ for minimum-energy or $\vz$ is in cartesian coordinates for operation space control for a robot manipulator.
Regardless of form, it is assumed to be fixed and part of the objective.
The inference problem is therefore computing the most likely state-action distribution that reconstructs this desired trajectory, which can be done through optimizing the joint likelihood (eq:i2cmodel).
As $\mZ$ is designed, rather than sampled data, the belief in  is with respect to optimality, rather than typical statistical uncertainty.

§.§ Cost Functions and Constraints as Likelihoods

In practice, this desired trajectory $\mZ$ will not be perfectly reconstructed.
Therefore, we augment the observation model with additive uncertainty to account for this mismatch (eq:observation).
In Section <ref>, we discussed LQR for its tractable quadratic the control costs.
These distances can be expressed more concisely as a Mahalanobis distance,
Assuming Gaussian distributions, we can draw a connection between LQR and  as the `energy' / log-likelihood of this distribution is also a Mahalanobis distance.
For Gaussian state space models, the log-likelihood takes an attractive quadratic form,

\begin{align}
\vy &= \vf(\vx) + \bm{\eta}, \hspace{0.7cm}\bm{\eta} \sim \mathcal{N}(\bm{0},\mSigma_{\eta}), \\
-\mathcal{L}_p(\vy, \vx)
&= -\log p(\vy, \vx) \notag\\
\lVert\vy-\vf(\vx)\rVert_{\mSigma_{\eta}}^2
+ d_x\log|2\pi\mSigma_{\eta}|.\label{eq:loglike}
\end{align}

However, while the control objective in LQG is affine-invariant, the loglikehood is not.
Therefore, for equivalence with the observation loglikehood, there is an unknown scale factor ($\alpha$) to relate the LQG cost weights ($\mQ,\mR$) to the precision (inverse covariance) of the disturbance $\vxi$ in ,
In Section <ref>, we discuss how it relates to risk sensitivity.
As the only unknown model parameter in the inference problem ($\vtheta{=}\{\alpha\}$),
a benefit of  is that $\alpha$ can be automatically tuned during inference.

§.§ Inference of the Graphical Model

From eq:i2cmodel, we see that the general  inference problem depends on two unknowns: the latent state-action distribution $\mX,\mU$ and model parameters $\vtheta$.
Performing inference with both these unknown quantities jointly is intractable.
Fortunately, for a Gaussian dynamical system computation can be achieved iteratively using expectation maximization (EM) <cit.>, which guarantees monotonic improvement.
In the EM convention, the E-Step corresponds to estimating the latent state-action distribution given the model, while the M-Step the optimizes the model parameters $\vtheta$ to maximize the expected log-likelihood given the estimated latent distribution.
Using eq:loglike, this is of the form

\begin{align}
-\mathbb{E}[\mathcal{L}_p(\vy, \vx)] =
&\text{tr}\{\sigEta{~}\inv\mathbb{E}[(\vy{-}\vf(\vx))(\vy{-}\vf(\vx))\tran]\} \notag\\
&+ d_x\log|2\pi\bm{\Sigma}_{\veta}|.\label{eq:expectedloglike}
\end{align}

This EM procedure is iterated till convergence.
While the M-Step improvement can be expressed in closed-form <cit.>, the E-Step requires closer consideration.

The factor graph of  for a single timestep.

To perform the E-Step efficiently, message passing methods are a flexible approach to designing such algorithms <cit.>.
Passing linear Gaussian messages is especially straightforward, following simple rules that can be expressed in closed-form.
For nonlinear Gaussian messages, approximate inference can be performed in several ways.
Akin to the various nonlinear Bayesian filters, Taylor series approximations, quadrature and Monte Carlo methods can be used to approximate the marginalization integral <cit.>.
As the graph is a Gaussian dynamical system, the E-Step relates closely to general recursive Bayesian smoothing, using Equations <ref> and <ref> as the Gaussian dynamics and `observation' model respectively.
However, the key separation from typical smoothing is the inclusion of the input as part of the probabilistic state.
With a prior on the state-action distribution at time $t$,
during the forward message passing we compute
$p(\overrightarrow{\vx}_t, \overrightarrow{\vu}_t)$
from the incoming $\overrightarrow{\vx}_t$ computed in the previous timestep by marginalizing over the conditional distribution,

\begin{align}
p(\overrightarrow{\vu}_t) &=
\int p(\vu_t|\vx_t{=}\overrightarrow{\vx}_t)p(\overrightarrow{\vx}_t)d\overrightarrow{\vx}_t,\\
\mur{\vu}{}{t} &= \mum{\vu}{}{t} + \mK_t(\mur{\vx}{}{t}-\mum{\vx}{}{t}),\\
\sigr{\vu\vu}{}{t} &= \sigm{\vu}{}{t} - \sigm{\vu\vx}{}{t}\sigm{\vx\vx}{}{t}\inv\sigm{\vu\vx}{}{t}\tran
+ \mK_t \sigr{\vx}{}{t} \mK_t\tran,\\
\sigr{\vu\vx}{}{t} &= \sigm{\vu\vx}{}{t},\;\;
\mK_t = \sigm{\vu\vx}{}{t}\sigm{\vx\vx}{}{t}\inv.
\end{align}

Where $\mK_t$ is the  feedback control gain, which has been shown to be equivalent to maximum entropy LQR in the deterministic setting <cit.>.
Note that if the state and action are independent (i.e. $\sigm{ux}{}{t}{=}\bm{0}$) or the state distribution is unchanged ($\overrightarrow{\vx}_t{=}\vx_t$), then no `feedback' is applied and the input distribution remains as the prior $\overrightarrow{\vu}_t{=}\vu_t$.
This feedback term also exhibits the `turn-off phenomena' from dual control <cit.>, as the gain is attenutated as the state uncertainty increases.

§ INFERENCE-BASED ADVANTAGES

In this section we introduce algorithmic improvements for  and explain existing qualities from the SOC perspective.

§.§ Optimizing Expert Linear Gaussian Controllers

When performing trajectory optimization, considering the open- or closed-loop setting is a crucial distinction.
The open-loop approach (i.e. using SQPs) is a simpler optimization problem, but may be brittle and is limited to shorter planning horizons and levels of stochasticity.
On the other hand, closed-loop optimization (i.e. DDP) provides extra stability with the local control law, however the numerical procedure has unfortunate side-effect in the nonlinear setting: the feedback fights exploration during optimization, resulting in highly-actuated, and therefore likely sub-optimal, solutions.
As  derives its controller from the joint state-action distribution, both open- and closed-loop optimization can be incorporated into the E-step by considering the independent ($\sigr{\vu\vx}{}{t}{=}\,\bm{0}$) or full joint distribution respectively.

In Table 1 we compare these variations on two stochasic, nonlinear swing-up tasks, looking at open- and closed-loop optimization for both  and baseline SOC solvers, which also vary between considering the actual stochastic problem or a certainty equivalent approximation.
For both tasks, open-loop methods resulted in better optima but were also high variance in the cost, while the closed-loop alternatives had much lower variance but sub-optimal performance on the simulated systems due to their over-actuation.
Reassuringly, the results of the  variant and equivalent baseline solver were generally similar due to the comparable computation.

These results then suggest that ideally we should be able to combine the benefits of open- and closed-loop optimization to improve performance during inference.
With  this is indeed possible, by adjusting the conditional distribution ($\pi_t$) to act as an `expert' controller, which utilizes the probability that the local controller applies to the current state $p(\vx_t{=}\vx)$,

\begin{align}
\pi_t(\vu|\vx) &= p(\vx_t{=}\vx)\,p(\vu_t\mid\vx_t{=}\vx) + p(\vx_t{\neq}\vx)\,p(\vu_t).
\end{align}

This weighting softens the controller to fall back to the open-loop controls $p(\vu_t)$ as $\vx$ moves away from estimated optimal state trajectory.
Therefore, the feedback effect is reduced during significant exploration in the E-step, and as the open-loop controls are well regularized this avoids highly actuated trajectories forming.
Table 1 demonstrates the effectiveness of this addition, where this expert controller matches the open-loop optima but with the closed-loop variance reduction.
For a Gaussian distribution, $p(\vx_t{=}\vx)$ is defined as the confidence interval $\lVert\vx-\mum{\vx}{}{t}\rVert_{\mSigma_{\vx_t}}^2\leq\mathcal{X}^2_k(p)$ for probability $p$, where $\mathcal{X}^2_k$ is the chi squared distribution.
Another interesting aspect of the expert controller results is how certainty equivalence outperforms optimizing the true stochastic setting.
In practice, we also found the stochastic setting required greater care regarding hyperparameter tuning.
This is because under greater uncertainty the inaccuracy in the approximate inference is more severe, thus requiring greater regularization during inference.
Moreover, from the probabilistic numerics view, greater uncertainty corresponds to greater numerical regularization e.g. the turn-off phenomena in the feedback gains, which may also be at play.
Studying this phenomena is greater detail, especially with more accurate approximate inference techniques, will be a topic of future work.

§.§ Adaptive Risk-Sensitivity

For Gaussian , the expected likelihood in eq:i2cmodel may be reformulated to expose the quadratic cost function from the observation likelihood,

\begin{equation}
\text{-}
{\frac{\alpha}{2}} %
\left[
\right]
\right)}
\right]
\end{equation}

where $A(.)$ contains the remaining likelihood and normalizing terms.
Comparing this form to eq:risk, $\alpha$ may be interpreted as equivalent to risk sensitivity.
For both LEQG and , $\sigma{=}\alpha{=}0$ corresponds to the deterministic LQR result.
However, for , $\alpha>0$ due to the probabilistic treatment, resulting in only risk-seeking behavior possible, with risk neural control as the limit of $\alpha{\rightarrow}0$.
Risk-seeking behavior is an inherent issue for control as inference methods, due to the probabilistic formulation naturally attenuating the effect of unlikely trajectories in the objective.
Related methods such as GPS and path integral control also share this property.
However, due to the other terms at play in , $\alpha$ and $\sigma$ do not have a direct correspondence, nor do their resultant trajectories exactly match.
This discrepancy can be identified by examining the linear Gaussian message passing equations.
First note in the deterministic setting (without $\bm{\eta}_t$),

\begin{align}
\lamm{\vx}{''}{t} &= \lamm{\vx}{}{t+1}, \text{where
$\vx_t{''}$ denotes the state after $\vf$.}
\intertext{Moreover in the LQR equivalent case, where there are sufficiently uninformative priors placed on $\mU_t$,}
\lamm{\vx}{}{t+1} &= \lamr{\vx}{}{t+1}{+}\laml{\vx}{}{t+1} \approx \laml{\vx}{}{t+1}\;\text{for large $\sigr{\vx}{}{t+1}$}.
\end{align}

This provides an inference-based perspective why LQR can be solved without considering the forward propagation.
However, considering disturbance $\bm{\eta}_t$,

\begin{align}
\sigm{\vx}{}{t+1} &= \sigm{\vx}{''}{t} + \sigEta{t},\;\;
\lamm{\vx}{}{t+1} = (\lamm{\vx}{''}{t}\inv + \sigEta{t})\inv,\\
\lamm{\vx}{''}{t}\inv &= (\lamr{\vx}{''}{t}{+}\laml{\vx}{''}{t})\inv,\\
\laml{\vx}{''}{t}\inv{-}\laml{\vx}{''}{t}\inv (\lamr{\vx}{''}{t}\inv{+}\laml{\vx}{''}{t}\inv)\inv\laml{\vx}{''}{t}\inv,\\
\lamm{\vx}{}{t+1} &\approx (\laml{\vx}{''}{t}\inv + \hat{\bm{\Sigma}}_{\veta_t})\inv.
\end{align}

Therefore, LEQG's risk sensitivity may be viewed as a means to approximate the prior state distribution while using LQR's Riccati equation, through the transformed disturbance term $\hat{\bm{\Sigma}}_{\veta_t}{=}\sigma\sigEta{t}$.
Moreover, we can see that setting $\sigma>0$ (i.e. risk avoidance) acts to reduce the uncertainty in the trajectory in a way that probabilistically invalid, as $\hat{\bm{\Sigma}}_{\veta_t}$ must be positive definite.

While risk-seeking behavior is typically seldom desired, $\alpha$'s influence on both exploration and risk suggests that risk-seeking control is useful (for nonlinear systems) by greedily seeking the desired trajectory.
Moreover, further consideration of this term may improve the robustness of the controller by limiting risky strategies.

§.§ Covariance Control as Inference

For , the idea of a terminal cost can be tackled in two ways.
Previously, we have used the idea of a terminal observation function $\vg_T(\vx)$, which allows for arbitrary terminal costs as described in Section <ref>.
However, another approach is to work with the terminal latent state distribution.
In other time-series inference applications, e.g. state estimation, the terminal state posterior is set to the prior as there is no additional information to use.
However, for control we can avoid the cost function design and the cost-to-likelihood translation by setting the terminal distribution directly, but are now limited to stipulating the desired state directly rather than a transformation (e.g. $\cos\theta{=}1$).
This approach is equivalent to covariance control (Section <ref>) as, like solving a Schrödinger bridge, iterations of forward and backward Riccati equations are performed until the boundary condition are satisfied.
In this setting, $\mSigma_{\vxi_T}$ now acts as the Lagrange multiplier.
Examining the expected log-likelihood (eq:expectedloglike) term for the terminal state ($\mathcal{L}_{\vz_T}$) for the direct state optimization case $\vz_T{=}\vx_T{+}\vxi_T$, we get

\begin{align}
\mathcal{L}_{\vz_T}
&\propto\lVert\mum{\vz}{}{T}{-}\mum{\vx}{}{T}\rVert_{\sigXi}^2{+} \text{tr}\{\mSigma_{\vxi_T}\inv\sigm{\vx}{~}{T}\}+\text{const.}\end{align}

which is equivalent the LQG correspondence proved in Goldshtein et al. (Equation 41) <cit.>, where $\mSigma_{\vxi_T}\inv$ is the terminal cost / Lagrange multiplier matrix.
However, rather than compute this term, using our probabilistic framework we can set the terminal distribution directly via the posterior of $\vx_T$.
In this case, the inference iteratively seeks to satisfy the boundary conditions on $\vx_0$ and $\vx_T$.
Figure <ref> demonstrates covariance control on a linear Gaussian system, where inference is exact.
It should be noted that  uses linear Gaussian controllers, whose uncertainty is required to achieve the desired state distribution, whereas the previous literature solves the task using deterministic control laws.
This variation of  naturally translates to nonlinear systems (Figure <ref>), avoiding the complexity of the additional forward sampling required for DDP-based covariance control <cit.>.
However, the terminal boundary constraint still requires a means of being applied in an gradual manner, due to the iterative aspect of the nonlinear optimization.
The terminal state distribution can be shifted from the prior to the constraint by `annealing' <cit.> the prior during inference

\begin{align}
p(\vx_T) = p(\vx_T^*)\,p(\overrightarrow{\vx}_T)^{\beta},\; 1 \geq \beta \geq 0,
\end{align}

where $\beta$ is the annealing temperature on the prior from the forward pass.
By decreasing $\beta$ over iterations, $p(\vx_T){\rightarrow}p(\vx_T^*)$ as $\beta{\rightarrow}0$.
This annealing strategy has previously been applied to regularize inference algorithms.

for minimum-energy exact linear Gaussian covariance control on an unstable system, with a fixed small $\alpha$ and $\sigEta{t}{=}\text{diag}(0.1, 0.1)$.
The KL divergence is between the terminal goal and closed-loop distributions.

Nonlinear minimum-energy covariance control on the pendulum swing-up task, using  with approximate inference. Plot depicts the inferred trajectory <ref> for target distribution <ref>, with simulated rollouts <ref>.

§ CONCLUSION

In this work, we have discussed control as inference and the  algorithm from the perspective of stochastic optimal control.
We have analysed  with respect to the risk-sensitive control, highlighting the similarities between LEQG and inference and how  utilizes adaptive risk during EM.
Moreover, by analyzing the numerical consequences of open- and closed-loop optimization with  and baseline solvers, we have motivated an expert linear Gaussian controller that leverages the state belief to balance exploration with stabilizing feedback, which achieved superior results on simulated nonlinear control tasks.
Finally, we demonstrated how covariance control can be implemented with a minor adjustment to , enabling exact and approximate solutions in the linear and nonlinear setting respectively.

§ ACKNOWLEDGMENTS

This project has received funding from the European Union’s Horizon 2020  research and innovation programme under grant agreement No. 713010 (GOAL Robots).
The paper was inspired in part by many insightful discussions with Hany Abdulsamad and Boris Belousov.

Experimental Details

For the experiment reported in Table 1, task parameters are similar to prior work <cit.>, but with the Cartpole environment now operating at a 4ms timestep for a time horizon of 500.
For , cubature quadrature was used for approximate inference <cit.>.
T-PFC was implemented using the transcription method (i.e. $(T{+}1)d_x{+}Td_u$ state space and $T{+}1$ constraints)