# Estimation of Spatially Varying Parameters with Application to Hyperbolic SPDES

Authors: David Angwenyi
Journal: Journal of Applied Mathematics
 (2023)

---

## Abstract

Parameter estimation is a growing area of interest in statistical signal processing. Some parameters in real-life applications vary in space as opposed to those that are static. Most common methods in estimating parameters involve solving an optimization problem where the cost function is assembled variously, for example, maximum likelihood and maximum a posteriori methods. However, these methods do not have exact solutions to most real-life problems. It is for this reason that Monte Carlo methods are preferred. In this paper, we treat the estimation of parameters which vary with space. We use the Metropolis-Hastings algorithm as a selection criterion for the maximum filter likelihood. Comparisons are made with the use of joint estimation of both the spatially varying parameters and the state. We illustrate the procedures employed in this paper by means of two hyperbolic SPDEs: the advection and the wave equation. The Metropolis-Hastings procedure registers better estimates.


---

## Body

## 1. Introduction

Most state-space models are characterized by, among other things, parameters—which can be constant or varying. A parameter is comprehended and signified as a measurable factor which defines a model and influences its operation. As the parameter changes, so does the model; expressed differently, a parameter is unique to a model it characterizes. The choice of a certain model, therefore, is achieved by choosing the right parameters. It occurs more often than not—in hidden Markov models, for instance—that measurements are available but the underlying signal is not readily apparent. This forms an example where parameter estimation is paramount: measurements are used to learn the model parameters, which, in turn, are used to fit the model.Consider a signal and measurement equations, with a spatially varying parameter signified by ad-dimensional vector, θ,
(1a)Signal:dxt=fxt,θdt+gxtdβt;t0≤t,(1b)Measurement:dyt=hxt,θdt+R1/2tdηt;t0≤t,where the terms are as described in Table 1.Table 1
Terms, their descriptions, and dimensions.
TermNameDimensionxtState vectorn×1fxt,tDrift functionn×1gxt,tDiffusion functionn×mβt,t>t0Brownian motion processm×1ytOutput vectorr×1hxt,tSensor functionr×1RtTime-function matrixr×rηt,t>t0Standard Brownian motion processr×1Parameter estimation problem concerns finding the optimal parameter so that the signal best fits the data [1, 2]. This is, classically, achieved by an optimization procedure where a cost function is minimized [3]. The cost function mostly defines the discrepancy between the state and the measurements. Intuitively, parameter estimation can be seen as a procedure for seeking a parameter value that gives the least discrepancy between the state and the corresponding measurements (also known as the algebraic distance or the residual). The method of least squares has been extensively used to define an objective function. Given the increment in measurement, dyt, of the state, xt, at time t, the objective function, Jθ, in the least-squares sense, is given by
(2)Jθ=∫t0twtdyt−hxt,θdt2,where wt is the weighting function.Most commonly used procedures in the framework of least squares include the following: linear least squares, orthogonal least squares, gradient-weighted least squares, and bias-corrected renormalization. This paper, however, attends not to the study of least-squares approaches, they being outside the scope of its design. Suffice it to only direct the interested reader to the article [4] for an elaborate explanation and application of least-squares methods in computer vision. Instead, we study parameter estimation by means of filtering. But before that, we mention a few merits and demerits of least squares and other methods defined by algebraic distances.The use of algebraic distances in defining a cost function is computationally efficient, and closed-form solutions are possible. The end result, however, is not satisfactory. This is due, in one part, to the fact that the objective function is mostly not invariant with respect to Euclidean transformations, for example, translation. This limits the coordinate systems to be used. In the other part, outliers may not contribute to the parameters the same way as inliers [4]. Other more satisfactory parameter estimation methodologies are highly desired. We consider, in this paper, the use of Bayesian inference techniques.Estimation of parameters by means of a filter can be achieved in a number of ways, one of them being the use of the filter evidence, or its near approximation, and selection criteria for parameters which give a reasonable estimate of the evidence. The second method involves updating the parameters and the state at the same time. This is known as dual estimation, which further subdivides into joint estimation and a dual filter. Joint estimation entails subjoining the vector of parameters to the state vector to form an extended state-space. The filter is then implemented and run forward in time with the hope of filter convergence to the optimal state and parameter values. A dual filter, on the other hand, involves implementing a filter for the state and parameters simultaneously. The filter provides a self-correcting mechanism which may lead to convergence of state and parameter estimates.The rest of this paper is arranged as follows. In the first place, we introduce some notions on Bayesian inference of parameters where we pass the limit in the mean in order to move from a discrete setting to a continuum in time. We introduce the different ways in which a filter can be used for parameter estimation. In anticipation of application in the later part of the paper, we introduce the stochastic advection and wave equation together with their spatial discretisation. We then illustrate how parameters are to be estimated by means of a filter likelihood and the dual filters. Results and discussions follow, and the conclusion forms the closing part of this paper.


## 2. Bayesian Parameter Inference

In Bayesian inference of parameters, the parameters are treated as a random variable. The parameter is assigned a prior,πt0θ, based on some initial belief. Let tn such that tn+1>tn∀n=0,1,2,⋯N be a partition of the interval t0,T and let δt=tn+1−tn. Bayes’ rule gives the joint posterior of parameters and the states;
(3a)πt0,Tx,θYT≈l.i.m.δt⟶0N⟶∞πt0:tNxt0:tN,θyt0:tN=l.i.m.δt⟶0N⟶∞πt0:tNyt0:tNxt0:tN,θπt0:tNxt0:tNθπt0θπt0:tNyt0:tN,where YT=yt0,T,
(3b)πt0:tNyt0:tNxt0:tN,θ=∏n=1Nπtnytnxtn,θ,πt0:tNxt0:tNθ=πt0xt0θ∏n=1Nπtnxtnxtn−1,θ.Now to arrive at the marginal posterior of parameters, we integrate out the states from the joint posterior of states and parameters, Equation (3a):
(4)πt0:tNθyt0:tN=∫πt0:tNyt0:tNxt0:tN,θπt0:tNxt0:tNθπt0θπt0:tNyt0:tNdxt0:tN.It turns out that direct computation of the integral in (4) is difficult, especially with the increase in measurements [2]. This challenge is circumvented through the use of recursive techniques which include the use of filters and smoothers, maximum a posteriori (MAP) estimates, and drawing samples from the posterior using Markov Chain Monte Carlo (MCMC) methods.To use recursive methods, we begin with the following parameter posterior(5a)πθYT≈l.i.m.δt⟶0N⟶∞πθyt0:tN∝l.i.m.δt⟶0N⟶∞πt0:tNyt0:tNθπt0θ,where
(5b)πt0:tNyt0:tNθ=∏n=1Nπtnytnyt1:tn−1,θ=∏n=1Nπtnytnxtn,θπtnxtnyt1:tn−1,θdxtn.The state’s predictive distribution,πtnxtnyt1:tn−1,θ, can be computed recursively as follows [2]:
(6a)πtnxtnyt1:tn−1,θ=∫πtnxtnxtn−1,θπtn−1xtn−1yt1:tn−1,θdxtn−1.Instead of working with the posterior,πθYT, it is quite convenient to use the negative log-posterior obtained by expressing the posterior as follows:
(7)πθYT≈l.i.m.δt⟶0N⟶∞πθyt0:tN∝l.i.m.δt⟶0N⟶∞exp−ψTθ,where the energy function, ψTθ, is given by
(8)ψTθ=−logπt0:tNyt0:tNθ−logπt0θ.The maximum a posteriori (MAP) estimate can then be obtained by the mode of the posterior distribution or, equivalently, the minimum of the energy function, the latter of which is easier to compute; that is,(9)θ^MAP=argmaxθπθYT=argminθψTθ.One demerit of the MAP estimate is that it yields a point estimate of the parameter posterior and therefore ignores the dispersion of the estimate. Setting the prior,πt0θ, to be a uniform density, (9) yields a maximum likelihood estimate.


## 3. Metropolis-Hastings Method

Metropolis-Hastings (named after Nicholas Constantine Metropolis (1915-1999) and Wilfred Keith Hastings (1930-2016)) [5] is a Markov Chain Monte Carlo sampling algorithm with optimal convergence. It is premised on detailed balance and ergodicity. Given a probability density, say πθ, from which it is difficult to sample (for instance, if the said distribution is known to a normalization constant), and given another distribution ρθ, say from which it is easy to sample, then detailed balance is the condition
(10)πθkρθkθk+1=πθk+1ρθk+1θk,where ρθkθk+1 is a transition distribution. The detailed balance condition is necessary for any random walk to asymptotically converge to a stationary distribution. Ergodicity is meant that there is a nonzero probability in moving from a state to any other state in a Markov Chain.The following algorithm summarises the Metropolis-Hastings procedure.Algorithm 1: Metropolis-Hastings.
1: Drawθ~ from ρθ~θk2: Setθk+1⟵θ~ with probability α=min1,πθ~ρθkθ~/πθkρθ~θk3: Otherwise, setθk+1⟵θk


## 4. Dual Estimation

Dual estimation comprehends simultaneous estimation of state and parameters by means of an appropriate filter. The self-correcting mechanism of the filter is taken advantage of to converge to both the true state and the true parameters. Depending on the initial parameter, the filter sooner or later converges to the true parameter value. Dual estimation can be achieved in two ways: joint estimation and by a dual filter [6–8].


### 4.1. Joint Estimation (Augmented State-Space)

In joint estimation, the state vector is augmented with the vector of parameters to form an extended state-space, and then, the filter is run forward in time for an update of both the state and the parameters. The parameters are induced with artificial dynamics or are made to assume a random walk; that is, respectively,(11)dzt=ζt;t0≤t,where
(12a)dzt=dxtdθt,ζt=fxt,θdt+gxtdβt0,or where
(12b)ζt=fxt,θdt+gxtdβtσdχt,in which χt,t>t0 is a d-dimensional standard Brownian motion vector process and σ is a small constant. A filter is then implemented with the augmented state zt in the place of xt. The demerit of this method is that the extended state-space has an increased degree of freedom owing to many unknowns, of both the state and the parameters, which renders the filter unstable and intractable, especially in nonlinear models [8].


### 4.2. Dual Filter

Dual filtering of the state and parameters is attained by use of two filters, one for state update and another for updating parameters, both run simultaneously. The two filters interact symbiotically in that one provides the update of the state to be used by the other, whilst the other provides an update of the parameters to be used by the former. A very good example in literature is the dual extended Kalman filter [9] used for estimating neural network models and the weights. In this case, the state is the model signal and the weights are parameters. Another example appears in [10] where a dual filter comprising the ensemble transform particle filter (ETPF) and the feedback particle filter (FPF) is used for simultaneous estimation of the state of a wave equation and its speed.The model for the dual filter of our consideration comprises ad-dimensional vector equation of artificial dynamics of parameters together with the state-space model, Equations (1a) and (1b):
(13a)Parameter:dθt=0,t0≤t,(13b)Signal:dxt=fxt,θtdt+gxtdβt,t0≤t,(13c)Measurement:dyt=hxtdt+R1/2tdηt,t0≤t,where the nomenclature and dimensions remain as stipulated for Equations (1a) and (1b).In the following, we introduce the equations to which we shall apply dual filters in estimating spatially varying parameters.


## 4.1. Joint Estimation (Augmented State-Space)

In joint estimation, the state vector is augmented with the vector of parameters to form an extended state-space, and then, the filter is run forward in time for an update of both the state and the parameters. The parameters are induced with artificial dynamics or are made to assume a random walk; that is, respectively,(11)dzt=ζt;t0≤t,where
(12a)dzt=dxtdθt,ζt=fxt,θdt+gxtdβt0,or where
(12b)ζt=fxt,θdt+gxtdβtσdχt,in which χt,t>t0 is a d-dimensional standard Brownian motion vector process and σ is a small constant. A filter is then implemented with the augmented state zt in the place of xt. The demerit of this method is that the extended state-space has an increased degree of freedom owing to many unknowns, of both the state and the parameters, which renders the filter unstable and intractable, especially in nonlinear models [8].


## 4.2. Dual Filter

Dual filtering of the state and parameters is attained by use of two filters, one for state update and another for updating parameters, both run simultaneously. The two filters interact symbiotically in that one provides the update of the state to be used by the other, whilst the other provides an update of the parameters to be used by the former. A very good example in literature is the dual extended Kalman filter [9] used for estimating neural network models and the weights. In this case, the state is the model signal and the weights are parameters. Another example appears in [10] where a dual filter comprising the ensemble transform particle filter (ETPF) and the feedback particle filter (FPF) is used for simultaneous estimation of the state of a wave equation and its speed.The model for the dual filter of our consideration comprises ad-dimensional vector equation of artificial dynamics of parameters together with the state-space model, Equations (1a) and (1b):
(13a)Parameter:dθt=0,t0≤t,(13b)Signal:dxt=fxt,θtdt+gxtdβt,t0≤t,(13c)Measurement:dyt=hxtdt+R1/2tdηt,t0≤t,where the nomenclature and dimensions remain as stipulated for Equations (1a) and (1b).In the following, we introduce the equations to which we shall apply dual filters in estimating spatially varying parameters.


## 5. Application Equations

### 5.1. Advection Equation

We take up an advection equation, excited with the space-time white noise process, with some diffusion term added to it, and on a periodic spatial domain of lengthL, which we write as follows:
(14)∂u∂t=∂Cxu∂x+μ∂2u∂x2+σβ̇x,t,0≤t≤Tt×0≤x≤L,where Cx is a spatially varying velocity (of which constant velocity is a special case), σ is a constant, and ux,t is the function to be obtained, whose function describes the state of the signal. μ is a constant whilst β̇x,t is the space-time white noise process where the dot denotes the singularity of the noise process.Equation (14) needs a remark owing to the roughness of the stochastic-forcing term β̇x,t, which is a mixed distributional derivative of the Brownian sheet. As is well known (see [11] for details), the Brownian sheet is nowhere differentiable. We, however, use the method introduced in [12]; that is, we approximate the noise term as follows. Let the domain 0≤t≤Tt×0≤x≤L be tessellated into rectangles tn,tn+1×xi,xi+1 of dimensions δt×δx each, for n=1,2,3,⋯,T and i=1,2,3,⋯,N so that δt=Tt/T and δx=L/N. Then,
(15)β̇x,t≈1δxδt∑i=1N∑n=1Tωi,nδxδt1/2χixχnt,where ωi,ni=1N are independent and identically distributed random variables of mean 0 and unit variance. χix and χnt are characteristic functions and are given by
(16)χnt=1,iftn≤t≤tn+1,0,otherwise,χix=1,ifxi≤x≤xi+1,0,otherwise.By a three-point upwind scheme in space [13], we discretise (14) and arrive at the following:
(17)duidt≈3Ciui−4Ci−1ui−1+Ci−2ui−22δx+μ3ui+2−16ui+1+26ui−16ui−1+3ui−24δx2+σ1δxω̇i,where δx is the spatial step size and ωi,t,t>t0 is the standard Brownian motion process. The ith grid point is represented by xi=iδx. With this notation, ui,n is understood to mean the value of u at the ith grid point at time tn.Time discretisation, by means of the Euler-Maruyama scheme, leads to(18)ui,tn+1=ui,tn+3Ciui,tn−4Ci−1ui−1,tn+Ci−2ui−2,tn2δxδt+μ3ui+2,tn−16ui+1,tn+26ui,tn−16ui−1,tn+3ui−2,tn4δx2δt+σδt1/2δxωi,tn,where ωi,tn is a random variable of mean 0 and variance 1. The time increment, δt>0, is such that the limit of δui as δt⟶0 is dui. Furthermore, n=1,2,3,⋯,T. We use the following initial value.
(19)ux,t0=sinx.Moreover,Cx=eλx where
(20)λx=sin2πx.Considering every grid point in (18) leads to a vector representation of the signal u. To do so requires the following shorthand for operations:
(21)D1ui≔3ui−4ui−1+ui−22δx,∀i=1,2,3,⋯,N,(22)D1Tui≔ui+2−4ui+1+3ui2δx,∀i=1,2,3,⋯,N,so that
(23)D1D1Tui≔3ui+2−12ui+1+9ui4δx2+−4ui+1+16ui−12ui−14δx2+ui−4ui−1+3ui−24δx2=3ui+2−16ui+1+26ui−16ui−1+3ui−24δx2,∀i=1,2,3,⋯,N.We finally have(24)utn+1=utn+Ftnutnδt+Gtnωtn,where utn is an N-dimensional column vector at time tn comprising of elements ui,tn,i=1,2,3,⋯,N, and
(25)Ftn=D1Cdiag−μD1D1T,Gtn=σδt1/2δxIN×N,in which Cdiag is a diagonal matrix made of the elements of C and IN×N is the Nth-order identity matrix.The surface and contour plots for the stochastic advection equation are shown below, that is, whenσ=0.1. The ruggedness evident in Figure 1 is consequent upon the addition of noise to the underlying dynamics.Figure 1
Contour and surface plots for a single realisation of the stochastic advection equation under the following setting:L=2π, N=100, δx=L/N, δt=0.007, T=1000, σ=0.1, μ=0.01, Cx=eλx where λx=sin2πx, and ut0=sinx.
(a)(b)In the next subsection, we introduce the wave equation.


### 5.2. Wave Equation

The wave equation—for our consideration—is given by(26)∂2u∂t2=∂Cx∂u/∂x∂x+μ∂3u∂x2∂t+σβ̇x,t,0≤t≤Tt×0≤x≤L,where Cx=eλx is the wave velocity and is for a wave travelling in a heterogeneous medium and ux,t is the function to be obtained, whose function, as in the previous subsection, describes the state of the signal. Moreover, σ is a constant and β̇x,t, as before, is the space-time white noise.We employ mixed difference schemes to discretise (26) in space, so that we have
(27a)duidt≈pi,(27b)dpidt≈Ci+1wi+1−Ciwiδx+μpi+1−2pi+pi−1δx2+1δxω̇i,where wi≔ui−ui−1/δx and δx is the spatial step size. For time integration, we use Verlet’s method, because of its geometric properties, namely, volume preservation, symplecticity, conservation of first integrals, and reversibility [14, 15]—whose method, applied to the deterministic part of Equations (27a) and (27b), yields
(28a)ui,tn+1=ui,tn+pi,tn+1/2δt,(28b)pi,tn+1/2=pi,tn+Ci+1wi+1,tn−Ciwi,tnδxδt2+μpi+1,tn−2pi,tn+pi−1,tnδx2δt2,(28c)pi,tn+1=pi,tn+1/2+Ci+1wi+1,tn+1−Ciwi,tn+1δxδt2+μpi+1,tn+1−2pi,tn+1+pi−1,tn+1δx2δt2+σδt1/2δxωi,tn,where δt is the time step and wi,tn≔ui,tn−ui−1,tn/δx. Substituting Equation (28b) into Equations (28a) and (28c), we get
(29a)ui,tn+1=ui,tn+pi,tn+Ci+1wi+1,tn−Ciwi,tnδxδt2+μpi+1,tn−2pi,tn+pi−1,tnδx2δt2,(29b)pi,tn+1=pi,tn+Ci+1wi+1,tn−Ciwi,tnδxδt2+μpi+1,tn−2pi,tn+pi−1,tnδx2δt2+Ci+1wi+1,tn+1−Ciwi,tn+1δxδt2+μpi+1,tn+1−2pi,tn+1+pi−1,tn+1δx2δt2+σδt1/2δxωi,tn.We use the following initial values:(30a)ux,0=exp−4x−0.5L.2,(30b)px,0=0,where L is the length of the domain. Now, Cx=eλx where
(31)λx=sinx.Considering every grid point leads to a vector representation of the signalu. We use the following shorthand:
(32)D2ui≔ui−ui−1δx,∀i=1,2,3,⋯,N.Equations (29a) and (29b) then become
(33)u¯tn+1=u¯tn+Ftnu¯tn+Gtnω¯tn,where
(34)Ftn=−I2N×2N+IN×N−δt22D2CdiagxD2TδtIN×N−δt22μD2D2T−δt2D2CdiagxD2TIN×N−δt2μD2D2T×IN×N0N×Nδt2D2CdiagxD2TIN×N−δt2μD2D2T−1,u¯=pu,Gtn=0N×N0N×N0N×Nσδt1/2δxIN×N,ω¯=ωω,where I2N×2N is the identity matrix of order 2N whilst 0N×N is an Nth-order null matrix. 0N is an Nth-dimensional null vector. Ftn∈ℝ2N×2N, Gtn∈ℝ2N×2N, and ω¯∈ℝN.The contour and surface plots for a single realisation of stochastic wave equation are in Figure2.Figure 2
Plots for a single realisation of the stochastic wave equation. The settings areL=2π, N=100, δx=L/N, δt=0.005, and T=1000, μ=0.01, σ=0.2, Cx=eλx where λx=sinx, and ut0=exp−4x−0.5L.2.
(a)(b)For the purpose of estimating the speed of the wave, which varies spatially, we give a little prelude to estimation of spatially varying parameters in the next section.


## 5.1. Advection Equation

We take up an advection equation, excited with the space-time white noise process, with some diffusion term added to it, and on a periodic spatial domain of lengthL, which we write as follows:
(14)∂u∂t=∂Cxu∂x+μ∂2u∂x2+σβ̇x,t,0≤t≤Tt×0≤x≤L,where Cx is a spatially varying velocity (of which constant velocity is a special case), σ is a constant, and ux,t is the function to be obtained, whose function describes the state of the signal. μ is a constant whilst β̇x,t is the space-time white noise process where the dot denotes the singularity of the noise process.Equation (14) needs a remark owing to the roughness of the stochastic-forcing term β̇x,t, which is a mixed distributional derivative of the Brownian sheet. As is well known (see [11] for details), the Brownian sheet is nowhere differentiable. We, however, use the method introduced in [12]; that is, we approximate the noise term as follows. Let the domain 0≤t≤Tt×0≤x≤L be tessellated into rectangles tn,tn+1×xi,xi+1 of dimensions δt×δx each, for n=1,2,3,⋯,T and i=1,2,3,⋯,N so that δt=Tt/T and δx=L/N. Then,
(15)β̇x,t≈1δxδt∑i=1N∑n=1Tωi,nδxδt1/2χixχnt,where ωi,ni=1N are independent and identically distributed random variables of mean 0 and unit variance. χix and χnt are characteristic functions and are given by
(16)χnt=1,iftn≤t≤tn+1,0,otherwise,χix=1,ifxi≤x≤xi+1,0,otherwise.By a three-point upwind scheme in space [13], we discretise (14) and arrive at the following:
(17)duidt≈3Ciui−4Ci−1ui−1+Ci−2ui−22δx+μ3ui+2−16ui+1+26ui−16ui−1+3ui−24δx2+σ1δxω̇i,where δx is the spatial step size and ωi,t,t>t0 is the standard Brownian motion process. The ith grid point is represented by xi=iδx. With this notation, ui,n is understood to mean the value of u at the ith grid point at time tn.Time discretisation, by means of the Euler-Maruyama scheme, leads to(18)ui,tn+1=ui,tn+3Ciui,tn−4Ci−1ui−1,tn+Ci−2ui−2,tn2δxδt+μ3ui+2,tn−16ui+1,tn+26ui,tn−16ui−1,tn+3ui−2,tn4δx2δt+σδt1/2δxωi,tn,where ωi,tn is a random variable of mean 0 and variance 1. The time increment, δt>0, is such that the limit of δui as δt⟶0 is dui. Furthermore, n=1,2,3,⋯,T. We use the following initial value.
(19)ux,t0=sinx.Moreover,Cx=eλx where
(20)λx=sin2πx.Considering every grid point in (18) leads to a vector representation of the signal u. To do so requires the following shorthand for operations:
(21)D1ui≔3ui−4ui−1+ui−22δx,∀i=1,2,3,⋯,N,(22)D1Tui≔ui+2−4ui+1+3ui2δx,∀i=1,2,3,⋯,N,so that
(23)D1D1Tui≔3ui+2−12ui+1+9ui4δx2+−4ui+1+16ui−12ui−14δx2+ui−4ui−1+3ui−24δx2=3ui+2−16ui+1+26ui−16ui−1+3ui−24δx2,∀i=1,2,3,⋯,N.We finally have(24)utn+1=utn+Ftnutnδt+Gtnωtn,where utn is an N-dimensional column vector at time tn comprising of elements ui,tn,i=1,2,3,⋯,N, and
(25)Ftn=D1Cdiag−μD1D1T,Gtn=σδt1/2δxIN×N,in which Cdiag is a diagonal matrix made of the elements of C and IN×N is the Nth-order identity matrix.The surface and contour plots for the stochastic advection equation are shown below, that is, whenσ=0.1. The ruggedness evident in Figure 1 is consequent upon the addition of noise to the underlying dynamics.Figure 1
Contour and surface plots for a single realisation of the stochastic advection equation under the following setting:L=2π, N=100, δx=L/N, δt=0.007, T=1000, σ=0.1, μ=0.01, Cx=eλx where λx=sin2πx, and ut0=sinx.
(a)(b)In the next subsection, we introduce the wave equation.


## 5.2. Wave Equation

The wave equation—for our consideration—is given by(26)∂2u∂t2=∂Cx∂u/∂x∂x+μ∂3u∂x2∂t+σβ̇x,t,0≤t≤Tt×0≤x≤L,where Cx=eλx is the wave velocity and is for a wave travelling in a heterogeneous medium and ux,t is the function to be obtained, whose function, as in the previous subsection, describes the state of the signal. Moreover, σ is a constant and β̇x,t, as before, is the space-time white noise.We employ mixed difference schemes to discretise (26) in space, so that we have
(27a)duidt≈pi,(27b)dpidt≈Ci+1wi+1−Ciwiδx+μpi+1−2pi+pi−1δx2+1δxω̇i,where wi≔ui−ui−1/δx and δx is the spatial step size. For time integration, we use Verlet’s method, because of its geometric properties, namely, volume preservation, symplecticity, conservation of first integrals, and reversibility [14, 15]—whose method, applied to the deterministic part of Equations (27a) and (27b), yields
(28a)ui,tn+1=ui,tn+pi,tn+1/2δt,(28b)pi,tn+1/2=pi,tn+Ci+1wi+1,tn−Ciwi,tnδxδt2+μpi+1,tn−2pi,tn+pi−1,tnδx2δt2,(28c)pi,tn+1=pi,tn+1/2+Ci+1wi+1,tn+1−Ciwi,tn+1δxδt2+μpi+1,tn+1−2pi,tn+1+pi−1,tn+1δx2δt2+σδt1/2δxωi,tn,where δt is the time step and wi,tn≔ui,tn−ui−1,tn/δx. Substituting Equation (28b) into Equations (28a) and (28c), we get
(29a)ui,tn+1=ui,tn+pi,tn+Ci+1wi+1,tn−Ciwi,tnδxδt2+μpi+1,tn−2pi,tn+pi−1,tnδx2δt2,(29b)pi,tn+1=pi,tn+Ci+1wi+1,tn−Ciwi,tnδxδt2+μpi+1,tn−2pi,tn+pi−1,tnδx2δt2+Ci+1wi+1,tn+1−Ciwi,tn+1δxδt2+μpi+1,tn+1−2pi,tn+1+pi−1,tn+1δx2δt2+σδt1/2δxωi,tn.We use the following initial values:(30a)ux,0=exp−4x−0.5L.2,(30b)px,0=0,where L is the length of the domain. Now, Cx=eλx where
(31)λx=sinx.Considering every grid point leads to a vector representation of the signalu. We use the following shorthand:
(32)D2ui≔ui−ui−1δx,∀i=1,2,3,⋯,N.Equations (29a) and (29b) then become
(33)u¯tn+1=u¯tn+Ftnu¯tn+Gtnω¯tn,where
(34)Ftn=−I2N×2N+IN×N−δt22D2CdiagxD2TδtIN×N−δt22μD2D2T−δt2D2CdiagxD2TIN×N−δt2μD2D2T×IN×N0N×Nδt2D2CdiagxD2TIN×N−δt2μD2D2T−1,u¯=pu,Gtn=0N×N0N×N0N×Nσδt1/2δxIN×N,ω¯=ωω,where I2N×2N is the identity matrix of order 2N whilst 0N×N is an Nth-order null matrix. 0N is an Nth-dimensional null vector. Ftn∈ℝ2N×2N, Gtn∈ℝ2N×2N, and ω¯∈ℝN.The contour and surface plots for a single realisation of stochastic wave equation are in Figure2.Figure 2
Plots for a single realisation of the stochastic wave equation. The settings areL=2π, N=100, δx=L/N, δt=0.005, and T=1000, μ=0.01, σ=0.2, Cx=eλx where λx=sinx, and ut0=exp−4x−0.5L.2.
(a)(b)For the purpose of estimating the speed of the wave, which varies spatially, we give a little prelude to estimation of spatially varying parameters in the next section.


## 6. Estimating Spatially Varying Parameters

Varying parameters exist naturally, an example being in the speed of a wave moving in a heterogeneous media. Varying parameters present a challenge owing to their varying nature as opposed to static (in time) parameters. There are three broad categories of varying parameters: spatially varying parameters, parameters that vary with time, and parameters that vary with both time and space. Estimation of time-varying parameters, albeit for deterministic models, and application to estimation of parameters in a HIV/AIDS model, appears in [16]—in which least-squares methods have been used. Although consistent estimates are realised, extension to nonlinear models remains to be realised. In this chapter, we study spatially varying parameters using filtering algorithms. We also provide applications in order to evaluate the performance of the algorithms introduced here.The strategy employed in this study comprises of three steps:(1)
Express the parameterθx as a Fourier series with a given number of modes, say, Nm and collect all the coefficients in a vector λ(2)
By means of an appropriate filtering algorithm, estimate the vector of hyperparameters,λ(3)
Substitute the estimated constant coefficients back in the Fourier series to obtain an estimate of the parameterθxTo illustrate this method, we employ it to estimate the velocity of a wave travelling in a heterogeneous media. Let the parameter,Cx, x∈ℝN, where N is the number of dimensions, be given by
(35)Cx=expfx,where
(36)fx=a0+∑k=1Nmakk2sinkx+∑k=1Nmbkk2coskx,k=1,2,3,⋯,Nm,where the coefficients a0, ak, and bk are drawn randomly from a normal distribution of mean 0 and variance 1. This way, the parameter Cx will be positive for all values of x. The aim of this section is to obtain an estimate Cxest of the wave velocity Cx by means of the Kalman-Bucy filter (KBF) and ensemble Kalman-Bucy filter (EnKBF) [17]. To this end, we use a function
(37)gx=A0+∑k=1ℵmAkk2sinkx+∑k=1ℵmBkk2coskx,k=1,2,3,⋯,ℵm,where ℵm≤Nm so that
(38)Cxest=expgx.Letλ∈ℝ2ℵm+1 be a vector whose elements are the coefficients of the function gx. That is,
(39)λ=A0A1B1A2B2⋮AℵmBℵm.Estimation of the spatially varying parameter,Cx, is equivalent to estimating the parameters that form the vector λ. We consider two methods: using the likelihood with the Metropolis-Hastings method and using a dual filter. Let us now consider each method in turn.


## 7. Using the Likelihood with Metropolis-Hastings Method

We now adapt Algorithm1 to estimating the vector of parameters, λ. We let each parameter λi have artificial dynamics with the transition density
(40)ζk+1λi,k+1λi,k=Nλi,kcosϕ,ωℵsinϕ,where ℵ is the mode number and the constants ϕ and ω are to be chosen. The density π is defined by the filter likelihood,
(41)πδyt0,Tu~t0,T,λk=l.i.m.δt⟶0N⟶∞πt0:tNδyt0:tNu~t0:tN,λk,where u~tn is the filter prediction of the state at time tn. Notice that the parameters, λ, are implicitly contained in the dynamics.Example 1.
Advection equation.
We take the advection equation, of Section5.1, and the given initial conditions. Furthermore, let there be time-continuous measurements of the state, u, given by
(42)dyt=Htux,tdt+Qtdηx,t,where ηt,x is the standard space-time Brownian motion process. The initial values of ut, βt, and ηt,x are uncorrelated. The aim is to estimate the spatially varying velocity, Cx, by means of filter likelihood and the Metropolis-Hastings algorithm.We follow the discretisation described in Section5.1 for the advection equation.The measurements’ equation, (42), is discretised in time using the Euler-Maruyama scheme to yield
(43)δytn=Htnux,tnδt+R1/2tnδηtn,upon substituting R=QQT/δx and where EδηtnδηtnT=Ir×rδt. The observation likelihood pdf for the KBF is Gaussian since the initial condition and the observation errors are Gaussian. So is the posterior pdf. With observation increments expressed as in (43), the observation increment likelihood pdf is
(44)πδyt0:tNu~tN,λk∝∏n=0Nexp−12δytn−Htnu~tnδtδtRtn2.Similarly, the filter forecast pdf is(45)πutNδytN−1:t0∝∏n=0Nexp−12utn−u~tnPtn2.The next step is to implement a KBF and EnKBF and to obtain the likelihood at each time step. This is followed by implementing a Metropolis-Hastings algorithm.Algorithm 2: KBF likelihood with MH.
Require:δt, ℵm, N, ut0, πtn, and λk.Ensure:λkk=1T.1:fork=1toNdo2: Drawλ~~Nλi,kcosϕ,ω/ℵmsinϕ∀i=1,2,⋯,2ℵm+1.3: ComputeCkx=expgx,λk.4:forn=1toT, δt>0do5:  Run a single-step KBF prediction meanu~tn+1=utn+Ftn,λkutnδt6:  Run a single-step KBF prediction covarianceP~tn+1=Ptn+Ftn,λkPtnδt+PtnFTtn,λkδt+GtnGTtnδt7:  Run a single-step KBF analysis meanutn+1=u~tn+1+PtnHTtnR−1tndytn−Htn,λku~tnδt8:  Run a single-step KBF analysis covariancePtn+1=P~tn+1+GtnGTtnδt−P~tnHTtnR−1tnHtnPtnδt9:end for10:Metropolis-Hastings11: Computeαratio=πδytN:t0u~tN,λk/πδytN:t0u~tN,λk−1ρλλ~/ρλλk12: Computeα=min1,αratio13:ifα>U0,1then14:λ~=λk15:else16:λk=λ~17:end if18:end forThe same is repeated but with EnKBF in the place of KBF. Algorithm3 is the pseudocode showing the basic steps.Algorithm 3: EnKBF likelihood with MH.
Require:δt, ℵm, M, N, ut0ii=1M, πtn, and λ0.Ensure:λkk=1T.1:fork=1toNdo2: Drawλ~~Nλi,kcosϕ,ω/ℵmsinϕ∀i=1,2,⋯,2ℵm+1.3: ComputeCkx=expgx,λk.4:forn=1toT, δt>0do5:fori=1toMdo6:   Run a single-step EnKBF prediction ensembleu~tn+1i=utni+Ftn,λkutniδt7:   Run a single-step EnKBF analysis ensembleutn+1i=u~tn+1i+PtnHTtnR−1tndytn+εi−Htnu~tniδt8:end for9:  Compute prediction ensemble mean:u¯tn=1/M∑i=1Mu~tni.10:  Compute analysis ensemble mean:u^tn=1/M∑i=1Mutni.11:  Compute covariance:Ptn=1/M−1∑i=1Mutni−u^tnutni−u∧tnT.12:end for13:Metropolis-Hastings14: Computeαratio=πδytN:t0u¯tN,λk/πδytN:t0u¯tN,λk−1ρλλ~/ρλλk15: Computeα=min1,αratio16:ifα>U0,1then17:λ~=λk18:else19:λk=λ~20:end if21:end forResults for the first 2 parameters are shown in Figure3. The results in Figure 3 show that the EnKBF performs like the KBF filter. This agrees with the theory (see, for example, [17]) that the EnKBF yields optimal results in the limit M⟶∞. It is also noteworthy that the Metropolis-Hastings algorithm converges to the true parameter estimate. This can be seen in Figure 3(a), for example, where the filter estimate converges after about 100 parameter draws.Figure 3
Plots for velocity parameters for the first three parameters inλ obtained using KBF and EnKBF and run for 1000 time steps and for 1000 Metropolis-Hastings cycles. The number of particles used for EnKBF is M=1000, the time step size used in both filters is dt=0.01, μ=0.001, and 100 grid points are used. The plots indicate that the estimates, for both filters, converge after about 100 iterations.
(a)(b)We now look at the errors in the parameter estimates, the better to see the performance of the filters for the21 hyperparameters.Figure4(b) plots the boxplots showing the dispersion of parameter estimates resulting from the use of EnKBF and the root mean square errors for parameter estimates for both the EnKBF and KBF. The RMSE values are computed as follows.
(46)RMSE=1N∑k=1Nλi,k−λitrue2,where λi,k is the estimate of the ith parameter at the Metropolis-Hastings cycle k and λitrue is the trueith parameter.Figure 4
(a) Boxplot for the21 velocity parameters for the EnKBF run for 1000 time steps with 1000 Metropolis-Hastings cycles. A burn-in of 500 parameter draws is discarded. The stochastic advection equation model is used with the following settings: L=2π, 100 grid points, δt=0.01, M=1000 particles, μ=0.001, and localization radius of 10 grid points. (b) A plot of the root mean square error for the 21 hyperparameter estimates obtained using EnKBF and KBF. The plot indicates that the performance of EnKBF matches that of KBF in this setting.
(a)(b)The RMSE for both the KBF and EnKBF, as shown in Figure4(b), indicate that the performance of EnKBF matches that of the KBF for the 21 parameters. These heuristic results corroborate the theoretical findings. The boxplot, Figure 4(a), shows the dispersion of parameter samples in the case when EnKBF is used. The result indicates that the estimates match the true parameter values, with not so many outliers. This is indicative of the performance of not only EnKBF but also the Metropolis-Hastings procedure in locating the true parameter values and ensuring that no large excursions are made from the true parameter values.We now implement Algorithms2 and 3 for the discretised wave equation, (33).Example 2.
Wave equation.
We take the wave equation of Section5.2 and the associated initial conditions, Equations (30a) and (30b). The measurements are given by (42). The initial values of ut, βt, and ηt are uncorrelated. The aim is to estimate the spatially varying velocity, Cx, by means of filter likelihood and the Metropolis-Hastings algorithm.The discretisation of the wave equation is as shown in Section5.2. Figures 5(a) and 5(b) show the results.Figure 5
Plots for velocity parameters for the first three parameters inλ obtained using KBF and EnKBF and run for 1000 time steps and 1000 Metropolis-Hastings cycles. The number of particles used for EnKBF is M=1000, the time step size used in both filters is dt=0.01, μ=0.001, and 100 grid points are used. The plots indicate that the estimates, for both filters, converge after about 100 iterations.
(a)(b)The results in Figure5 indicate a close match in the performance of the EnKBF and KBF. This is as was anticipated in the theoretical findings, some of which are found in [17]. Notice also that the two filters converge to the true parameter values after a few parameter draws (about 50 in Figure 5(a) and 100 in Figure 5(b))—which is indicative of the robustness of the Metropolis-Hastings algorithm atop the EnKBF and KBF filters. The results also show that there are no wide excursions from the true parameter values, which testifies to the good performance of Algorithms 2 and 3.In Figure6 are plotted the boxplots showing the dispersion of the 21 hyperparameter estimates resulting from the use of EnKBF and the root mean square errors for parameter estimates for both the EnKBF and KBF.Figure 6
(a) Boxplot for velocity parameters for the EnKBF run for1000 time steps and 1000 Metropolis-Hastings cycles. A burn-in of 500 parameter draws is discarded. The stochastic wave equation model is used with the following settings: L=2π, 100 grid points, δt=0.01, M=1000 particles, μ=0.001, and localization radius of 10 grid points. (b) A plot of the root mean square error for the parameter estimates obtained using EnKBF and KBF. The plot indicates that the performance of EnKBF matches that of KBF in this setting.
(a)(b)The EnKBF elicits an optimal performance as can be seen in Figure6(b) where the RMSEs for both the EnKBF and the optimal filter (KBF) match for all the 21 hyperparameters. This also, as in the advection equation above, is in agreement with the theoretical findings that the EnKBF attains an optimal estimate in the limit M⟶∞. The boxplot, Figure 6(a), shows that the mean of EnKBF parameter estimates matches the true parameter values. What is more, there are not many outliers in the estimates. All these show that the EnKBF-Metropolis-Hastings algorithm is robust.In the next section, we apply the concepts on dual estimation to stochastic hyperbolic PDEs, which, in this case, are the advection and the wave equations.


## 8. Simultaneous Estimation of the State and Spatially Varying Parameters

The dual filter (see Section4.2 for details) can be adapted to allow for estimating both the state and spatially varying parameters, contemporaneously. The idea here is to replace the parameters in the dual filter with hyperparameters of the varying parameters to be approximated. The hyperparameters are then updated simultaneously and in parallel, with the state at each iteration, where one filter estimates the state and the other filter updates the hyperparameters—with each filter making use of the outcome of the other. To illustrate this argument, we turn to the advection and wave equation described in Examples 1 and 2, respectively, and use the KBF-EnKBF dual filter—in which the state is propagated and updated by means of the KBF whilst the hyperparameters are updated using EnKBF—and ENKBF dual filter—where both the state and the hyperparameters are propagated and updated using two EnKBFs running in parallel. The spatially varying velocity is as shown in Section 6.In the KBF-EnKBF dual filter, we update, for everyjth parameter particle λtj∈λtjj=1L, the state estimate, u^tj, using the KBF; that is,
(47a)du^tj=Ftu^tjdt+PtHTtR−1tdyt−Htu^tjdt,(47b)dPt=FtPtdt+PtFTtdt+GtGTtdt−PtHTtR−1tHtPtdt.The parameters are updated using the EnKBF; that is, each parameter hypothesis,λtj, is updated using
(48)dλtj=DtLHtR−1tdyt−0.5Htx^tj+Htx^tdt;t0≤t,where
(49a)λ^t=1L∑j=1Lλtj;t0≤t,(49b)DtL=1L−1∑i=1Lλtj−λ^tu∧tj−u∧tT;t0≤t,where
(49c)u^tj=1M∑i=1Muti,j,(49d)u^t=1L∑i=1Lu^tj.Algorithm4 gives a summary of the KBF-EnKBF dual filter.Algorithm 4: KBF-EnKBF dual filter.
Require:ut0j, λt0j, wt0j=1/L∀j∈1,2,⋯,L, Pt0, and δyt0,tT.Ensure:u^t0,tN, λ^t0,tN.1:forn=1toN, δt>0do2:forj=1toLdo3:  Updateu^tnj using Equations (47a) and (47b)4:  Update parametersλtnj using (48)5:end for6:  Computeλ^tn=1/L∑j=1Lλtnj7:  Computeu^tn=1/L∑i=1Lu^tnj8:end forThe EnKBF dual filter consists of an update ofM particles of the state, uti,j∈uti,ji,j=1M,L, for every parameter particle, λtj∈λtjj=1L, using the EnKBF; that is,
(50)duti,j=Ftuti,jdt+Gtdβti,j+PtMHTtR−1tdyt+R1/2tηti,j−Htuti,jdt,where ηti,j,t0≤t and βti,j,t0≤t are standard Brownian motion vector processes. The parameters are updated using the EnKBF given by (48). The summary of the EnKBF dual filter is given in Algorithm 5.Algorithm 5: EnKBF dual filter.
Require:ut0i,j, λt0j, wt0j=1/L∀i∈1,2,⋯,Mj∈1,2,⋯,L, Pt0, and δyt0,tT.Ensure:u^t0,tN, λ^t0,tN.1:forn=1toN, δt>0do2:forj=1toLdo3:fori=1toMdo4:   Calculateutni,j using (50)5:end for6:  Updateu^tnj using Equation (49c)7:  Update parametersλtnj using (48)8:end for9: Computeλ^tn=1/L∑j=1Lλtnj10: Computeu^tn=1/L∑i=1Lu^tnj11:end forFigures7(a) and 7(b) show the results for the first two parameters in (39) when the dual filters are applied to the advection equation.Figure 7
(a, b) are plots for velocity parameters for the first two parameters inλ obtained using KBF-EnKBF and EnKBF dual filters, applied to the advection equation, both run for 1000 time steps. The number of particles, for the state and hyperparameters, used in EnKBF is M=1000 and L=1000; the time step size used in both filters is dt=0.01. μ=0.001 and 100 grid points are used. Compared to the results obtained using the Metropolis-Hastings algorithm (Section 7), the dual filters register a better performance, at least in this example.
(a)(b)Evidently, from Figure7, the performance of the EnKBF and KBF-EnKBF dual filters almost matches. Both filters converge to the true estimate after a few iterations (about 100 in Figure 7(a)). This is another testament to the fact that EnKBF attains optimal estimates at high ensemble values. Moreover, both the EnKBF and KBF-EnKBF parameter estimates are not much spread as compared to the previous case where Metropolis-Hastings was used. This is more evident in the following results for the 21 hyperparameters estimated.In the following panels are plotted the root mean square errors for parameter estimates for both the EnKBF and KBF-EnKBF dual filters and the boxplots showing the dispersion of parameter estimates resulting from the use of the EnKBF dual filter applied to the advection equation.That the boxplots of EnKBF dual filter parameter estimates have short whiskers (see Figure8(a)) and few outliers and that the estimates match the true parameter values indicate that the EnKBF dual filter is robust in this setting. The EnKBF dual filter registers a slight variation in RMSE from the KBF-EnKBF dual filter. This indicates that the EnKBF, in this setting, performs optimally.Figure 8
(a, b) are, respectively, the boxplots showing the distribution of hyperparameter estimates of the EnKBF dual filter, applied to the advection equation, after a burn-in of500 iterations and the plots of the root mean square error for the hyperparameter estimates obtained using KBF-EnKBF and EnKBF (with different ensemble sizes) dual filters.
(a)(b)We now repeat the same procedure but with the wave equation described in Section5.2 in place of the advection equation. The following panels show the results. The results shown in Figure 9 are indicative of a dismal performance of the two dual filters—KBF-EnKBF and EnKBF dual filters—when applied to the wave equation as compared to the results obtained when the dual filters are applied to the advection equation (see Figure 7). We note that the wave equation is partially observed; that is, onlyuis observed in the discretised wave equation, (33), whereas the advection equation is fully observed. Furthermore, the number of unknowns in the state-parameter system of the wave equation is 300 whilst that in the advection equation is 200. These account for the dismal performance of the KBF-EnKBF and EnKBF dual filters when applied to the wave equation.Figure 9
(a, b) are plots for velocity parameters for the first two parameters inλ obtained using KBF-EnKBF and EnKBF dual filters, applied to the wave equation (Example 2), with both filters run for 1000 time steps. The number of particles, for the state and hyperparameters, used in the EnKBF filter is M=1000; the time step size used in both filters is dt=0.01. μ=0.001 and 100 grid points are used. Compared to the results obtained using the Metropolis-Hastings algorithm (Section 7), the dual filters register a dismal performance, at least in this example.
(a)(b)In Figure10 are plotted the root mean square errors for parameter estimates for both the EnKBF and KBF and the boxplots showing the dispersion of parameter estimates resulting from the use of EnKBF.Figure 10
(a, b) are, respectively, the boxplots showing the distribution of the21 hyperparameter estimates of the EnKBF dual filter after a burn-in of 500 iterations and plots of the root mean square error for the hyperparameter estimates obtained using KBF-EnKBF and EnKBF (with different ensemble sizes) dual filters, applied to the wave equation (Example 2).
(a)(b)


## 9. Conclusions

We have considered a special case of parameter estimation where the parameter to be estimated is spatially varying. Such a case arises, for example, in the velocity of a wave travelling through an inhomogeneous media. We proposed and studied two approaches: the use of filter likelihood and the Metropolis-Hastings procedure and joint estimation of state and parameters. The parameter is expressed as a Fourier series with constant coefficients. The coefficients are approximated and then substituted back to the Fourier series to obtain an approximation of the velocity. The Kalman-Bucy filter and the ensemble Kalman-Bucy filter are used. The filter likelihood with the Metropolis-Hastings procedure registers a better performance compared to the joint estimation procedure in both advection and wave equations. From the foregoing, Metropolis-Hastings with the filter evidence performs well in estimation of parameters compared to the dual filters—especially when the number of unknowns is large. This is indicative of the robustness of the Metropolis-Hastings algorithm in searching, and remaining in, the high probability region of the state-space.



---

*Source: 7909668-2023-01-10.xml*