This document will describe the new approach and prove its sampling
from the same distribution, followed by a description of the
computational changes.
Equivalence of the samplers
Let \(\mathbf{\phi}\) be a vector
containing the fixed parameters, specifically the (vectorised)
covariance matrices, the (vectorised) AR(1) parameter matrices as well
as the shared long-term discrepancy \(\delta\): \[
\mathbf{\phi} =
(\Lambda_{1},\Lambda_{2},\ldots,\Lambda_{m},\Lambda_{\eta},\Lambda_{y},C_{\gamma},R_{1},
R_{2},\ldots R_{m},R_{\eta},\mathbf{\delta})\,.
\] If we have a time series of length \(T\), the aim is to fit these parameters
given the data \(\mathbf{w}^{(1)},\mathbf{w}^{(2)},\ldots,\mathbf{w}^{(T)}\).
When the option sampler = "explicit" is chosen in
fit_ensemble_model, the sampler generates samples from the
distribution
\[
p(\mathbf{\phi},\mathbf{\theta}^{(1)},\ldots,\mathbf{\theta}^{(T)}\,|\,\mathbf{w}^{(1)},\mathbf{w}^{(2)},\ldots,\mathbf{w}^{(T)})\,.
\] A sample from the posterior of \(\mathbf{\phi}\) is then produced by simple
marginalisation
\[
p(\mathbf{\phi}\,|\,\mathbf{w}^{(1)},\ldots,\mathbf{w}^{(T)}) = \int
p(\mathbf{\phi},\mathbf{\theta}^{(1)},\ldots,\mathbf{\theta}^{(T)}\,|\,\mathbf{w}^{(1)},\mathbf{w}^{(2)},\ldots,\mathbf{w}^{(T)})
\,d\mathbf{\theta}^{(1)}\cdots d\mathbf{\theta}^{(T)}\,.
\] Alternatively, directly sampling from the joint distribution
of \(\mathbf{\phi}\) and \(\mathbf{\theta}^{(1)},\ldots,\mathbf{\theta}^{(T)}\)
can be avoided by noting that
\[
\begin{split}
&p(\mathbf{\phi}\,|\,\mathbf{w}^{(1)},\ldots,\mathbf{w}^{(T)}) =
\int
p(\mathbf{\phi},\mathbf{\theta}^{(1)},\ldots,\mathbf{\theta}^{(T)}\,|\,\mathbf{w}^{(1)},\mathbf{w}^{(2)},\ldots,\mathbf{w}^{(T)})
\,d\mathbf{\theta}^{(1)}\cdots d\mathbf{\theta}^{(T)}\\
&\propto \int
p(\mathbf{w}^{(1)},\mathbf{w}^{(2)},\ldots,\mathbf{w}^{(T)}\,|\,\mathbf{\phi},\mathbf{\theta}^{(1)},\ldots,\mathbf{\theta}^{(T)})p(\mathbf{\phi},\mathbf{\theta}^{(1)},\ldots,\mathbf{\theta}^{(T)})
\,d\mathbf{\theta}^{(1)}\cdots d\mathbf{\theta}^{(T)}\\
&=p(\mathbf{\phi})\int \bigg(\prod_{t =
1}^{T}p(\mathbf{w}^{(t)}\,|\,\mathbf{\phi},\mathbf{\theta}^{(t)})\bigg)\bigg(\prod_{t
= 1}^{T - 1}p(\mathbf{\theta}^{(t +
1)}\,|\,\mathbf{\phi},\mathbf{\theta}^{(t)})\bigg)p(\mathbf{\theta}^{(1)}\,|\,\mathbf{\phi})\,d\mathbf{\theta}^{(1)}\cdots
d\mathbf{\theta}^{(T)}\\
&= p(\mathbf{\phi})\int \bigg(\prod_{t = 1}^{T-1}p(\mathbf{w}^{(t +
1)}\,|\,\mathbf{\phi},\mathbf{\theta}^{(t + 1)})p(\mathbf{\theta}^{(t +
1)}\,|\,\mathbf{\phi},\mathbf{\theta}^{(t)})\bigg)p(\mathbf{\theta}^{(1)}\,|\,\mathbf{\phi})p(\mathbf{w}^{(1)}\,|\,\mathbf{\phi},\mathbf{\theta}^{(1)})\,d\mathbf{\theta}^{(1)}\cdots
d\mathbf{\theta}^{(T)}\,.\\
\end{split}
\] The integral can be evaluated iteratively using the Kalman
filter and the fixed parameters \(\mathbf{\phi}\) can be sampled using
MCMC.
Output
The four currently available configurations of a
EcoEnsemble model were tested with the new
explicit option and compared with the kalman
option. The runs used the non-hierarchical model, specified by
lkj priors, the hierarchical model, specified by
hierarchical priors, the drivers model and the
hierarchical drivers model, both specified by
drivers = T. Code and data is taken from the vignettes, IncludingDrivers and EcoEnsemble.
Posterior Comparison
The plot below compares posterior summaries between the
explicit and kalman options. It shows the
difference in posterior means and distributions for the shared
state-space parameters used in posterior state reconstruction for 14
chains. As each chain is independent, each chain would have persistent
discrepancies between explicit and kalman if
the posterior distributions are different.
The parameters with the greatest difference is plotted below.
There are no considerable differences in the posterior output between
the explicit or kalman sampler.
Sampler Warnings
The two options induce different posterior geometries, so the stan
(Stan Development Team (2020)) No-U-Turn
Sampler’s diagnostics can differ. The table below summarises the
diagnostics for each model formulation.
|
|
Drivers
|
Model
|
Model + hierarchical
|
Drivers + hierarchical
|
|
Diagnostic
|
E
|
K
|
E
|
K
|
E
|
K
|
E
|
K
|
|
Divergences
|
1
|
57
|
4
|
13
|
0
|
9
|
8
|
37
|
|
Treedepth hits
|
0
|
0
|
3996
|
0
|
4000
|
0
|
0
|
0
|
|
Max treedepth seen
|
9
|
8
|
10
|
8
|
10
|
8
|
9
|
9
|
|
Min E-BFMI
|
0.828
|
0.825
|
0.726
|
0.777
|
0.769
|
0.797
|
0.847
|
0.82
|
|
Median E-BFMI
|
0.863
|
0.847
|
0.81
|
0.819
|
0.852
|
0.812
|
0.946
|
0.838
|
In these examples, the explicit option frequently
reaches the maximum tree depth. This is because FishSUMS, one of the
simulators included in the model, does not have estimates for sole,
causing a flat posterior geometry. Because this warning primarily
reflects sampling efficiency rather than invalid inference, it is
suppressed when it occurs in isolation. However, if multiple warnings
are given by stan, all warnings are printed. The explicit
option results in a lower ESS, which may lead to some warnings, but
these warnings may not be concerning the variable of interest in your
model. Therefore it may not have an influence on your model. See this vignette for further information.
