Joint state-parameter estimation of a nonlinear stochastic energy balance model from sparse noisy data

TL;DR

This paper develops a regularized Bayesian approach for joint state-parameter estimation in nonlinear stochastic energy balance models, effectively handling sparse noisy data and ill-posedness to improve inference accuracy in paleoclimate studies.

Contribution

It introduces a strongly regularized posterior with physical priors and a particle Gibbs sampler tailored for high-dimensional SEBM, enhancing estimation robustness and uncertainty quantification.

Findings

Regularized posteriors overcome ill-posedness and stay within physical ranges.

Posterior of states concentrates near the true states, reducing noise impact.

Parameter posteriors exhibit large uncertainty due to ill-posedness and regularization.

Abstract

While nonlinear stochastic partial differential equations arise naturally in spatiotemporal modeling, inference for such systems often faces two major challenges: sparse noisy data and ill-posedness of the inverse problem of parameter estimation. To overcome the challenges, we introduce a strongly regularized posterior by normalizing the likelihood and by imposing physical constraints through priors of the parameters and states. We investigate joint parameter-state estimation by the regularized posterior in a physically motivated nonlinear stochastic energy balance model (SEBM) for paleoclimate reconstruction. The high-dimensional posterior is sampled by a particle Gibbs sampler that combines MCMC with an optimal particle filter exploiting the structure of the SEBM. In tests using either Gaussian or uniform priors based on the physical range of parameters, the regularized posteriors overcome the ill-posedness and lead to samples within physical ranges, quantifying the uncertainty in estimation. Due to the ill-posedness and the regularization, the posterior of parameters presents a relatively large uncertainty, and consequently, the maximum of the posterior, which is the minimizer in a variational approach, can have a large variation. In contrast, the posterior of states generally concentrates near the truth, substantially filtering out observation noise and reducing uncertainty in the unconstrained SEBM.

Tables

Table 1. Table 1 : The physical upper and lower bounds of the parameters in the SEBM.

	${\theta }_0$	${\theta }_1$	${\theta }_4$
upper bound	32.57	-22.70	-4.80
lower bound	27.64	-25.46	-6.00

Table 2. Table 2 : The priors of θ = ( θ 0 , θ 1 , θ 4 ) 𝜃 subscript 𝜃 0 subscript 𝜃 1 subscript 𝜃 4 \theta=(\theta_{0},\theta_{1},\theta_{4}) based on the physical constraints in Table 1 .

Uniform prior	[27.64, 32.57]$\times$[-25.46, -22.70]$\times$[-6.00, -4.80]
Gaussian prior	$\mathrm{mean}=(30.11,-24.08,-5.40)$
	$\mathrm{covariance}=\mathrm{Diag}$(${0.82}^2$, ${0.46}^2$, ${0.20}^2$)

Table 3. Table 3 : The settings of the stochastic energy balance model and its discretization.

$\nu$	$=0.1$	Diffusion constant
${\sigma }_f$	$=0.1$	Scale of the stochastic forcing
$\Delta t$	$=0.01$	Time step size
$d_b$	$=12$	Number of total nodes
$d_o$	$=6$	Number of observed nodes
${\sigma }_{ϵ}$	$=0.01$	Std of the observation noise

Table 4. Table 4 : The settings of the particle MCMC using SMC with optimal importance densities.

$M$	$=5$	Number particles in SMC
$L$	$={10}^4$	Length of the Markov chain
$N$	$=100$	Number of time steps of observations.

Equations

$${\partial }_{t}u(t,\xi )-\nu \Delta u(t,\xi )=g_{\theta }(u)+f(t,\xi ),$$

$$g_{\theta }(u)={\theta }_0+{\theta }_{1}u+{\theta }_4{u}^4,$$

$$\mathbb{E}\left[f(t,\xi )f(s,\eta )\right]=\delta (t-s)C(|\xi -\eta |),$$

$$C_{\alpha }(r)={\sigma }_f^2\frac{2^{1-\alpha }}{\Gamma (\alpha )}{\left(\sqrt{2\alpha }\frac{r}{\rho }\right)}^{\alpha }{K}_{\alpha }\left(\sqrt{2\alpha }\frac{r}{\rho }\right),$$

$$y_i(t)=H_i(u(t))+{ϵ}_i(t)=u(t,{\xi }_i)+{ϵ}_i(t),$$

$$U_{n+1}={\mu }_{\theta }(U_n)+W_n$$

$$p_{\theta }(u_{n+1}|u_n)=\det (2\pi \mathbf{R}{)}^{-1/2}\exp \left(-\frac{(u_{n+1}-{\mu }_{\theta }(u_n){)}^T{\mathbf{R}}^{-1}(u_{n+1}-{\mu }_{\theta }(u_n))}{2}\right).$$

$$Y_n=\mathbf{H}{U}_n+{ϵ}_n,y_n\in {\mathbb{R}}^{d_o},$$

$$p(y_n|U_n)=\det (2\pi \mathbf{Q}{)}^{-1/2}\exp \left(-\frac{(y_n-\mathbf{H}{U}_n{)}^T{\mathbf{Q}}^{-1}(y_n-\mathbf{H}{U}_n)}{2}\right).$$

$$p(\theta ,u_{1:N}|y_{1:N})=p(\theta )\frac{p_{\theta }(u_{1:N})p_{\theta }(y_{1:N}|u_{1:N})}{p_{\theta }(y_{1:N})},$$

$$p_{\theta }(y_{1:N}|u_{1:N})=p(y_{1:N}|u_{1:N})=\prod _{n}p(y_n|u_n),$$

$$p_{\theta }(u_{1:N})=p_{\theta }(u_1)\prod _{n=1}^{N-1}{p}_{\theta }(u_{n+1}|u_n),$$

$$p(u_{1:N}|\theta ,y_{1:N})=\frac{p_{\theta }(u_{1:N})p_{\theta }(y_{1:N}|u_{1:N})}{p_{\theta }(y_{1:N})},$$

$$p(\theta |u_{1:N},y_{1:N})=p(\theta |u_{1:N})=p(\theta )p_{\theta }(u_{1:N}).$$

$$l(\theta ,u_{1:N})=c-\frac{1}{2}\sum _{n=1}^N(u_{n+1}-{\mu }_{\theta }(u_n){)}^T{\mathbf{R}}^{-1}(u_{n+1}-{\mu }_{\theta }(u_n)),$$

$${\mathbf{F}}_N=\frac{1}{N}{\left(\sum _{n=1}^N{G}_{\theta ,k}(u_n{)}^T{\mathbf{R}}^{-1}{G}_{\theta ,l}(u_n)\right)}_{k,l=0,1,4},$$

$$p^N(\theta ,u_{1:N}|y_{1:N})=\frac{1}{Z}p(\theta ){\left[\frac{p^c(u_{1:N})p_{\theta }(u_{1:N})p_{\theta }(y_{1:N}|u_{1:N})}{p_{\theta }(y_{1:N})}\right]}^{1/N},$$

$$p^c(u_{1:N}):=\prod _{i=1}^N\frac{1}{2\pi {\sigma }_c^{d_b}}\exp \left(-\frac{|u_i-u_c{|}^2}{2{\sigma }_c^2}\right)$$

$$-\log p^N(\theta ,u_{1:N}|y_{1:N})=\frac{1}{N}{C}_{y_{{}_{1:N}}}(\theta ,u_{1:N})+\log Z-\frac{1}{N}\log p_{\theta }(y_{1:N}),$$

$$C_{y_{{}_{1:N}}}(\theta,u_{1:N})\;=\;-\sum_{n=1}^{N}\log\big{[}p(u_{n}|u_{n-1},\theta)p(y_{n}|u_{n})\big{]}-N\log p(\theta)-\log p^{c}(u_{1:N}).$$

$$p^N(\theta |u_{1:N},y_{1:N})\propto p(\theta )[p_{\theta }(u_{1:N}){]}^{1/N}.$$

$$u_{d_b}(t,\xi )=\sum _{i=1}^{d_b}{\hat{u}}_i(t){\varphi }_i(\xi ).$$

$$⟨u_{d_b}(t,\cdot ),\varphi ⟩=⟨u_0,\varphi ⟩-\nu {\int }_0^t⟨\nabla u_{d_b}(s,\cdot ),\nabla \varphi ⟩ds+{\int }_0^t⟨g_{\theta }(u_{d_b}(s,\cdot )),\varphi ⟩ds+{\int }_0^t⟨f(s,\cdot ),\varphi ⟩,$$

$$U(t)$$

$$\Phi (\xi )$$

$$u_{d_b}(t,\xi )$$

$$⟨\Phi ,{\Phi }^T⟩U(t)=⟨\Phi ,{\Phi }^T⟩U(0)-\nu ⟨\nabla \Phi ,\nabla {\Phi }^T⟩{\int }_0^{t}U(s)ds+{\int }_0^t⟨g_{\theta }(U_n^T\Phi ),\Phi ⟩ds+{\int }_0^t⟨f(s,\cdot ),\Phi ⟩.$$

$${\mathbf{M}}_0=⟨\Phi ,{\Phi }^T⟩,{\mathbf{M}}_1=\nu ⟨\nabla \Phi ,\nabla {\Phi }^T⟩,$$

$$G_{\theta }(U(t)):=⟨g_{\theta }(U^T(t)\Phi ),\Phi ⟩.$$

$${\mathbf{M}}_{0}U(t)={\mathbf{M}}_{0}U(0)-{\mathbf{M}}_1{\int }_0^{t}U(s)ds+{\int }_0^t{G}_{\theta }(U(t))ds+{\int }_0^t⟨f(s,\cdot ),\Phi ⟩.$$

$$\int g_{\theta }(u(t,\xi )){\varphi }_l(\xi )d\xi \approx \sum _{{\mathcal{T}}_k\subset \mathrm{supp}({\varphi }_l)}\frac{\mathrm{Area}({\mathcal{T}}_k)}{3}{g}_{\theta }\left(\sum _i{U}_i(t){\varphi }_i({\xi }_k^c)\right),$$

Full text

We investigate joint parameter-state estimation by the regularized posterior in a physically motivated nonlinear stochastic energy balance model (SEBM) for paleoclimate reconstruction. The high-dimensional posterior is sampled by a particle Gibbs sampler that combines MCMC with an optimal particle filter exploiting the structure of the SEBM. In tests using either Gaussian or uniform priors based on the physical range of parameters, the regularized posteriors overcome the ill-posedness and lead to samples within physical ranges, quantifying the uncertainty in estimation. Due to the ill-posedness and the regularization, the posterior of parameters presents a relatively large uncertainty, and consequently, the maximum of the posterior, which is the minimizer in a variational approach, can have a large variation. In contrast, the posterior of states generally concentrates near the truth, substantially filtering out observation noise and reducing uncertainty in the unconstrained SEBM.

1. 1 Introduction

2. 2 State-space model formulation

3. 2.1 The stochastic energy balance model

4. 2.2 State-space model representation

5. 2.2.1 The state model

6. 2.2.2 The observation model

7. 2.3 Bayesian inference for SSM

8. 2.4 Sampling the posterior by particle MCMC methods

9. 3 Ill-posedness and regularized posteriors

10. 3.1 Model settings and tests

11. 3.2 Ill-posedness of the standard Bayesian inference of parameters

12. 3.3 Regularized posteriors

13. 4 Bayesian inference with regularized posteriors

14. 4.1 Diagnosis of the Markov Chain

15. 4.2 Parameter estimation

16. 4.3 State estimates

17. 5 Discussion

18. 5.1 Observing fewer nodes

19. 5.2 Observing a longer trajectory.

20. 5.3 Estimates of the nonlinear function

21. 6 Conclusions and future work

22. 7 Appendix: Technical details of the estimation procedure

23. 7.1 Discretization of the SEBM

24. 7.2 SMC with optimal importance sampling

25. 7.3 Particle Gibbs and PGAS

1 Introduction

Physically motivated nonlinear stochastic (partial) differential equations (SDEs and SPDEs) are natural models of spatiotemporal processes with uncertainty in geoscience. In particular, such models arise in the problem of reconstructing geophysical fields from sparse and noisy data (see e.g. [42, 20, 45] and the references therein). The nonlinear differential equations, derived from physical principles, often come with unknown but physically constrained parameters also to be determined from data. This promotes the problem of joint state-parameter estimation from sparse and noisy data. When the parameters are interrelated, which is often the case in nonlinear models, their estimation can be an ill-posed inverse problem. Physical constraints on the parameters must then be taken into account. In variational approaches, physical constraints are imposed using a regularization term in a cost function, whose minimizer provides an estimator of the parameters and states. In a Bayesian approach, the physical constraints are encoded in prior distributions, extending the regularized cost function in the variational approach to a posterior and quantifying the estimation uncertainty. When the true parameters are known, the Bayesian approach has demonstrated great success in state estimation, thanks to the developments in Monte Carlo sampling and data assimilation techniques (see e.g. [8, 25, 49]). However, the problem of joint state-parameter estimation, especially when the parameter estimation is ill-posed, has had relatively little success in nonlinear cases and remains a challenge [23].

In this paper, we investigate a Bayesian approach for joint state and parameter estimation of a non-linear two-dimensional stochastic energy balance model (SEBM) in the context of spatial-temporal paleoclimate reconstructions of temperature fields from sparse and noisy data [46, 44, 16, 19]. In particular, we consider a model of the energy balance of the atmosphere similar to those often used in idealized climate models (see e.g. [17, 50, 41]) to study climate variability and climate sensitivity. The use of such a model in paleoclimate reconstruction aims at improving the physical consistency of temperature reconstructions during e.g. the last deglaciation and the Holocene by combining indirect observations, so called proxy data, with physically-motivated stochastic models.

The SEBM models surface air temperature, explicitly taking into account sinks, sources, and horizontal transport of energy in the atmosphere, with an additive stochastic forcing incorporated to account for unresolved processes and scales. The model takes the form of a nonlinear SPDE with unknown parameters to be inferred from data. These unknown parameters are associated with processes in the energy budget (e.g. radiative transfer, air-sea energy exchange) that are represented in a simplified manner in the SEBM, and may change with a changing climate. The parameters must fall in a prescribed range such that the SEBM is physically meaningful. Specifically, they must be in sufficiently close balance for the stationary temperature of the SEBM to be within a physically realistic range. As we will show, the parametric terms arising from this physically-based model are strongly correlated, leading to a Fisher information matrix that is ill-conditioned. Therefore, the parameter estimation is an ill-posed inverse problem, and the maximum likelihood estimators of individual parameters have large variations and often fall out of the physical range.

To overcome the ill-posedness in parameter estimation, we introduce a new strongly regularized posterior by normalizing the likelihood and by imposing the physical constraints through priors on the parameters and the states, based on physical constraints and the climatological distribution. In the regularized posterior, the prior has the same weight as the normalized likelihood to enforce the support of the posterior to be in the physical range. Such a regularized posterior is a natural extension of the regularized cost function in a variational approach: the maximum of the posterior (MAP) is the same as the minimizer of the regularized cost function, but the posterior quantifies the uncertainty in the estimator.

The regularized posterior of the states and parameters is high-dimensional and non-Gaussian. It is represented by its samples, which provide an empirical approximation of the distribution and allow efficient computation of quantities of interest such as posterior means. The samples are drawn using a particle Gibbs sampler with ancestor sampling (PGAS ) [28], a special sampler in the family of particle Markov chain Monte Carlo (MCMC) methods [2] that combines the strengths of both MCMC and sequential Monte Carlo methods (see e.g. [15]) to ensure the convergence of the empirical approximation to the high-dimensional posterior. In the PGAS, we use an optimal particle filter that exploits the forward structure of the SEBM.

We consider two priors for the parameters, each based on their physical ranges: a uniform prior and a Gaussian prior with three standard deviations inside the range. We impose a prior for the states based on their overall climatological distribution. Tests show that the regularized posteriors overcome the ill-posedness and lead to samples of parameters and states within the physical ranges, quantifying the uncertainty in their estimation. Due to the regularization, the posterior of the parameters is supported on a relatively large range. Consequently, the MAP of the parameters has a large variation, and it is important to use the posterior to assess the uncertainty. In contrast, the posterior of the states generally concentrates near the truth, substantially filtering out the observational noise and reducing the uncertainty in state reconstruction.

Tests also show that the regularized posterior is robust to spatial sparsity of observations, with sparser observations leading to larger uncertainties. However, due to the need for regularization to overcome ill-posedness, the uncertainty in the posterior of the parameters can not be eliminated by increasing the number of observations in time. Therefore, we suggest alternative approaches, such as re-parametrization of the nonlinear function according to the climatological distribution or nonparametric Bayesian inference (see e.g. [36, 18]), to avoid ill-posedness.

The rest of the paper is organized as follows. Section 2 introduces the SEBM and its discretization, and formulates a state-space model. We also outline in this section the Bayesian approach to the joint parameter-state estimation and the particle MCMC samplers. Section 3 analyzes the ill-posedness of the parameter estimation problem and introduces the regularized posterior. The regularized posterior is sampled by PGAS and numerical results are presented in Section 4. Discussions and conclusions are presented in Sections 5 and 6. Technical details of the estimation procedure are described in Appendix 7.

2 State-space model formulation

After providing a brief physical introduction to the SEBM, we present its discretization and the observation model by representing them as a state-space model suitable for application of sequential Monte Carlo methods in Bayesian inference.

2.1 The stochastic energy balance model

The SEBM describes the evolution in space (both latitude and longitude) and time of the surface air temperature $u(t,\xi)$:

[TABLE]

where $\xi\in[-\pi,\pi]\times[-\pi/2,\pi/2]$ is the two-dimensional coordinate on the sphere and the solution $u(t,\xi)$ is periodic in longitude. Horizontal energy transport is represented as diffusion with diffusivity $\nu$, while sources and sinks of atmospheric internal energy are represented by the nonlinear function $g_{\theta}(u)$

[TABLE]

with the unknown parameters $\theta$. Upper and lower bounds of these three parameters, shown in Table 1, are derived from the energy balance model in [17], adjusted to current estimates of the Earth’s global energy budget from [47] using appropriate simplifications. The equilibrium solution of the SEBM for the average values of the parameters approximates the current global mean temperature closely, and the magnitude of sinks and sources approximates the respective magnitudes in [47] well. The physical ranges of the parameters are very conservative and cover current estimates of the global mean temperature during the Quaternary [43]. The state variable and the parameters in the model have been nondimensionalized so that the equilibrium solution of Eqn. (2.1) with $f=0$ is approximately equal to one.

The quartic nonlinearity of the function $g_{\theta}(u)$ arises from the Stefan-Boltzmann dependence of long-wave radiative fluxes on atmospheric temperature, while a linear feedback is included to represent state dependence of e.g. surface energy fluxes and albedo. Inclusion of quadratic and cubic nonlinarities in $g_{\theta}(u)$ (to account for nonlinearities in the feedbacks just noted) was found to exacerbate the ill-posedness of the model without qualitatively changing the character of the model dynamics within the parameter range appropriate for the study of Quaternary climate variability (e.g. without admitting multiple deterministic equilibria associated with the ice-albedo feedback). In reality, the diffusivity $\nu$ and the parameters $\theta_{j}$, $j=(0,1,4)$ will depend on latitude, longitude, and time. We will neglect this complexity in our idealized analysis.

The stochastic term $f(t,\xi)$, which models the net effect of unresolved or oversimplified processes in the energy budget, is a centered Gaussian field that is white in time and colored in space, specified by an isotropic Matérn covariance function with order $\alpha=1$ and scale $\rho>0$. That is,

[TABLE]

with the covariance kernel $C(r)$ being the Matérn covariance kernel given by

[TABLE]

where $\Gamma$ is the gamma function, $\rho$ is a scaling factor, and $K_{\alpha}$ is the modified Bessel function of the second kind. We focus on the estimation of the parameters $\theta$ and assume that $\nu$ and the parameters of $f$ are known. Estimating $\nu$ in energy balance models with data assimilation methods is studied in [3], whereas estimation of parameters of $f$ in the context of linear SPDEs is covered for example in [26].

In a paleoclimate context, temperature observations are sparse (in space and time) and derived from climatic proxies, such as pollen assemblages, isotopic compositions, and tree rings, that are indirect measures of the climate state. To simplify our analysis, we neglect the potentially nonlinear transformations associated with the proxies and focus on the effect of observational sparseness. This is a common strategy in the testing of climate field reconstruction methods [51]. As such, we take the data to be noisy observations of the solution at $d_{o}$ locations:

[TABLE]

for $i=1,\dots,d_{o}$, where each $\xi_{i}\in[-\pi,\pi]\times[-\pi/2,\pi/2]$ is a location of observation, $H$ is the observation operator, and $\epsilon_{i}(t)\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$ are iid Gaussian noise. The data are sparse in the sense that only a small number of the spatial locations are observed.

2.2 State-space model representation

In practice, the differential equations are represented by their discretized systems and the observations are discrete in time, therefore we consider only the state space model based on a discretization of the SEBM. We refer the reader to [40, 4, 21, 34, 30] for studies about inference of SPDEs in a continuous-time setting.

2.2.1 The state model

We discretize the SPDE (2.1) using linear finite elements in space and a semi-backward Euler method in time, using the computationally efficient Gaussian Markov random field approximation of the Gaussian field by [26] (see details in Section 7.1). We write the discretized equation as a standard state space model:

[TABLE]

where $\mu_{\theta}:\mathbb{R}^{d_{b}}\to\mathbb{R}^{d_{b}}$ is the deterministic function and $\{W_{n}\}$ is a sequence of iid Gaussian noise with mean zero and covariance $\mathbf{R}$ described in more detail in Section (7.19). Therefore, the transition probability density $p_{\theta}(u_{n+1}|u_{n})$, the probability density of $U_{n+1}$ conditional on $U_{n}$ and $\theta$, is

[TABLE]

2.2.2 The observation model

In discrete form, we assume that the locations of observation are the nodes of the finite elements. Then the observation function in (2.5) is simply $H_{i}(U_{n})=U_{n,k_{i}}$ with $k_{i}\in\{1,\dots,d\}$ denoting the index of the node under observation, for $i=1,\dots,d_{0}$, and we can write the observation model as

[TABLE]

where $\mathbf{H}\in\mathbb{R}^{d_{o}\times d_{b}}$ is called the observation matrix, and $\{\epsilon_{n}\}$ is a sequence of iid Gaussian noise with distribution $\mathcal{N}(0,\mathbf{Q})$, where $\mathbf{Q}=\textrm{Diag}\{\sigma_{i}^{2}\}$. Equivalently, the probability of observing $y_{n}$ given state $U_{n}$ is

[TABLE]

2.3 Bayesian inference for SSM

Given observations $y_{1:N}:=(y_{1},\dots,y_{N})$, our goal is to jointly estimate the state $U_{1:N}:=(U_{1},\dots,U_{N})$ and the parameter vector $\theta:=(\theta_{0},\theta_{1},\theta_{4})$ in the state-space model (2.6)-(2.9). The Bayesian approach estimates the joint distribution of $(U_{1:N},\theta)$ conditional on the observations by drawing samples to form an empirical approximation of the high-dimensional posterior. The empirical posterior efficiently quantifies the uncertainty in the estimation. Therefore, the Bayesian approach has been widely used (see the review [23] and the references therein).

Following Bayes’ rule, the joint posterior distribution of $(U_{1:N},\theta)$ can be written as

[TABLE]

where $p(\theta)$ is the prior of the parameters and $p_{\theta}(y_{1:N})=\int p_{\theta}(u_{1:N})p_{\theta}(y_{1:N}|u_{1:N})du_{1:N}$ is the unknown marginal probability density function of the observations. In the importance sampling approximation to the posterior, we do not need to know the value of $p_{\theta}(y_{1:N})$, because as a normalizing constant, and it will be cancelled out in the importance weights of samples. The quantity $p_{\theta}(y_{1:N}|u_{1:N})$ is the likelihood of the observations $y_{1:N}$ conditional on the state $U_{1:N}$ and the parameter $\theta$, which can be explicitly derived from the observation model (2.8):

[TABLE]

with $p(y_{n}|u_{n})$ given in (2.9). Finally, the probability density function of the state $U_{1:N}$ given parameter $\theta$ can be derived from the state model (2.6):

[TABLE]

with $p_{\theta}(u_{n+1}|u_{n})$ specified by (2.7).

2.4 Sampling the posterior by particle MCMC methods

In practice, we are interested in the expectation of quantities of interest or the probability of certain events. These computations involve integrations of the posterior that can neither be computed analytically nor by numerical quadrature methods due to the curse of dimensionality: the posterior is a high-dimensional non-Gaussian distribution involving variables with a dimension at the scale of thousands to millions. Monte Carlo methods generate samples to approximate the posterior by the empirical distribution, so that quantities of interest can be computed efficiently.

Markov Chain Monte Carlo (MCMC) methods are popular Monte Carlo methods (see e.g. [29]) that generate samples along a Markov chain with the posterior as the invariant measure. For joint distributions of parameters and states, a standard MCMC method is Gibbs sampling which consists of alternatively updating the state variable $U_{1:N}$ conditional on $\theta$ and $y_{1:N}$ by sampling

[TABLE]

and then updating the parameter $\theta$ conditional on $U_{1:N}=u_{1:N}$ by sampling the marginal posterior of $\theta$:

[TABLE]

Due to the high-dimensionality of $U_{1:N}$, a major difficulty in sampling $p(u_{1:N}|\theta,y_{1:N})$ is the design of efficient proposal densities that can effectively explore the support of $p(u_{1:N}|\theta,y_{1:N})$.

Another group of rapidly-developing MC methods are sequential Monte Carlo (SMC) methods [7, 15] that exploit the sequential structure of state space models to approximate the posterior densities $p(u_{1:n}|\theta,y_{1:N})$ sequentially. SMC methods are efficient but suffer from the well-known problem of depletion (or degeneracy), in which the marginal distribution $p(u_{n}|\theta,y_{1:N})$ becomes concentrated on a single sample as $N-n$ increases (see Section 7.2 for more details).

The particle MCMC methods introduced in [2] provide a framework for systematically combining SMC methods with MCMC methods, exploiting the strengths of both techniques. In the particle MCMC samplers, SMC algorithms provide high-dimensional proposal distributions, and Markov transitions guide the SMC ensemble to sufficiently explore the target distribution. The transition is realized by a conditional SMC technique, in which a reference trajectory from the previous step is kept throughout the current step of SMC sampling.

In this study, we sample the posterior by PGAS [28], a particle MCMC method that enhances the mixing of the Markov chain by sampling the ancestor of the reference trajectory. For the SMC, we use an optimal particle filter, which takes advantage of the linear Gaussian observation model and the Gaussian transition density of the state variables in our current SEBM. More generally, when the observation model is nonlinear and the transition density is non Gaussian, the optimal particle filter can be replaced by implicit particle filters [12, 35] or local particle filters [38, 39]; we refer to [8, 25, 49] for other data assimilation techniques. The details of the algorithm are provided in Section 7.3.

3 Ill-posedness and regularized posteriors

In this section, we first demonstrate and then analyze the failure of standard Bayesian inference of the parameters with the posteriors in (2.10). The standard Bayesian inference of the parameters fails in the sense that the posterior (2.10) tends to have a large probability mass at non-physical parameter values. In the process of approximating the posterior by samples, the values of these samples often either hit the (upper or lower) bounds in Table 1 when we use a uniform prior or exceed these bounds when we use a Gaussian prior. As we shall show next, the standard Bayesian inverse problem is numerically ill-posed because the Fisher information matrix is ill-conditioned, which makes the inference numerically unreliable. Following the idea of regularization in variational approaches, we propose to use regularized posteriors in the Bayesian inference. This approach unifies the Bayesian and the variational approaches: the MAP is the minimizer of the regularized cost function in the variational approach, but the Bayesian approach quantifies the uncertainty of the estimator by the posterior.

3.1 Model settings and tests

Based on the physical upper and lower bounds in Table 1, we consider two priors for the parameters: a uniform distribution on these intervals and a Gaussian distribution centered at the median and with three standard deviations in the interval, as listed in Table 2.

Throughout this study, we shall consider a relatively small numerical mesh for the SPDE with only 12 nodes for the finite elements. Such a small mesh provides a toy model that can neatly represent the spatial structure on the sphere, while allowing for systematic assessments of statistical properties of the Bayesian inference with moderate computational costs. Numerical tests show that the above FEM semi-backward Euler scheme is stable for a time step size $\Delta t=0.01$ and a stochastic forcing with scale $\sigma_{f}=0.1$ (see Section 7.1 for more details about the discretization). A typical realization of the solution is shown in Figure 1 (left and middle), where we present the solution on the sphere at a fixed time with the 12-node finite element mesh, as well as the trajectories of all 12 nodes.

The standard deviation of the observation noise is set to $\sigma_{\epsilon}=0.01$, i.e. one order of magnitude smaller than the stochastic forcing and two orders of magnitude smaller than the climatological mean.

We first assume that six out of the 12 nodes are observed; we discuss results obtained using sparser or denser observations in the discussion section. Figure 1 also shows the climatological probability histogram of the true state variables and the partial noisy observations. The climatological distribution of the observations is close to that of the true state variables (with a slightly larger variance due to the noise). The histograms show that the state variables are centered around $1$ and vary mostly in the interval $[0.92,1.05]$. We shall use a Gaussian approximation based on the climatological distribution of the partial noisy observations as a prior to constrain the state variables.

We summarize the settings of numerical tests in Table 3.

3.2 Ill-posedness of the standard Bayesian inference of parameters

By the Bernstein-von Mises theorem (see e.g. [48, Chaper 10]), the posterior distribution of the parameters conditional on the true state data approaches the likelihood distribution as the data size increases. That is, $p(\theta|u_{1:N})$ in (2.14) becomes close to the likelihood distribution $p(u_{1:N}|\theta)$ (which can be viewed as a distribution of $\theta$) as the data size increases. Therefore, if the likelihood distribution is numerically degenerate (in the sense that some components are undetermined), then the Bayesian posterior will also become close to degenerate, so that the Bayesian inference for parameter estimation will be ill-posed. In the following, we show that for this model the likelihood is degenerate even if the full states are observed with zero observation noise and that the maximum likelihood estimators have large nonphysical fluctuations (particularly when the states are noisy). As a consequence, the standard Bayesian parameter inference fails by yielding nonphysical samples.

We show first that the likelihood distribution is numerically degenerate because the Fisher information matrix is ill-conditioned. Following the transition density (2.7), the log-likelihood of the state $\{u_{1:N}\}$ is

[TABLE]

where $c$ is a constant independent of $(\theta,u_{1:N})$. Since $\mu_{\theta}(\cdot)$ is linear in $\theta$ (cf. Equation (7.19)), the likelihood function is quadratic in $\theta$ and the corresponding scaled Fisher information matrix is

[TABLE]

where the vectors $G_{\theta,k}(u_{n})\in\mathbb{R}^{d_{b}}$ are defined in (7.20). Figure 2 shows the means and standard deviations of the condition numbers (the ratio between the maximum and the minimum singular values) of the Fisher information matrices from 100 independent simulations. Each of these simulations generates a long trajectory of length $10^{5}$ using a parameter drawn randomly from the prior, and computes the Fisher information matrices using the true trajectory of all 12 nodes, for subsamples of lengths $N$ ranging from $10^{2}$ to $10^{5}$. For both Gaussian and uniform priors, the condition numbers are on the scale of $10^{8}-10^{11}$ and therefore the Fisher information matrix is ill-conditioned. In particular, the condition number increases as the data size increased, due to the ill-posedness of the inverse problem of parameter estimation.

The ill-conditioned Fisher information matrix leads to highly variable maximum likelihood estimators (MLE), computed from $\mathbf{F}_{N}\theta=b_{N}$ with $b_{N}=\frac{1}{N}\left(\sum_{n=1}^{N}G_{\theta,k}(u_{n})^{T}\mathbf{R}^{-1}(u_{n+1}-\mathbf{M}_{\Delta t}^{-1}\mathbf{M}_{0}u_{n})\right)_{k=0,1,4}$, which follows from (7.20).

The ill-posedness is particularly problematic when $\{u_{1:N}\}$ is observed with noise, as the ill-conditioned Fisher information matrix amplifies the noise in observations and leads to nonphysical estimators. Figure 3 shows the means and standard deviations of errors of MLEs computed from true and noisy trajectories in 100 independent simulations. In each of these simulations, the “noisy” trajectory is obtained by adding a white noise with standard deviation $\sigma_{\epsilon}=0.01$ to a “true” trajectory generated from the system with a true parameter randomly drawn from the prior. For both Gaussian and uniform priors, the standard deviations and means of the errors of the MLE from the noisy trajectories are one order of magnitude larger than those from true trajectories. In particular, the variations are large when the data size is small. For example, when $N=100$, the standard deviation of the MLE for $\theta_{0}$ from noisy observations is on the order of $10^{3}$, two orders of magnitude larger than its physical range in Table 2. The standard deviations decrease as the data size increases, at the expected rate of $1/\sqrt{N}$. However, the errors are too large to be practically reduced by increasing the size of data: for example, a data size $N=10^{10}$ is needed to reduce the standard deviation of $\theta_{4}$ to less than $0.1$ (which is about 10% the size of the physical range $[-6.00,-4.80]$ as specified in Table 2). In summary, the ill-posedness leads to parameter estimators with large variations that are far outside the physical ranges of the parameters.

3.3 Regularized posteriors

To overcome the ill-posedness of the parameter estimation problem, we introduce strongly regularized posteriors by normalizing the likelihood function. In addition, to prevent unphysical values of the states, we further regularize the state variables in the likelihood by an uninformative climatological prior. That is, consider the regularized posterior:

[TABLE]

where $Z:=\int p(\theta)\left[\frac{p^{c}(u_{1:N})p_{\theta}(u_{1:N})p_{\theta}(y_{1:N}|u_{1:N})}{p_{\theta}(y_{1:N})}\right]^{1/N}d\theta du_{1:N}$ is a normalizing constant and $p^{c}(u_{1:N})$ is the prior of the states estimated from a Gaussian fit to climatological statistics of the observations, neglecting correlations. That is, we set $p^{c}(u_{1:N})$ as

[TABLE]

with $\sigma_{c}=2\sqrt{\sigma_{o}^{2}-\sigma_{\epsilon}^{2}}$, where $u_{c}$ and $\sigma_{o}$ are the mean and standard deviation of the observations over all states. Here the multiplicative factor 2 aims for a larger band to avoid an overly narrow prior for the states.

This prior can be viewed as a joint distribution of the state variables assuming all components are independent identically Gaussian distributed with mean $u_{c}$ and variance $\sigma_{c}^{2}$. Clearly, it uses the minimum amount of information about the state variables, and we expect it can be improved by taking into consideration spatial correlations or additional field knowledge in practice.

The regularized posterior can be viewed as an extension of the regularized cost function in the variational approach. In fact, the negative logarithm of the regularized posterior is the same (up to a multiplicative factor $\frac{1}{N}$ and an additive constant $\log Z-\frac{1}{N}\log p_{\theta}(y_{1:N})$) as the cost function in variational approaches with regularization. More precisely, we have

[TABLE]

where $C_{y_{{}_{1:N}}}(\theta,u_{1:N})$ is the cost function with regularization:

[TABLE]

When the prior is Gaussian, the regularization corresponds to Tikhonov regularization. Therefore, the regularized posterior extends the regularized cost function to a probability distribution, with the maximum of the posterior (MAP) being the minimizer of the regularized cost function.

The regularized posterior normalizes the likelihood by an exponent $1/N$. This normalization allows for a larger weight (more trust) on the prior, which can then sufficiently regularize the singularity in the likelihood and therefore reduces the probability of nonphysical samples. Intuitively, it avoids the shrinking of the likelihood as the data size increases. When the system is ergodic, the sum \frac{1}{N}\sum_{n=1}^{N}\log\big{[}p_{\theta}(u_{n}|u_{n-1})p(y_{n}|u_{n})\big{]} converges to the spatial average \mathbb{E}[\log\big{[}p_{\theta}(U_{n}|U_{n-1})p(y_{n}|U_{n})] with respect to the invariant measure as $N$ increases. While being effective, this factor may not be optimal [37] and we leave the exploration of optimal regularization factors to future work.

In the sampling of the regularized posterior, we update the state variable $U_{1:N}$ conditionally on $\theta$ and $y_{1:N}$ by sampling $p^{c}(u_{1:N})p_{\theta}(u_{1:N}|\theta,y_{1:N})$ (with $p_{\theta}(u_{1:N}|\theta,y_{1:N})$ specified in (2.13)) using SMC methods. Compared to the standard PMCMC algorithm outlined in Section 2.4, the only difference occurs when we update the parameter $\theta$ conditional on the estimated states $u_{1:N}$. Instead of (2.14), we draw a sample of $\theta$ from the regularized posterior

[TABLE]

4 Bayesian inference with regularized posteriors

The regularized posteriors are approximated by the empirical distribution of samples drawn using particle MCMC methods, specifically particle Gibbs with ancestor sampling (PGAS, see Section 7.3) in combination with SMC using optimal importance sampling (see Section 7.2). In the following section, we first diagnose the Markov chain and choose a reasonable chain length for subsequent analyses. We then present the results of parameter estimation and state estimation.

In all the tests presented in this study, we use only $M=5$ particles for the SMC, as we can be confident of the Markov chain produced by the particle MCMC methods converging to the target distribution based on theoretical results [2, 28]. In general, the more particles are used, the better the SMC algorithm (and hence the particle MCMC methods) will perform, at the price of increased computational cost.

4.1 Diagnosis of the Markov Chain

To ensure that the Markov Chain generated by PGAS is well-mixed and to find a length for the chain such that the posterior is acceptably approximated, we shall assess the Markov chain by three criteria: the update rate of states; the correlation length of the Markov chain; and the convergence of the marginal posteriors of the parameters. These empirical criteria are convenient and, as we discuss below, have found to be effective in our study. We refer to [13] for a detailed review of various criteria for diagnosing MCMC.

The update rate of states is computed at each time of the state trajectory $u_{1:N}$ along the Markov chain. That is, at each time, we say the state is updated from the previous step of of the Markov chain if any entry of the state vector changes. The update rate measures the mixing of the Markov chain. In general, an update rate above 0.5 is preferred, but a high rate close to 1 is not necessarily the best. Figure 4 shows the update rates of typical simulations for both the Gaussian prior and the uniform prior. For both priors, the update rates are above 0.5, indicating a fast mixing of the chain. The rates tend to increase with time (except for the first time step) to a value close to 1 at the end of the trajectory. This phenomenon agrees with the particle depletion nature of the SMC filter: when tracing back in time to sample the ancestors, there are fewer particles and therefore the update rate is lower. The high update rate at time $t=1$ step is due to our initialization of the particles near the equilibrium, which increases the possibility of ancestor updates in PGAS. We also note that the uniform prior has update rates close to 1 at all times, much higher than the rates of the Gaussian prior. Higher update rates occur for the uniform prior because the deviations of parameter samples from the previous values are larger, resulting in an increased probability of updating the reference trajectory in the conditional SMC.

We test the correlation length of the Markov chain by finding the smallest lag at which the empirical autocorrelation functions (ACF) of the states and the parameters are close to zero. Figure 5 shows the empirical ACFs of the parameters and states at representative nodes, computed using a Markov chain with length 10000. The ACFs approach zero within a time lag of around 40 ( based on a threshold value of 0.1) for the Gaussian prior, and within a time lag of around 5 for the uniform prior. The smaller correlation length in the uniform prior case is again due to the larger parameter variation in the uniform prior case than the Gaussian prior case.

The relatively small decorrelation length of the Markov chain indicates that we can accurately approximate the posterior by a chain of a relatively short length. This result is demonstrated in Figure 6, where we plot the empirical marginal posteriors of the parameters, using Markov chains of three different lengths: $L=1000,5000,10000$. The marginal posteriors with $L=1000$ are reasonably close to those with $L=10^{4}$, and those with $L=5000$ are almost identical to those with $L=10^{4}$. In particular the marginal posteriors with $L=10^{3}$ capture the shape and spread of the distributions for $L=10^{4}$. Therefore, a Markov chain with length $L=10^{4}$ provides a reasonably accurate approximation of the posterior. Hence, we use Markov chains with length $L=10^{4}$ in all simulations from here on. This choice of chain length may be longer than necessary, but allows for confidence that the results are robust.

In summary, based on the above diagnosis of the Markov chain generated by PMCMC, to run many simulations for statistical analysis of the algorithm within a limited computation cost, we use chains with length $L=10^{4}$ to approximate the posterior. For the SMC algorithm, we use only five particles. The number of observations in time is $N=100$.

4.2 Parameter estimation

One of the main goals in Bayesian inference is to quantify the uncertainty in the parameter-state estimation by the posterior. We access the parameter estimation by examining the samples of the posterior in a typical simulation, for which we consider the scatter plots and marginal distributions, the maximum of the posterior (MAP) and the posterior mean. We also examine the statistics of the MAP and the posterior mean in 100 independent simulations. In each simulation, the parameters are drawn from the prior distribution of $\theta$. Then, a realization of the SEBM is simulated. Finally, observations are created by applying the observation model to the SEBM realization.

The empirical marginal posteriors of the parameters $\theta=(\theta_{0},\theta_{1},\theta_{4})$ in two typical simulations, for the Gaussian and the uniform priors, are shown in Figure 7. The top row presents scatter plots of samples along with the true values of the parameters (asterisks), and the bottom row presents the marginal posteriors for each parameter in comparison with the priors.

In the case of the Gaussian prior, the scatter plots show a posterior that is far from Gaussian, with clear nonlinear dependence between $\theta_{0}$ and the other parameters. The marginal posteriors of $\theta_{0}$ and $\theta_{1}$ are close to their priors, with larger tails (to the left for $\theta_{0}$ and to the right for $\theta_{1}$). The marginal distribution of $\theta_{4}$ concentrates near the center of the prior with larger tail to the right. The posterior has the most probability mass near the true values of $\theta_{0}$ and $\theta_{1}$, which are in the high probability region of the prior. However, it has no probability mass near the true value of $\theta_{4}$ – which is of a low probability in the prior.

In the case of the uniform prior, the scatter plots show a concentration of probability near the boundaries of the physical range. The marginal posteriors of $\theta_{0}$ and $\theta_{1}$ clearly deviate from the priors, concentrating near the parameter bounds (the upper bound for $\theta_{0}$ and the lower bound for $\theta_{1}$ in this realization); the marginal posterior of $\theta_{4}$ is close to the prior with slightly more probability mass for large values.

Further tests show that the posterior is not sensitive to changes in the true values of the parameters. This fact is demonstrated in Figure 8, which presents the marginal distributions for another set of true values of the parameters (but without changing the priors). Though the data change when the true parameters change, the posteriors, in comparison with those Figure 7, change little for both cases of Gaussian and uniform prior.

The non-Gaussianity of the posterior (including the concentration near the boundaries), its insensitivity to changes in the true parameter, and its limited reduction of uncertainty from the prior (Figures 7 - 8) are due to the degeneracy of the likelihood distribution and to the strong regularization. Recall that the degenerate likelihood leads to MLEs with large variations and biases, with the standard deviation of the estimators of $\theta_{0}$ and $\theta_{1}$ being about 10 times larger than those of $\theta_{4}$ (see Figure 3). As a result, when regularized by the Gaussian prior, the components $\theta_{0}$ and $\theta_{1}$, which are more under-determined by the likelihood, are constrained mainly by the Gaussian prior and therefore their marginal posteriors are close to their marginal priors. In contrast, the component $\theta_{4}$ is forced to concentrate around the center of the prior but with a large tail. While dramatically reducing the large uncertainty of $\theta_{0}$ and $\theta_{1}$ in the ill-conditioned likelihood, the regularized posterior still exhibits a slightly larger uncertainty than the prior for these two components.

In the case of the uniform prior, it is particularly noteworthy that the marginal posteriors of $\theta_{0}$ and $\theta_{1}$ differ more from their priors than the parameter $\theta_{4}$. These results are the opposite of what was found for the Gaussian prior. Such differences are due to the different mechanism of “regularization” by the two priors. The Gaussian prior eliminates the ill-posedness by regularizing the ill-conditioned Fisher information matrix with the covariance of the prior. So, the information in the likelihood, e.g. the bias and the correlations between $(\theta_{0},\theta_{1})$ and $\theta_{4}$, are preserved in the regularized posterior. The uniform prior, on the other hand, cuts the support of the degenerate likelihood and rejects out-of-range samples. As a result, the correlation between $\theta_{0}$ and $\theta_{1}$ is preserved in the regularized posterior because they feature similar variations, but the correlations between $(\theta_{0},\theta_{1})$ and $\theta_{4}$ are weakened (Figure 7).

In practice, one is often interested in a point estimate of parameters. Commonly used point estimators are the MAP and the posterior mean. Figures 7-8 show that both the MAP and the posterior mean can be far away from the truth for Gaussian as well as uniform priors. In particular, in the case of the uniform prior, the MAP values are further away from the truth than the posterior mean. In the case of the Gaussian prior, the MAP values do not present a clear advantage or disadvantage over the posterior mean.

Table 5(a) shows the means and standard deviations of the errors of the posterior mean and MAP from 100 independent simulations. In each simulation and for each prior, we drew a parameter sample from the prior and generated a trajectory of observations, and then estimated jointly the parameters and states. The table shows that both posterior mean and MAP estimates are generally biased, consistent with the biases in Figures 7 and 8. More specifically, in the case of the Gaussian prior, the MAP has slightly smaller biases than the posterior mean, but the two have almost the same variances. Both are negatively biased for $\theta_{0}$ and slightly positively biased for $\theta_{1}$ and $\theta_{4}$. In the case of the uniform prior, the MAP features biases and standard deviations which are about 50% larger than those of the posterior mean. Both estimators exhibit large positive biases in $\theta_{0}$, large negative biases in $\theta_{1}$, and small positive biases in $\theta_{4}$.

4.3 State estimates

The state estimation aims both to filter out the noise from the observed nodes and to estimate the states of unobserved nodes. We access the state estimation by examining the ensemble of the posterior trajectories in a typical simulation, for which we consider the marginal distributions and the coverage probability of 90% credible intervals. We also examine the statistics of these quantities in 100 independent simulations.

We present the ensemble of posterior trajectories at an observed node in Figure 9 and at an unobserved node in Figure 10. In each of these figures, we present the ensemble mean with a one-standard-deviation band, in comparison with the true trajectories, superimposed on the ensembles of all sample trajectories at these nodes. We also present histograms of samples at three instants of time: $t=20$, $t=60$ and $t=100$.

Figure 9 shows that the trajectory of the observed node is well estimated by the ensemble mean, with a relative error of 0.7%. Recall that the observation noise leads to a relative error of about 1%, so the posterior filters out 30% of the noise. Also note that the ensemble quantifies the uncertainty of the estimation, with the true trajectory being mostly enclosed within a one-standard-deviation band around the ensemble mean. Further, the histograms of samples at the three time instants show that the ensemble generally concentrates near the truth. In the Gaussian prior case, the peak of the histogram decreases as time increases. partially due to the degeneracy of SMC when we trace back the particles in time. In the uniform prior case, the ensembles are less concentrated than those in the Gaussian case, due to the wide spread of the parameter samples (Figure 7).

Figure 10 shows sample trajectories of an unobserved node. Despite the fact that the node is unobserved, the posterior means have relative errors of 0.8% and 3.3% in cases of Gaussian and uniform priors respectively, with a one-standard-deviation band covering the true trajectory at most times. While the sparse observations do cause large uncertainties for both priors, the histograms of samples show that the ensembles concentrate near the truth. Particularly, in the case of Gaussian prior, the peaks of the histogram are close to the true states, even when the histograms form a multi-modal distribution due to the degeneracy of SMC.

We find that the posterior is able to filter out the noise in the observed nodes and reduce the uncertainty in the unobserved nodes from the climatological distribution. In particular, in the case of the Gaussian prior, the ensemble of posterior samples concentrates near the true state at both observed and unobserved nodes and substantially reduces the uncertainty. In the case of the uniform prior, the ensemble of posterior samples spreads more widely, and only slightly reduces the uncertainty.

The coverage probability (CP), the proportion of the states whose 90% credible intervals contain the true values, is 95% in the Gaussian prior case and 92% for the uniform prior in the above simulation. The target probability is 90% as in this case 90% of the true values would be covered by 90% credible intervals. The values indicate statistically meaningful uncertainty estimates, for example larger uncertainty ranges at nodes with higher mean errors. The slight over-dispersiveness, i.e. higher CPs than the target probabilities, might be a result of the large uncertainty in the parameter estimates.

Table 6 shows the means and standard deviations of the relative errors and CPs in state estimation by the posterior mean in 100 independent simulations, averaging over observed and unobserved notes. The relative errors at each time $t$ are computed by averaging the error of the ensemble mean (relative to the true value) over all the nodes. The relative error of the trajectory is the average over all times along the trajectory. The relative errors are 1.14% and 2.39% respectively for the cases of Gaussian and uniform prior. These numbers are a result of averaging over the observed and unobserved nodes. Note that the relative errors are similar at different times $t=(20,60,100)$, indicating that the MCMC is able to ameliorate the degeneracy of the SMC to faithfully sample the posterior of the states.

In the Gaussian prior case, the CPs are above the target probability in the 100 independent simulations with a mean of 96%. This supports the finding from above that the posteriors are slightly over-dispersive due to the large uncertainty in the parameter estimates. The standard deviation is very small with 2% which indicates the robustness of the Gaussian prior model. In the uniform prior case, the CPs are much lower with a mean of 73%. This might be a result of larger biases compared to the Gaussian prior case which are not compensated by larger uncertainty estimates. In addition, the standard deviation is much higher in the uniform prior case with 31%. This shows that this case is less robust than the Gaussian prior case.

5 Discussion

5.1 Observing fewer nodes

We tested the consequences of having sparser observations in space, e.g. observing only two out of the 12 nodes. In the Gaussian prior case, in a typical simulation with the same true parameters and observation data as in Section 4.2, the relative error in state estimation increases slightly, from 0.7% to 0.8% for the observed node and from 0.8% to 1.1% for the unobserved node. As a result, the overall error increases. The parameter estimates show small but noticeable changes (see Figure 11): the posteriors of the parameters have slightly wider support and the posterior means and MAPs exhibit slightly larger errors than those in Section 4.2.

We also ran 100 independent simulations to investigate sampling variability in the state and parameter estimates. Table 6(b) reports the means and standard deviations of the relative errors of the posterior mean trajectory, and CPs for state estimation in these simulations. The Gaussian prior case shows small increases in both the means and the standard deviations of errors, as well as slightly lower and less robust CPs. This confirms the results quoted above for a typical simulation. The uniform prior case shows almost negligible error and CP increases. Table 5(b) reports the mean and standard deviations of the posterior means and MAP for parameter estimation in these simulations. Small changes in comparison to the results in Table 5(a) are found. These small changes are due to the strong regularization that has been introduced to overcome the degeneracy of the likelihood.

5.2 Observing a longer trajectory.

When the length $N$ of the trajectory of observation increases, the exponent of the regularized posterior (3.5), viewed as a function of $\theta$ only, tends to its expectation with respect to the ergodic measure of the system, i.e. $\frac{1}{N}C_{y_{{}_{1:N}}}(\theta,u_{1:N})\xrightarrow[]{N\to\infty}\mathbb{E}[C_{y_{{}_{1:N}}}(\theta,u_{1:N})]$ almost surely. As a result, the marginal posterior tends to be stable as $N$ increases. This result indicates that an increase of data size has a limited effect on the regularized posterior of parameters. This fact is verified by numerical tests with $N=1000$, in which the marginal posteriors only have a slightly wider support than those in Figure 7 with $N=100$.

In general, the number of observations needed for the posterior to reach a steady state depends on the dimension of the parameters and the speed of convergence to the ergodic measure of the system. Here we have only three parameters and the SEBM converges to its stationary measure exponentially (in fewer than 10 time steps), therefore $N=100$ is large enough to make the posterior be close to the steady state.

When the trajectory is long, a major issue is the computational cost from sampling the posterior of the states. Note that as $N$ increases, the dimension of the states in the posterior increases, demanding a longer Markov chain to explore the target distribution. In numerical tests with $N=1000$, the correlation length of the Markov chain is at least 100, about four times the correlation length found for $N=100$. Therefore, to obtain the same number of effective samples as before, we would need a Markov chain with length at least four times the previous length, say, $L=4\times 10^{4}$. The computational cost increases linearly in $NL$, with each step requiring an integration of the SPDE. The high computational cost, an instance of the well-known “curse of dimensionality”, renders the direct sampling of the posterior unfeasible. Two groups of methods could reduce the computational cost and make the Bayesian inference feasible. The first group of methods, dynamical model reduction, exploits the low-dimensional structure of the stochastic process to develop low-dimensional dynamical models which efficiently reproduce the statistical-dynamical properties needed in the SMC (see e.g. [10, 32, 9, 24] and the references therein). The other group of methods approximates the marginal posterior of the parameter by reduced order models for the response of the data to parameters (see e.g. [33, 6, 14, 11, 31, 22]). In a paleoclimate reconstruction context, the number of observations will generally be determined by available observations and the length of the reconstruction period rather than by computational considerations. We leave these further development of efficient sampling methods for long trajectories as a direction of future research.

5.3 Estimates of the nonlinear function

One goal of parameter estimation is to identify the nonlinear function $g_{\theta}$ (specified in (2.2)) in the SEBM. The posterior of the parameters also quantifies the uncertainty in the identification of $g_{\theta}$. Figure 12 shows the nonlinear function $g_{\theta}$ associated with the true parameters and with the MAPs and posterior means presented in Figure 7, superposed on an ensemble of the nonlinear function evaluated with all the samples. Note that in the Gaussian prior case, the true and estimated functions $g_{\theta}$ are close even though $\theta_{4}$ is estimated with large biases by either the posterior mean or by the MAP. In the uniform prior case, the posterior mean has a smaller error than the MAP and leads to a better estimate of the nonlinear function. In either case, the large band of the ensemble represents a large uncertainty in the estimates.

For the Gaussian prior, neither the posterior distribution of the equilibrium state $u_{e}$ (for which $g_{\theta}(u_{e})=0$) nor of the feedback strength $dg_{\theta}/du(u_{e})$ are substantially changed from the corresponding priors. Both experience only a small reduction of uncertainty. In contrast, the posterior distributions are narrower than the priors for the uniform prior case - although the posterior means and MAPs are both biased.

6 Conclusions and future work

We have investigated the joint state-parameter estimation of a nonlinear stochastic energy balance model (SEBM) motivated by the problem of spatial-temporal paleoclimate reconstruction from sparse and noisy data, for which parameter estimation is an ill-posed inverse problem. We introduced strongly regularized posteriors to overcome the ill-posedness by restricting the parameters and states to physical ranges and by normalizing the likelihood function. We considered both a uniform prior and a more informative Gaussian prior based on the physical ranges of the parameters. We sampled the regularized high-dimensional posteriors by a Particle Gibbs with Ancestor Sampling (PGAS) sampler that combines Markov Chain Monte Carlo (MCMC) with an optimal particle filter to exploit the forward structure of the SEBM.

Results show that the regularization overcomes the ill-posedness in parameter estimation and leads to physical posteriors quantifying the uncertainty in parameter-state estimation. Due to the ill-posedness, the posterior of the parameters features a relatively large uncertainty. This result implies that there can be a large uncertainty in point estimators such as the posterior mean or the maximum a posteriori (MAP), the latter of which corresponds to the minimizer in a variational approach with regularization. Despite the large uncertainty in parameter estimation, the marginal posteriors of the states generally concentrate near the truth, reducing the uncertainty in state reconstruction. In particular, the more informative Gaussian prior leads to much better estimations than the uniform prior: the uncertainty in the posterior is smaller, the MAP and posterior mean have smaller errors in both state and parameter estimates, and the coverage probabilities are higher and more robust.

Results also show that the regularized posterior is robust to spatial sparsity of observations, with sparser observations leading to slightly larger uncertainties due to less information. However, due to the need of regularization to overcome ill-posedness, the uncertainty in the posterior of the parameters cannot be eliminated by increasing the number of observations in time. Therefore, we suggest alternative approaches, such as re-parametrization of the nonlinear function according to the climatological distribution or nonparametric Bayesian inference (see e.g. [36, 18]) to avoid ill-posedness.

The ill-posedness of the parameter estimation problem for the model we have considered is of particular interest because the form of the nonlinear function $g_{\theta}(u)$ is not arbitrary but is motivated by the physics of the energy budget of the atmosphere. The fact that wide ranges of the parameters $\theta_{i}$ are consistent with the “obserations” even in this highly idealized setting indicates that surface temperature observations themselves may not be sufficient to constrain physically-important parameters such as albedo, graybody thermal emissivity, or air-sea exchange coefficients separately. While state-space modeling approaches allow reconstruction of past surface climate states, it may be the case that the associated climate forcing may not contain sufficient information to extract the relative contributions of the individual physical processes that produced it.

Acknowledgments.

This research started in a working group supported by the Statistical and Applied Mathematical Sciences Institute (SAMSI). FL thanks Prof. Peter Jan van Leeuwen, Prof. Kayo Ide, Prof. Mauro Maggioni, Prof. Xuemin Tu, and Dr. Wenjun Ma for helpful discussions. FL is supported by the National Science Foundation under Grant DMS-1821211. NW thanks Andreas Hense and Douglas Nychka for inspiring discussions. NW was supported by the German Federal Ministry of Education and Research (BMBF) through the Palmod project (FKZ: 01LP1509D). NW thanks the German Research Foundation (code RE3994-2/1) for funding. AM acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC), and thanks SAMSI for hosting him in the autumn of 2017.

7 Appendix: Technical details of the estimation procedure

7.1 Discretization of the SEBM

Finite element representation in space We discretize the SEBM in space by finite element methods (see e.g. [1])see e.g.. Denote by $\{\phi_{i}(\xi)\}_{i=1}^{d_{b}}$ the finite element basis functions, and approximate the solution $u(t,\xi)$ by

[TABLE]

The coefficients $\widehat{u}_{i}$ are determined by the following weak Galerkin projection of the SEBM (2.1)

[TABLE]

where $\phi$ is a continuously differentiable compactly supported test function and the integral $\int_{0}^{t}\langle f(s,\cdot),\phi\rangle$ is an Itô integral.

For convenience, we write this Galerkin approximate system in vector notation. Denote

[TABLE]

Taking $\phi=\phi_{j}$, $j=1,\dots,d_{b}$ in equation (7.2) and using the symmetry of the inner product, we obtain a stochastic integral equation for the coefficient $U(t)\in\mathbb{R}^{d_{b}}$:

[TABLE]

To simplify notation, we denote the mass and stiffness matrices by

[TABLE]

which are symmetric, tri-diagonal, positive definite matrices in $\mathbb{R}^{d_{b}\times d_{b}}$, and we denote the nonlinear term as

[TABLE]

The above stochastic integral equation can then be written as

[TABLE]

The mesh on the sphere and the matrices $\mathbf{M}_{0}$ and $\mathbf{M_{1}}$ are computed with the R package INLA [27, 5].

Representation of the nonlinear term. The parametric nonlinear functional $G_{\theta}(U(t))$ is approximated using the finite elements. We approximate each spatial integration over an element-triangle in $\langle g_{\theta}(U_{n}^{T}\Phi),\Phi\rangle$ by the volume of the triangular pyramid whose height is the value of the nonlinear function at the center of the element-triangle $\mathcal{T}_{k}$, i.e.

[TABLE]

where $\xi_{k}^{c}$ is the center of the triangle $\mathcal{T}_{k}$. In the discretized system, we assume that this approximation has a negligible error and take it as our nonlinear functional. In vector notation, it reads

[TABLE]

where $A_{\mathcal{T}}=\left(\frac{\mathrm{Area}(\mathcal{T}_{k})}{3}\right)\in\mathbb{R}^{d_{b}\times d_{e}}$ with $d_{e}$ denoting the number of triangle elements and the matrix $\mathbf{A}=\left(\phi_{i}(\xi_{k}^{c})\right)\in\mathbb{R}^{d_{e}\times d_{b}}$, such that the function $g_{\theta}(\mathbf{A}U(t))$ is interpreted as element-wise evaluation. For the nonlinear function $g_{\theta}$ in (2.2), we can write the above nonlinear term as

[TABLE]

where $\circ k$ denotes entry-wise product of the array.

Representation of the stochastic forcing. Following [26], the stochastic forcing $f(t,\xi)$ is approximated by its linear finite element truncation,

[TABLE]

with the stochastic processes $\{f_{i}(t),i=1,\dots,d_{b}\}$ being spatially correlated and white in time. Note that for $\nu=0.1$ and $\rho>0$ in the Matérn covariance (2.4), the process $f(t,\xi)$ is the stationary solution of the stochastic Laplace equation

[TABLE]

where $W$ is a spatio-temporal white noise [52, 53]. Computationally efficient approximations of the forcing process are obtained using the GMRF approximation of [26] which generates $F(t)\equiv\left(f_{1}(t),f_{2}(t),\dots,f_{d_{b}}(t)\right)$ by solving (7.14). That is, using the above finite element notation, we solve for each time $t$ the linear system

[TABLE]

where the random vector $\langle\Phi,W(t,\cdot)\rangle:=\left(\langle\phi_{1},W(t,\cdot)\rangle,\dots,\langle\phi_{d_{b}},W(t,\cdot)\rangle\right)$ is Gaussian with mean 0 and covariance $\mathbf{M}_{0}$. Solving (7.15) yields

[TABLE]

where $\mathbf{M}_{\rho}:=(\rho^{-2}\mathbf{M}_{0}+\mathbf{M}_{1})$.

Semi-backward Euler time integration. Equation (7.9) is integrated in time by a semi-backward Euler scheme

[TABLE]

where $U_{n}$ is the approximation of $U(t_{n})$ with $t_{n}=n\Delta t$, and $\{F_{n}\}$ is a sequence of iid random vectors with distribution $\mathcal{N}\left(0,\sigma_{f}^{2}\mathbf{M}_{\rho}^{-1}\mathbf{M}_{0}\mathbf{M}_{\rho}^{-1}\right)$, with the matrix $\mathbf{M}_{\Delta t}$ denoting

[TABLE]

Efficient generation of the Gaussian field. It follows from (7.15) that $\mathbf{M}_{0}F_{n}$ is Gaussian with mean zero and covariance $\mathbf{M}_{0}\mathbf{M}_{\rho}^{-1}\mathbf{M}_{0}\mathbf{M}_{\rho}^{-1}\mathbf{M}_{0}$. Note that while $\mathbf{M}_{\rho}$ is a sparse matrix, its inverse matrix $\mathbf{M}_{\rho}^{-1}$ is not. To efficiently use the sparseness of $\mathbf{M}_{\rho}$, following [26], we approximate $\mathbf{M}_{0}$ by $\widehat{\mathbf{M}}_{0}:=\text{diag}(\langle\phi_{i},1\rangle)$ and compute the noise $\mathbf{M}_{0}F_{n}$ by $\mathbf{C}^{-1}\mathcal{N}(0,I_{d})$, where $\mathbf{C}$ is the Cholesky factorization of the inverse of the covariance matrix (called precision matrix) $\widehat{\mathbf{M}}_{0}^{-1}\mathbf{M}_{\kappa}\widehat{\mathbf{M}}_{0}^{-1}\mathbf{M}_{\kappa}\widehat{\mathbf{M}}_{0}^{-1}$.The precision matrix is a sparse representation of the inverse of the covariance. Therefore, the matrix $\mathbf{C}$ is also sparse and the noise sequence can be efficiently generated.

In summary, we can write the discretized SEBM in the form

[TABLE]

where the deterministic function $\mu_{\theta}(\cdot)$ is given by

[TABLE]

with $G_{\theta,k}(U_{n}):=\Delta t\mathbf{M}_{\Delta t}^{-1}A_{\mathcal{T}}(\mathbf{A}U(t))^{\circ k}$, and $\{W_{n}\}$ is a sequence of iid Gaussian noise with mean 0 and covariance $\mathbf{R}$:

[TABLE]

7.2 SMC with optimal importance sampling

SMC methods approximate the target density $p_{\theta}(u_{1:N}|y_{1:N})$ sequentially by weighted random samples called particles (hereafter we drop the subindex $\theta$ to simplify notation)

[TABLE]

with $\sum_{m=1}^{M}w_{n}^{m}=1$. These weighted samples are drawn sequentially by importance sampling based on the recurrent formation

[TABLE]

More precisely, suppose that at time $n$, we have weighted samples $\{U_{1:n-1}^{m},w_{n-1}^{m}\}_{m=1}^{M}$. One first draws a sample $U_{n}^{m}$ from an easy to sample importance density $q(u_{n}|y_{n},U_{n-1}^{m})$ that approximates the “incremental density” which is proportional to $p\left(y_{n}|u_{n}\right)p(u_{n}|U^{m}_{n-1})$ for each $m=1,\dots,M$, and computes incremental weights

[TABLE]

which account for the discrepancy between the two densities. One then assigns normalized weights $\{w_{n}^{m}\,\propto\,w_{n-1}^{m}\alpha_{n}^{m}\}_{m=1}^{M}$ to the concatenated sample trajectories $\{U_{1:n}^{m}\}_{m=1}^{M}$.

A clear drawback of the above procedure is that all but one of the weights $\{w_{n}^{m}\}$ will become close to zero as the number of iterations increases, due to the multiplication and normalization operations. To avoid this, one replaces the unevenly weighted samples $\{(U_{n-1}^{m},w_{n-1}^{m})\}$ by uniformly weighted samples from the approximate density $\widehat{p}_{\theta}(u_{n-1}|y_{1:N-1})$. This is the well-known resampling technique. In summary, the above operations are carried out as follows:

• (i)

draw random indices $\{A_{n-1}^{m}\}_{m=1}^{M}$ according to the discrete probability distribution $\mathbb{F}(\cdot|w^{1:M}_{n-1})$ on the set $\{1,\dots,M\}$, which is defined as

[TABLE]

• (ii)

for each $m$, draw a sample $U_{n}^{m}$ from $q(u_{n}|y_{n},U_{n-1}^{A_{n-1}^{m}})$ and set $U_{1:n}^{m}:=(U_{n-1}^{A_{n-1}^{m}},U_{n}^{m})$;

• (iii)

compute and normalize the weights

[TABLE]

The above SMC sampling procedure is called sequential importance sampling with resampling (SIR) (see e.g.[15]) and is summarized in Algorithm 1.

**Optimal importance sampling. ** Note that the conditional transition density of the states $p_{\theta}(u_{n+1}|u_{n})$ in (2.7) is Gaussian and the observation model in (2.8) is linear and Gaussian. These facts allow for a Gaussian optimal importance density $q(u_{n}|y_{n},U_{n-1}^{m})$ that is proportional to $p\left(y_{n}|u_{n}\right)p(u_{n}|U^{m}_{n-1})$ for each $m=1,\dots,M$:

[TABLE]

with the mean $\mu_{n}^{m}$ and the covariance $\mathbf{\Sigma}$ given by

[TABLE]

Drawbacks of SMC. While the resampling technique prevents $w_{n}^{m}$ from being degenerate at each current time $n$, SMC algorithms suffer from the degeneracy (or particle depletion) problem: the marginal distribution $\widehat{p}(u_{n}|(y_{1:N}))$ becomes concentrated on a single particle as $N-n$ increases because each resampling step reduces the number of distinct particles of $u_{n}$. As a result, the estimate of the joint density $p(u_{1:N}|y_{1:N})$ of the trajectory deteriorates as time $N$ increases.

7.3 Particle Gibbs and PGAS

The framework of particle MCMC introduced in [2] is a systematic combination of SMC and MCMC methods, exploiting the strengths of both techniques. Among the various particle MCMC methods, we focus on the particle Gibbs sampler (PG) that uses a novel conditional SMC update [2], as well as its variant, the particle Gibbs with ancestor sampling (PGAS) sampler [28], because they are best fit for sampling our joint parameter and state posterior.

The PG and PGAS samplers use a conditional SMC update step to realize the transition between two steps of the Markov chain while ensuring that the target distribution will be the stationary distribution of the Markov chain. The basic procedure of a PG sampler is as follows:

• •

Initialization: draw $\theta(1)$ from the prior distribution $p(\theta)$. Run an SMC algorithm to generate weighted samples $\{U_{1:N}^{m},w_{N}^{m}\}_{m=1}^{M}$ for $p_{\theta(1)}(u_{1:N}|y_{1:N})$ and draw $U_{1:N}(1)$ from these weighted samples.

• •

Markov chain iteration: for $l=1,\cdots,L-1$,

Sample $\theta(l+1)$ from the marginal posterior $p(\theta|y_{1:N},U_{1:N}(l))$ given by (2.14).

• 2.

Run a conditional SMC algorithm, conditioned on $U_{1:N}(l)$, which is called the reference trajectory. That is, in the SMC algorithm, the $M$-th particle is required to move along the reference trajectory by setting $U_{n}^{M}=U_{n}(l)$. Draw other samples from the importance density, and normalize the weights and resample all the particles as usual. This leads to weighted samples $\{U_{1:N}^{m},w_{N}^{m}\}_{m=1}^{M}$ with $U_{1:N}^{M}=U_{1:N}(l)$.

• 3.

Draw $U_{1:N}(l+1)$ from the above weighted samples.

• •

Return the Markov chain $\{\theta(l),U_{1:N}(l)\}_{l=1}^{L}$.

The conditional SMC algorithm is the core of PG samplers. It retains the reference path throughout the resampling steps by deterministically setting $U_{1:N}^{M}=U_{1:N}(l)$ and $A_{n}^{M}=M$ for all $n$, while sampling the remaining $M-1$ particles according to a standard SMC algorithm. The reference path interacts with the other paths by contributing a weight $w_{n}^{M}$. This is the key to ensuring that the PG Markov chain converges to the target distribution. A potential risk of the PG sampler is that it yields a poorly mixed Markov chain, because the reference trajectory tends to dominate the SMC ensemble trajectories.

The PGAS sampler increases the mixing of the chain by connecting the reference path to the history of other particles by assigning an ancestor to the reference particle at each time. This is accomplished by drawing a sample for the ancestor index $A^{M}_{n-1}$ of the reference particle, which is referred to as ancestor sampling. The distribution of the index $A^{M}_{n-1}$ is determined by the likelihood of connecting $U_{n}(l)$ to the particles $\{U_{n-1}^{m}\}_{m=1}^{M}$, in other words, according to weights

[TABLE]

The above weight $\widetilde{\alpha}_{n-1|n}^{m}$ can be seen as a posterior probability, where the importance weight $w_{n-1}^{m}$ is the prior probability of the particle $U_{n-1}^{m}$, and the product $p_{{}_{\theta(l+1)}}(U_{n}(l)|U^{m}_{n-1})p(y_{n}|U_{n}(l))$ is the likelihood that $U_{n}(l)$ originates from $U_{n-1}^{m}$ conditional on observation $y_{n}$. In short, the PGAS sampler assigns the reference particle $U_{n}(l)$ an ancestor $A^{M}_{n-1}$ that is drawn from the distribution $\mathbb{F}(A^{M}_{n-1}=k|\widetilde{w}_{n-1|n}^{1:M})=\widetilde{w}_{n-1|n}^{k}.$

The above conditional SMC with ancestor sampling within PGAS is summarized in Algorithm 2.