Component-wise iterative ensemble Kalman inversion for static Bayesian models with unknown measurement error covariance

The ensemble Kalman filter (EnKF) is a Monte Carlo approximation of the Kalman filter for high dimensional linear Gaussian state space models. EnKF methods have also been developed for parameter inference of static Bayesian models with a Gaussian likelihood, in a way that is analogous to likelihood tempering sequential Monte Carlo (SMC). These methods are commonly referred to as ensemble Kalman inversion (EKI). Unlike SMC, the inference from EKI is only asymptotically unbiased if the likelihood is linear Gaussian and the priors are Gaussian. However, EKI is significantly faster to run. Currently, a large limitation of EKI methods is that the covariance of the measurement error is assumed to be fully known. We develop a new method, which we call component-wise iterative ensemble Kalman inversion (CW-IEKI), that allows elements of the covariance matrix to be inferred alongside the model parameters at negligible extra cost. This novel method is compared to SMC on three different application examples: a model of nitrogen mineralisation in soil that is based on the Agricultural Production Systems Simulator (APSIM), a model predicting seagrass decline due to stress from water temperature and light, and a model predicting coral calcification rates. On all of these examples, we find that CW-IEKI has relatively similar predictive performance to SMC, albeit with greater uncertainty, and it has a significantly faster run time.


Introduction
Consider the following statistical model where y ∈ Y ⊆ R dy×1 are the observations, θ ∈ Θ ⊆ R d θ ×1 are the parameters, G : Θ → Y is a deterministic mathematical model, and η ∈ Y are the measurement errors.The number of observations is d y , the number of parameters is d θ and Γ ∈ R dy×dy is the covariance matrix characterising the measurement errors.Our interest is in the posterior distribution of the static parameters θ conditional on the observed data y, p(θ|y) ∝ N (y | G(θ), Γ)p(θ), where N (y | G(θ), Γ) is the likelihood function and p(θ) is the prior density of θ.
Markov chain Monte Carlo (MCMC; Robert and Casella, 1999) or sequential Monte Carlo methods (SMC; Del Moral et al., 2006) can be used for asymptotically exact parameter inference of θ and Γ.These methods generally require many evaluations of G(•), however, which limits their feasibility when G(•) is computationally expensive to evaluate.
If the elements of the covariance matrix Γ are known, methods based on the ensemble Kalman filter (EnKF; Evensen, 1994a;Burgers et al., 1998) are a fast but inexact alternative to MCMC and SMC.
EnKF is a Monte Carlo approximation of the Kalman filter for state estimation of high dimensional linear Gaussian state space models (LG-SSMs).Unlike the Kalman filter, EnKF can also be applied to non-linear and non-Gaussian state space models, but the resulting inference is only asymptotically unbiased for LG-SSMs (Le Gland et al., 2009;Roth et al., 2017).Iglesias et al. (2013) extend EnKF for inverse problems, and the resulting methods are generally known as ensemble Kalman inversion (EKI).In the Bayesian setting, an initial ensemble is simulated from the prior, and iteratively updated to capture statistical properties of the EKI approximation p(θ | y) to the posterior p(θ | y).If the EKI algorithm is iterated long enough, i.e. as the number of iterations J approaches infinity, the ensemble will collapse to a single point which is a minimiser of the loss function ||y − G(θ)|| Γ .For uncertainty quantification, i.e. to obtain samples from p(θ | y), early stopping of the algorithm is essential (Iglesias et al., 2013).If G(•) is a linear function and the prior is Gaussian, EKI provides exact samples from the posterior in a single iteration (J = 1) (Iglesias, 2014;Duffield and Singh, 2021), and the ensemble converges to the maximum a posteriori (MAP) estimate as J → ∞ (Iglesias et al., 2013;Duffield and Singh, 2021).
As inverse problems are often ill-posed, regularisation is required.For EKI, regularisation is induced by the subspace property, i.e. that the final ensemble is in the linear span of the initial ensemble (Iglesias et al., 2013).Additional regularisation can be used to improve the robustness and stability of EKI for sampling from p(θ | y), and to avoid overfitting the data.To that end, Iglesias (2014) introduces an iteratively regularised extension of EKI.In the linear case, this method targets the power posterior π j (θ) ∝ p(y|θ) αj p(θ), where 0 ≤ α j ≤ 1 is the regularisation parameter at iteration j and α j+1 > α j for j = 0, . . ., J − 1 (Iglesias et al., 2018).Other forms of regularisation can also be applied (Chada et al., 2020).Iglesias et al. (2018) propose an adaptive version of the iterative EKI method of Iglesias (2014), which is analogous to density tempering SMC (Del Moral et al., 2006).An alternative adaptation method is given by Iglesias and Yang (2021).The ensemble Kalman sampler (Garbuno-Inigo et al., 2020;Ding and Li, 2021) is a variation of EKI, which perturbs the ensemble instead of the observations as in regular EKI.Duffield and Singh (2021) extend EKI for non-Gaussian likelihoods through a Gaussian approximation of the likelihood, and Wu et al. (2022) use EKI as the forward kernel in data annealing SMC.The latter method is exact, and requires much fewer evaluations of G(•) than SMC with a Metropolis-Hastings forward kernel.Also of note is the extension given by Rammay et al. (2020) to account for model misspecification.
Currently, a limiting factor of all these methods is that the noise η is assumed to be characterised by known covariance matrix Γ.In this paper we develop a new adaptive iterative EKI method, which we call component-wise iterative ensemble Kalman inversion (CW-IEKI).This new method extends that of Iglesias et al. (2018), and it can handle the common situation where Γ contains unknown elements φ that require estimation.The method that we develop is completely analogous to density tempering SMC, which permits a direct comparison with the latter in terms of posterior accuracy and computation time.A comparison of a similar iterative EKI method (following the method of Rammay et al., 2020) and density tempering SMC is provided in Vilas et al. (2021); however, for the iterative EKI method they estimate φ as a function of the remaining parameters and they only consider one model example.In this paper we illustrate the new CW-IEKI method and perform comparisons on three model examples: a model of nitrogen mineralisation in soil that is based on the Agricultural Production Systems Simulator (APSIM) with partially known noise (Vilas et al., 2021), a model predicting seagrass decline due to cumulative stress from water temperature and light (Adams et al., 2020), and a model for predicting coral calcification rates (Galli and Solidoro, 2018).
The rest of the paper is organized as follows.Section 2 gives the background on EnKF, particle filters, likelihood tempering SMC, and iterative EKI methods.Section 3 describes our novel CW-IEKI method and Section 4 compares the performance of our method on three ecological model examples.Section 5 concludes.

Background
This section describes the filtering problem and the solutions given by the EnKF (Evensen, 1994b) and the bootstrap particle filter (Gordon et al., 1993).It also describes how likelihood tempering SMC and EKI can be used for parameter inference of the static model given in equation ( 1).We use the notation z i:j := {z i , z i+1 , . . ., z j } for j ≥ i throughout.

Ensemble Kalman Filter
The Kalman filter and the EnKF were developed to solve the filtering problem for state space models -this is often referred to as state estimation or data assimilation.Consider a state space model of the form where y 1:T are the observed data, x 0:T are the unobserved or latent states and T is the number of observations.The filtering distribution, p(x t | y 1:t ), can be solved by recursively applying the time update and the measurement update Note that p(x 0 | y 1:0 ) = µ(x 0 ) is assumed to be known.See Chapter 3 of Schön and Lindsten (2017) for more detail on the filtering problem and its solution.The Kalman filter solves (3)-( 4) analytically for LG-SSMs, are assumed to be known.Here, the time (3) and measurement (4) updates are p( respectively, where ỹt = H xt and I ∈ R dx×dx is the identity matrix.The Kalman gain at time t is where C xỹ t is the cross covariance between xt and ỹt , C ỹ ỹ t is the covariance of ỹt , and The EnKF is a Monte Carlo approximation of the Kalman filter for LG-SSMs.For non-linear, non-Gaussian SSMs, EnKF is asymptotically biased.The EnKF simulates the initial ensemble from the prior xn 0 ∼ µ(x 0 ) for n = 1, . . ., N , then for each iteration, the ensemble is updated as follows.First, the state and observation predictions are simulated , and xn t ∼ p(x t | y 1:t−1 ).Then, the ensemble is given by where C xỹ t is the sample cross covariance between xt ∈ R dx×N and ỹt ∈ R dy×N , and For LG-SSMs, xn t ∼ p(x t | y 1:t ).If the observation density is Gaussian, i.e. g(y t | x t ) = N (y t | G(x t ), Γ) with known covariance Γ, the Monte Carlo error in the calculation of the covariance matrices can be reduced (Roth et al., 2017).Let gn t = G(x n t ), then the sample cross covariance and sample covariance defined in equations ( 5) and (6) become

Particle Filters
For non-linear, non-Gaussian state space models, sequential Monte Carlo (SMC) methods give an exact solution to the filtering problem.SMC methods for dynamic models are often referred to as particle filters.As with EnKF, the bootstrap particle filter (Gordon et al., 1993) draws an initial ensemble from the prior µ(x 0 ).The particle filter then transforms the prior ensemble to samples from the filtering distribution through a sequence of reweighting, resampling and mutation steps.Given a set of weighted samples, {x n t−1 , W n t−1 } N n=1 ∼ p(x t−1 | y 1:t−1 ), the bootstrap particle filter maps these to p(x t | y 1:t ) as follows: 1. Resample the particles according to their weights, W 1:N t−1 and set W n t−1 = 1/N for n = 1, . . ., N .This gives a set of evenly weighted particles {x n t−1 , 1 N } N n=1 distributed according to p(x t−1 | y 1:t−1 ).

Simulate the state predictions using the transition density xn
), which gives a set of unweighted particles distributed according to p(x t | y 1:t−1 ).
3. Reweight the particles using the observation density w n t = g(y t | xn t ) for n = 1, . . ., N and normalise the weights to get W 1:N t .The final set of weighted particles is distributed according to p(x t | y 1:t ).
Given a set of evenly weighted samples from π j−1 (θ), likelihood tempering SMC transforms these to samples from π j (θ) as follows: 1. Reweight the particles using the ratio of the current target to the previous target, w n j = π j (θ n j−1 )/ π j−1 (θ n j−1 ) for n = 1, . . ., N and normalise the weights to get W 1:N j .This gives a set of weighted particles that are distributed according to π j (θ).
2. Resample the particles according to their weights, and set W n j−1 = 1/N for n = 1, . . ., N .3. Mutate the resampled particles using a Markov chain Monte Carlo (MCMC) kernel which targets the distribution π j (θ).
Step 2 removes the negligible weight particles and duplicates the high weight particles, and Step 3 diversifies the particles to mitigate the duplication.A common approach to mutate the particles is to use M iterations of an MCMC algorithm with π j (θ) as its invariant distribution.
The tempering parameter α j can be adapted at each iteration by setting α j such that a pre-specified effective sample size (ESS) threshold is achieved (Jasra et al., 2010).This will be some proportion of N .While the ESS cannot be calculated exactly, it can be approximated at each iteration j using the normalised weights Likelihood tempering SMC can also be used if the noise parameter Γ has unknown elements φ.In this case, the method is applied to {θ, φ} instead of θ.
If the function G(•) in equation ( 1) is expensive to compute, SMC may be prohibitively expensive to run.Each iteration j = 1, . . ., J requires a minimum of N M j evaluations of G(•), where M j is the number of MCMC repeats in iteration j.The entire algorithm requires a minimum of N J j=1 M j +1 evaluations, where the extra evaluation comes from the initial calculation of the likelihood.A less expensive, but asymptotically biased alternative to SMC for static models is ensemble Kalman inversion.

Iterate
Step 2 as desired.
Since the likelihood p(y | θ) is Gaussian, equations ( 7) and ( 8) are used for the covariance calculations.
In equation ( 7), the ensemble xt ∈ R dx×N is replaced with θj ∈ R d θ ×N .
While the prior induces regularisation through the subspace property, additional regularisation is often required to properly explore regions of high posterior support without overfitting the data (Iglesias, 2014).An iteratively regularised extension of the EKI method of Iglesias et al. (2013) is the algorithm of Iglesias (2014): 1. Sample θ n 0 ∼ p(θ) for n = 1, . . ., N .
The parameter h j for j = 1, . . ., J can be chosen adaptively using the method of Iglesias et al. (2018).At iteration j, assume that the particles must be reweighted from πj−1 (θ) to πj (θ).Analogously to likelihood tempering SMC, these weights are given by , where h j = α j − α j−1 , and the obtained w 1:N j are thereafter normalised to give W 1:N j .The parameter α j can be set so that the ESS, estimated using (9), matches some target threshold.Once α j is chosen, h j is given by α j − J−1 i=1 h i .The function G(•) is evaluated once per particle at every iteration, so that the total number of evaluations for IEKI is JN , where J is the total number of iterations and N is the number of particles or the ensemble size.This is much less than the computation required for SMC, but the IEKI assumes that Γ is known.In the next section we develop a new adaptive IEKI method that can estimate unknown parameters associated with Γ.

Component-Wise Iterative Ensemble Kalman Inversion
A strong limitation of IEKI is that Γ must be known.We extend IEKI to the case where Γ depends on some unknown parameter or parameters φ.For example, in the simplest case, this might be Γ(φ) = φI, where I ∈ R dy×dy is the identity matrix, although our method does not require this assumption to hold.The target distribution at iteration j is πj (θ, φ), which approximates the power posterior At each iteration j = 1, . . ., J, the model parameters θ and the noise parameters φ are updated component-wise conditional on the other.We refer to our method as component-wise IEKI (CW-IEKI).Our proposed procedure for CW-IEKI is as follows: 1. Sample {θ n 0 , φ n 0 } ∼ p(θ, φ) for n = 1, . . ., N .

Iterate Steps 2 and 3 until
In Step 3, the noise parameters φ are updated from the exact conditional posterior π j (φ | θ).We propose to use M iterations of a Metropolis-Hastings MCMC kernel, where the ensemble φ 1:N j−1 can be used to inform the proposal distribution for φ.If it is possible to independently sample from π j (φ | θ), then Gibbs sampling can also be used to update φ 1:N j−1 .Note that Step 3 does not require evaluation of G(•) since θ is fixed.Consequently, the total number of evaluations of G(•) for our method is the same as for standard IEKI, i.e.JN , which again is typically much less than the N J j=1 M j + 1 evaluations required for likelihood tempering SMC.See Algorithm 1 for more details.
To adapt h j , the weights are calculated in a similar way to IEKI, Note that the likelihood covariance Γ does not uniquely define the measurement error of the data in CW-IEKI as it does for standard IEKI.Since elements of Γ are estimated, it may also capture aspects arising from model misspecification.
Algorithm 1 MCMC update of the noise parameters φ.
Input: data y, ensembles θ 1:N j and φ 1:N j−1 , model evaluations g n j = G(θ n j ) for all n = 1, . . ., N and α j Output: updated ensemble of noise parameters φ

Implementation of CW-IEKI and Likelihood Tempering SMC
We compare our novel CW-IEKI method to likelihood tempering SMC on three model examples.
The first is a model of nitrogen mineralisation in soil (Vilas et al., 2021) that has relatively few parameters.The second model predicts seagrass decline due to cumulative water temperature and light stress (Adams et al., 2020), and the final model predicts coral calcification rates (Galli and Solidoro, 2018).The seagrass model has more parameters than the first model, and its marginal posteriors are roughly Gaussian.The coral model also has a relatively large number of parameters, but relatively uninformative data -the marginal posteriors of the parameters are close to the priors.
All code is implemented in MATLAB.For both CW-IEKI and SMC, the ensemble size is fixed at 1000, and the tempering schedule is adapted to achieve a target ESS of N/2 = 500 unless otherwise specified.To mutate the noise ensemble (φ 1:N j ) in CW-IEKI and the particles ({θ, φ} 1:N j ) in SMC we use a random walk Metropolis-Hastings kernel (Hastings, 1970), where the covariance of the random walk proposal is set to the covariance of the samples being mutated, i.e. cov(φ 1:N j ) for CW-IEKI and cov({θ, φ} 1:N j ) for SMC.Due to the higher number of parameters for the seagrass and coral models, the covariance of the random walk is scaled by 2.38 2 /(d θ + d φ ) for SMC, where d θ is the number of parameters in θ and d φ is the number of parameters in φ (Roberts and Rosenthal, 2001).For CW-IEKI, the number of MCMC iterations is fixed at a conservative 1000 -as no extra model evaluations are required, the cost of these iterations is relatively small.For SMC, the number of MCMC iterations is adapted at each iteration j = 1, . . ., J as follows (South et al., 2019): 1. Run S j MCMC iterations and estimate the acceptance rate p.

Adapt the total number of MCMC iterations as
3. Complete the remaining M j − S j MCMC iterations.
4. Calculate S j+1 for the next iteration as S j+1 = M j /2 .The value 1 − c is the target acceptance rate, • denotes the ceiling function and • denotes the floor function.For all models, S 1 = 5 and the target acceptance rate is 1 − 0.01 = 0.99.
We assess the performance of CW-IEKI based on its accuracy, predictive performance and computation time relative to SMC.The marginal posterior density plots of the model parameters are used to compare the accuracy of CW-IEKI to the SMC solution.As these plots do not account for parameter interdependencies however, we also compare the marginal densities of parameter combinations that greatly influence the model fit (Monsalve-Bravo et al., 2022).These combinations are identified through the eigendecomposition of a sensitivity matrix that captures key characteristics of the posterior distribution.Unless otherwise specified, we calculate the sensitivity matrix as the inverse of the sample covariance of the natural logarithm of the posterior samples from SMC (Monsalve-Bravo et al., 2022).The logarithm of the kth parameter combination is where (v k ) j is the jth element of the kth normalised eigenvector, and θ j is the jth parameter.Following the terminology of Monsalve-Bravo et al. ( 2022), we refer to (10) as the logarithm of the kth eigenparameter.The stiffest and sloppiest eigenparameters are those associated with the highest and lowest eigenvalues respectively.For all examples, the noise parameters φ are treated as nuisance parameters in the analysis of model sloppiness and are excluded when calculating the sensitivity matrix (Monsalve-Bravo et al., 2022).
Posterior predictive plots are used to assess the predictive performance of CW-IEKI relative to SMC.The posterior predictive distribution is given by which can be sampled by first sampling from the posterior distribution {θ, φ} * J ∼ p(θ, φ | y 1:T ), then sampling from the likelihood y * | {θ, φ} * J ∼ N (G(θ * J ), Γ(φ * J )).We compare the posterior predictive distribution estimated using the biased posterior samples from CW-IEKI to the posterior predictive distribution using the asymptotically exact samples from SMC.
Since SMC is asymptotically unbiased, it is always expected to outperform CW-IEKI in terms of accuracy and predictive performance.The main advantage of CW-IEKI is a significant speed-up in computation time compared to SMC.We assume that the expense of evaluating the function G(•), i.e. the deterministic mean of the likelihood function, dominates the computation time.SMC has N J j=1 M j + 1 evaluations of G(•), while CW-IEKI only has N J.The value of N is fixed for both methods, while J and M j , j = 1, . . ., J are adapted.In general,

Model Example 1: Predicting Nitrogen Mineralisation
The first model predicts cumulative nitrogen mineralisation, assuming a measurement error distributed according to (Vilas et al., 2021): for j = 1, . . ., T and r = 1, . . ., R, where T is the number of timepoints, R is the number of replicates per timepoint, and x t1 , . . ., x t T are deterministic predictions of cumulative nitrogen mineralisation from version 7.10 of the APSIM model (Holzworth et al., 2014) configured with soil water and nitrogen modules (Probert et al., 1998).The function G(•) has numerous parameters, most of which are fixed at measured values (Probert et al., 1998), apart from the model parameters we seek to obtain improved estimates for.Following the approach of Rammay et al. (2020), the model error is separated into two parts, where the first term (ζ r tj ) is known and accounts for measurement error, and the second term (σ) is unknown and accounts for all other sources of error such as model misspecification.At each timepoint t j , j = 1, . . ., T and replicate r = 1, . . ., R, ζ r tj is set to 4% of the observation y r tj (APHA and AWWA, 2012).
We consider two versions of this model.The first estimates three parameters (fbiom, finert, σ) and is the one considered in Vilas et al. (2021).For the second model, three additional model parameters are estimated (ef biom = ef hum, rd biom and rd hum) -in the first model these parameters are fixed at ef biom = ef hum = 0.4, rd biom = 0.0081 and rd hum = 0.00015.As a shorthand, in the present work we refer to these two models as the three parameter and six parameter APSIM models respectively.See Probert et al. (1998) for more detail about the model parameters and the values of the remaining parameters.The models are applied to data from Allen et al. (2019) measuring changes in inorganic nitrogen in soil from the Mackay Whitsundays region of North Queensland.The data is obtained from four 301 day laboratory incubations (i.e.R = 4).The second model is also fitted to a dataset simulated using θ = {fbiom, finert, ef biom = ef hum, rd biom, rd hum} = {0.1,0.6, 0.3, 0.0025, 0.0005} and φ = σ = 8.To enable simulation from the model, the known portion of the error (ζ r tj ) is set to 4% of the mean at time t j , i.e. for the synthetic dataset, ζ r tj = 0.04 • x tj for all r = 1, . . ., 4, matching the number of replicates in the data from Allen et al. (2019).

Three parameter APSIM model applied to the real data
Figures 1 and 2 show the marginal posterior densities of the parameters and the eigenparameters of the three parameter APSIM model applied to the real data.Figure 3 shows the posterior predictive densities using CW-IEKI and SMC.On this example, both CW-IEKI and SMC have very similar results for accuracy and predictive performance.However, CW-IEKI is almost 11 times faster than SMC with 5000 evaluations of G(•) compared to 54000 for SMC.(The number of evaluations of G(•) in our study is always a multiple of 1000 because our chosen ensemble sizes for both CW-IEKI and SMC are N = 1000.)

Six parameter APSIM model applied to the real data
Figure 4 shows the marginal posterior density plots for the six parameter APSIM model applied to the real data.Figures 5 and 6 show the marginal densities and eigenvectors of the three stiffest eigenparameters, and Figure 7 shows the posterior predictive distributions.Unlike the three parameter  10).The logarithm of the eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed red-orange) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.Based on the eigenvectors in Figure 6, the parameters ef biom = ef hum and rd hum do not contribute significantly to the model fit.Interestingly, the CW-IEKI marginal posteriors for these two parameters show the greatest bias compared to the SMC results.Overall, CW-IEKI gives a reasonably good fit for the model.It is also around 33 times faster than SMC with 9000 evaluations of G(•) compared to 301000.
Six parameter APSIM model applied to the simulated data

SMC CW-IEKI prior
Figure 5: Marginal posterior density plots of the natural logarithm of the three stiffest eigenparameters for the six parameter APSIM model applied to the real data.Note that the uncertainty parameter σ is excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation ( 10).The logarithm of the eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed red-orange) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.

SMC CW-IEKI prior
Figure 9: Marginal posterior density plots of the natural logarithm of the three stiffest eigenparameters for the six parameter APSIM model applied to the simulated data.Note that the uncertainty parameter σ is excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation ( 10).The logarithm of the eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed red-orange) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.

Model Example 2: Predicting Seagrass Decline
The second model predicts shoot density decline in seagrass due to cumulative stress from water temperature and light (Adams et al., 2020).The model takes, as input, light, temperature and time period of stress, and outputs photosynthesis rates and changes in shoot density over time.The model has 18 model parameters and 5 noise parameters.Several of these parameters have different values for specific instantaneous temperatures T and mean daily temperatures T (see Adams et al., 2020, for full model and parameter details).Uniform priors are used for all parameters.See Table 1 for the parameter units and prior bounds.
The model is calibrated to net photosynthesis data (Collier et al., 2018) and shoot density data (Collier et al., 2016) separately for three species of tropical seagrass from the Great Barrier Reef --Cymodocea serrulata, Halodule uninervis and Zostera muelleri.In the likelihood function for modeldata calibration it is assumed that measurement noise present in net photosynthesis observations at a given temperature T are normally distributed with standard deviation σ P (T ).Similarly, measurement noise in shoot density observations is assumed to be normally distributed with standard deviation σ S (albeit with some modifications to account for when observed shoot density declines to zero, see Appendix B of Adams et al., 2020 shoots/pot -0 100 (Zm), 20 (Cs) and 50 (Hu) σS shoots/pot -0 20 Table 1: Units and prior bounds of all 23 parameters of the seagrass model.The noise parameters are φ = {σ P (21), σ P (25), σ P (30), σ P (35), σ S }, and θ is comprised of the remaining parameters.In the final column, Cs, Hu and Zm refers to the species C. serrulata, H. uninervis and Z. muelleri respectively.
For brevity, all results shown in this section are for C. serrulata.Results for H. uninervis and Z. muelleri are provided in Appendix A. On this model, we also test the impact of the target ESS threshold on the accuracy of CW-IEKI. Figure 12 shows the marginal posterior densities of the parameters for SMC and CW-IEKI with the different ESS targets.For the majority of the parameters, the SMC and CW-IEKI densities are very similar.The target ESS threshold therefore appears to have little impact on the results.
As the parameters µ net,max (21.9) and µ net,max (27.9) are bounded between −0.02 and 0.02 (see Table 1), the log-transform cannot be used when performing the analysis of model sloppiness.Instead, we rescale all the model parameter to be between [0, 1] using the prior bounds, and then apply a logit transformation to map these values back to [−∞, ∞].The sensitivity matrix is given by the inverse of the covariance of the logit-transformed posterior samples from SMC, and the eigenparameters are given by where (v k ) j is the jth element of the kth normalised eigenvector, θ j is the jth parameter, a j is the prior lower bound of parameter j and b j is the prior upper bound of parameter j.Figures 13 and 14 show the marginal densities and eigenvectors of the six stiffest eigenparameters.The densities of these eigenparameters are similar for SMC and CW-IEKI, although again, the CW-IEKI results have greater uncertainty.Based on the eigenvectors in Figure 14, the parameters µ net,max (21.9), µ net,max (27.9),C other loss (21.9) and C other loss (27.9) do not significantly influence the model fit.As with the six parameter APSIM model, the CW-IEKI marginal posteriors for less influential parameters show the greatest bias.
Figure 15 shows the posterior predictive plots of the net carbon fixation using SMC and CW-IEKI with an ESS target of 50%, and Figure 16 shows the posterior predictive plots of the shoot density decline.
The predictive performance of SMC and CW-IEKI are fairly similar for this example, except that the CW-IEKI predictions have greater uncertainty.As with the density plots, there is little difference between the posterior predictive plots for an ESS target threshold of 50% and higher ESS targets (not shown).Table 2 shows the computation cost for SMC and CW-IEKI.For an ESS target of 50%, CW-IEKI is approximately 40 times faster than SMC.
The results suggest that CW-IEKI gives a good fit for predicting shoot density decline and carbon fixation for C. serrulata.CW-IEKI also gives a good fit for H. uninervis, but not for Z. muelleri (see Appendix A).The relatively poor fit for the latter may be a result of the likelihood not being strictly Gaussian due to the modifications that ensure the predicted shoot density remains greater than or equal to 0. As a result of these modifications, the likelihood function is close to Gaussian for higher shoot density values, but deviates strongly when the shoot density declines to 0, which is more often the case for Z. muelleri than for the other seagrass species.

Model Example 3: Predicting Coral Calcification Rates
The final model predicts coral calcification rates by simulating the transport and reaction of relevant chemical species and metabolic fluxes from seawater to the coral skeleton.It is assumed that there are two layers between the seawater and the coral skeleton: the coelenteron and the extracellular calcifying medium (ECM).
The main reactions considered are photosynthesis and respiration (seawater ↔ coelenteron), passive transport processes (seawater ↔ coelenteron ↔ ECM), membrane transport processes (coelenteron ↔ ECM) and aragonite precipitation and dissolution (ECM ↔ coral skeleton).The two membrane transport pumps modelled as part of the membrane transport processes are a Ca-ATPase pump and a bicarbonate anion transport (BAT) pump.
The reactions are modelled by a system of ordinary differential equations (ODEs), and measurement error is assumed to be Gaussian with standard deviation σ.The calcification rate predictions of the model are obtained from the steady state solution of the ODEs -these are compared to the data for calibration.There are a total of 21 unknown parameters which correspond to the passive transport processes, the membrane transport processes and the measurement error variance.Uniform priors are used for all parameters.Table 3 shows the parameter units and prior bounds.See Galli and Solidoro (2018) for more detail about the model and the values of the remaining parameters, and Vollert et al. ( 2022) for an application of SMC and analysis of model sloppiness to this model-data calibration problem. Figure 13: Marginal density plots of the six stiffest eigenparameters (calculated using equation ( 11)) for the seagrass model applied to the C. serrulata data.Note that the uncertainty parameters in φ are excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation ( 11).The eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.Figure 17 shows the marginal densities of the three stiffest eigenparameters, Figure 18 shows the posterior predictive distribution and Table 4 shows the computation cost for SMC and CW-IEKI.Due to the limited data available for this model, the marginal posterior densities for SMC and CW-IEKI are close to the prior (see Appendix B).In contrast, the marginal densities of the stiffest eigenparameters shown in Figure 17 are much more informative.As with previous examples, the predictive performance of SMC and CW-IEKI are similar, except that the CW-IEKI results have greater uncertainty (Figure 18).Again, changing the ESS target threshold for CW-IEKI makes very little difference to the results.For a target threshold of 50%, CW-IEKI is almost 24 times faster than SMC (Table 4).
Figure 16: Comparison of the C. serrulata data to the median and 95% central credible intervals for the posterior predictive distribution of shoot density obtained from the seagrass model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).SI = surface irradiance and "hot" or "cold" indicates the temperature conditions under which the seagrass data was collected (Adams et al., 2020).
Table 3: Units and prior bounds of all 21 parameters of the coral calcification model.For this model, the uncertainty parameter is φ = σ.10).The logarithm of the eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.

Discussion
In this paper, we have introduced and tested a new method (CW-IEKI) which extends the IEKI method of Iglesias et al. (2018) to the case where the covariance matrix has unknown elements φ.Our component-wise IEKI approach is completely analogous to likelihood tempering SMC, and is a useful alternative to both MCMC and SMC for static Bayesian models of the form N (G(θ), Γ(φ)), where θ and φ are unknown, and when G(θ) is expensive to compute.Note that CW-IEKI can also be applied when the covariance matrix is a function of both θ and φ, as is the case for the six parameter APSIM model applied to the simulated data in Section 4.2.That is, CW-IEKI can be applied to models of the form N (G(θ), Γ(θ, φ)), where θ | φ is updated using EKI and φ | θ is updated using MCMC.Even though the inference from CW-IEKI is only unbiased for models with a linear Gaussian likelihood and Gaussian prior, we find in practice that it provides reasonable inference even if the model's likelihood is non-linear Gaussian and its prior is non-Gaussian.Additionally, CW-IEKI generally requires much fewer evaluations of G(•) than MCMC or SMC.
We compared our method to SMC on three ecological models, all of which have a non-linear Gaussian likelihood and a non-Gaussian prior.The accuracy, predictive performance and computation time, the latter of which is measured by the number of evaluations of the function G(•), were used to assess the performance of our method relative to the unbiased solution from SMC.In the three parameter APSIM model and the coral model, the accuracy of CW-IEKI and SMC were similar, but for the remaining models there was clear bias in the marginal posteriors for some of the parameters.Based on the stiffest eigenparameters however, the model parameters showing the most bias also had little impact on the model fit.Across all models, CW-IEKI had relatively similar predictive performance to SMC -except that the uncertainty of the predictions was consistently overestimated -but advantageously required 11-40 times less evaluations of G(•).We also found that increasing the ESS target threshold for CW-IEKI made little difference to its accuracy and predictive performance.
In all of the examples we found that the point predictions from our novel CW-IEKI method was quite accurate, but the uncertainty intervals were inflated relative to SMC, especially when the number of parameters was increased.Therefore if highly accurate uncertainty quantification or parameter inferences are needed for a given application, then SMC or MCMC may be worth the wait if they are computationally feasible.If exact inferences are desired, CW-IEKI proposals could potentially be used to speed up exact SMC, for example by incorporating them in the delayed-acceptance SMC algorithm of Bon et al. (2021).The inferences from CW-IEKI could also potentially be improved by following the approach of Lan et al. (2022) to build an emulator G(•) of G(•) using all evaluations of G(•) from CW-IEKI.An MCMC or SMC algorithm can then be used to target the approximate posterior distribution based on this emulator, i.e.N (y | G(θ), Γ(φ))p(θ, φ).
An area of future work is to improve the updates in CW-IEKI for the noise parameters φ.Currently, a fixed number of random-walk MCMC iterations are used.Adapting the number of MCMC iterations and using more efficient updates for φ, such as the Metropolis-adjusted Langevin algorithm (Girolami and Calderhead, 2011) or Hamiltonian Monte Carlo (Betancourt, 2017), may improve the performance of the method, especially if φ is high-dimensional or its elements are highly correlated.Another extension is to apply the CW-IEKI method to the hierarchical setting explored in Chada et al. (2018).
Another avenue of future work is to investigate how our approach can be incorporated into the IEKI method of Duffield and Singh (2021) for general likelihoods.In this case, it may be possible to update some of the model parameters with IEKI and some with MCMC, depending on the form of the likelihood function.The potential advantage of such a hybrid approach is that it may efficiently improve the accuracy of the final samples, given that the MCMC update targets the exact conditional posterior, while the IEKI portion targets some approximation to the conditional posterior.
It would also be interesting to incorporate the CW-IEKI method into the data annealing SMC algorithm of Wu et al. ( 2022), which currently requires the covariance of the likelihood function to be known.In particular, one extension here is developing a likelihood tempering SMC algorithm with our CW-IEKI method as the forward kernel.It may also be possible to extend their SMC algorithm to general likelihood models, such that a subset of the parameters are updated using the method of Duffield and Singh (2021), and the rest are updated using an MCMC forward kernel.
A Extra results for the seagrass model   11)) for the seagrass model applied to the H. uninervis data.Note that the uncertainty parameters in φ are excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation ( 11).The eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.11)) for the seagrass model applied to the Z. muelleri data.Note that the uncertainty parameters in φ are excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation ( 11).The eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.Figure 28: Comparison of the Z. muelleri data to the median and 95% central credible intervals for the posterior predictive distribution of shoot density obtained from the seagrass model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).SI = surface irradiance and "hot" or "cold" indicates the temperature conditions under which the seagrass data was collected (Adams et al., 2020

Figure 1 :Figure 2 :
Figure 1: Marginal posterior density plots for the three parameter APSIM model applied to the real data.

Figure 3 :
Figure 3: Comparison of the real data to the median and 95% central credible intervals for the posterior predictive distribution of cumulative nitrogen mineralisation obtained from the three parameter APSIM model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).

Figure 8
Figure8shows the marginal posterior densities of the six parameter APSIM model applied to the simulated data.As before, SMC and CW-IEKI have similar results, except that CW-IEKI has posterior support for larger values of σ.Thus, this simulation demonstrates that larger support for σ from CW-IEKI is not an artefact of model misspecification, as the data used here is simulated from the six parameter APSIM model.Figures9 and 10show the densities and eigenvectors of the three stiffest eigenparameters respectively.The eigenparameter densities are very similar for SMC and CW-IEKI, indicating that CW-IEKI gives a relatively good fit for the model, and the eigenvectors again show that ef biom = ef hum and rd hum have little influence on the model fit.The posterior predictive distribution in Figure11also shows similar performance between SMC and CW-IEKI, except that CW-IEKI has much greater uncertainty.On this example, CW-IEKI is around 37 times faster than SMC with 11000 evaluations of G(•) compared to 411000.

Figure 4 :
Figure 4: Marginal posterior density plots for the six parameter APSIM model applied to the real data.

Figure 6 :
Figure 6: Eigenvectors of the three stiffest eigenparameters for the six parameter APSIM model applied to the real data.These results are based on the SMC posterior samples.The labels on the left-hand side correspond to v k (λ k /λ 1 ), where v k and λ k are the eigenvector and associated eigenvalue of eigenparameter k, respectively.The shade of the cells in row k indicate the relative contribution (v k ) j of the jth parameter to eigenparameter k -parameters with darker colours have the greatest contribution.

Figure 7 :
Figure 7: Comparison of the real data to the median and 95% central credible intervals for the posterior predictive distribution of cumulative nitrogen mineralisation obtained from the six parameter APSIM model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).

Figure 8 :
Figure 8: Marginal posterior density plots for the six parameter APSIM model applied to the simulated data.

Figure 10 :
Figure 10: Eigenvectors of the three stiffest eigenparameters for the six parameter APSIM model applied to the simulated data.These results are based on the SMC posterior samples.The labels on the left-hand side correspond to v k (λ k /λ 1 ), where v k and λ k are the eigenvector and associated eigenvalue of eigenparameter k, respectively.The shade of the cells in row k indicate the relative contribution (v k ) j of the jth parameter to eigenparameter k -parameters with darker colours have the greatest contribution.

Figure 11 :
Figure 11: Comparison of the simulated data to the median and 95% central credible intervals for the posterior predictive distribution of cumulative nitrogen mineralisation obtained from the six parameter APSIM model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).

Figure 12 :
Figure 12: Marginal posterior density plots for the seagrass model applied to the C. serrulata data.

Figure 14 :
Figure 14: Eigenvectors of the six stiffest eigenparameters for the seagrass model applied to the C. serrulata data.These results are based on the SMC posterior samples.The labels on the lefthand side correspond to v k (λ k /λ 1 ), where v k and λ k are the eigenvector and associated eigenvalue of eigenparameter k, respectively.The shade of the cells in row k indicate the relative contribution (v k ) j of the jth parameter to eigenparameter k -parameters with darker colours have the greatest contribution.

Figure 15 :
Figure 15: Comparison of the C. serrulata data to the median and 95% central credible intervals for the posterior predictive distribution of net carbon fixation obtained from the seagrass model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).P-I = photosynthesis-irradiance and PAR = photosynthetically active radiation (see Adams et al., 2020).

Figure 17 :
Figure17: Marginal posterior density plots of the natural logarithm of the three stiffest eigenparameters for the coral model.CW-IEKI results use different ESS target thresholds.Note that the uncertainty parameter σ is excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation (10).The logarithm of the eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.

Figure 18 :
Figure 18: Comparison of the coral data to the median and 95% central credible intervals for the posterior predictive distribution of coral calcification obtained from the coral model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).
Figure 19: Marginal posterior density plots for the seagrass model applied to the H. uninervis data.

Figure 20 :
Figure 20: Marginal density plots of the six stiffest eigenparameters (calculated using equation (11)) for the seagrass model applied to the H. uninervis data.Note that the uncertainty parameters in φ are excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation (11).The eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.

Figure 21 :
Figure 21: Eigenvectors of the six stiffest eigenparameters for the seagrass model applied to the H. uninervis data.These results are based on the SMC posterior samples.The labels on the lefthand side correspond to v k (λ k /λ 1 ), where v k and λ k are the eigenvector and associated eigenvalue of eigenparameter k, respectively.The shade of the cells in row k indicate the relative contribution (v k ) j of the jth parameter to eigenparameter k -parameters with darker colours have the greatest contribution.

Figure 22 :
Figure22: Comparison of the H. uninervis data to the median and 95% central credible intervals for the posterior predictive distribution of net carbon fixation obtained from the seagrass model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).P-I = photosynthesis-irradiance and PAR = photosynthetically active radiation (seeAdams et al., 2020).

Figure 23 :
Figure23: Comparison of the H. uninervis data to the median and 95% central credible intervals for the posterior predictive distribution of shoot density obtained from the seagrass model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).SI = surface irradiance and "hot" or "cold" indicates the temperature conditions under which the seagrass data was collected(Adams et al., 2020).

Figure 24 :
Figure 24: Marginal posterior density plots for the seagrass model applied to the Z. muelleri data.

Figure 25 :
Figure25: Marginal density plots of the six stiffest eigenparameters (calculated using equation (11)) for the seagrass model applied to the Z. muelleri data.Note that the uncertainty parameters in φ are excluded from the analysis of sloppiness, and λ k is the eigenvalue associated with eigenvector v k in equation (11).The eigenparameters are calculated based on samples from the prior (black), CW-IEKI (dashed) and SMC (blue) using a sensitivity matrix calculated using the SMC samples.

Figure 26 :
Figure 26: Eigenvectors of the six stiffest eigenparameters for the seagrass model applied to the Z. muelleri data.These results are based on the SMC posterior samples.The labels on the lefthand side correspond to v k (λ k /λ 1 ), where v k and λ k are the eigenvector and associated eigenvalue of eigenparameter k, respectively.The shade of the cells in row k indicate the relative contribution (v k ) j of the jth parameter to eigenparameter k -parameters with darker colours have the greatest contribution.

Figure 27 :
Figure 27: Comparison of the Z. muelleri data to the median and 95% central credible intervals for the posterior predictive distribution of net carbon fixation obtained from the seagrass model fitted to this data.Models were fitted using CW-IEKI (red-orange) and likelihood tempering SMC (blue).P-I = photosynthesis-irradiance and PAR = photosynthetically active radiation (see Adams et al., 2020).

Figure 29 :
Figure 29: Marginal posterior density plots for the coral model.
for further details).

Table 2 :
Total and relative number of evaluations of G(•) for SMC and CW-IEKI with different ESS target thresholds.Results are for the seagrass model applied to the C. serrulata data.Note that the total number of evaluations of G(•) is a multiple of the number of samples N = 1000.

Table 4 :
Total and relative number of evaluations of G(•) for SMC and CW-IEKI with different ESS target thresholds for the coral model.Note that the total number of evaluations of G(•) is a multiple of the number of samples N = 1000.

Table 5 :
).Total and relative number of evaluations of G(•) for SMC and CW-IEKI with different ESS target thresholds.Results are for the seagrass model applied to the H. uninervis data.Note that the total number of evaluations of G(•) is a multiple of the number of samples N = 1000.

Table 6 :
).Total and relative number of evaluations of G(•) for SMC and CW-IEKI with different ESS target thresholds.Results are for the seagrass model applied to the Z. muelleri data.Note that the total number of evaluations of G(•) is a multiple of the number of samples N = 1000.