# Embedded Estimation Sequential Bayes Parameter Inference for the Ricker Dynamical System

Authors: Muhammad Javvad Ur Rehman; Raheel Zafar; Hammad Dilpazir; Muhammad Farhan Sohail; Muhammad Arslan Khan; Rifat Jahan
Journal: Journal of Sensors
 (2022)

---

## Abstract

The dynamical systems are comprised of two components that change over time: the state space and the observation models. This study examines parameter inference in dynamical systems from the perspective of Bayesian inference. Inference on unknown parameters in nonlinear and non-Gaussian dynamical systems is challenging because the posterior densities corresponding to the unknown parameters do not have traceable formulations. Such a system is represented by the Ricker model, which is a traditional discrete population model in ecology and epidemiology that is used in many fields. This study, which deals with parameter inference, also known as parameter learning, is the central objective of this study. A sequential embedded estimation technique is proposed to estimate the posterior density and obtain parameter inference. The resulting algorithm is called the Augmented Sequential Markov Chain Monte Carlo (ASMCMC) procedure. Experiments are performed via simulation to illustrate the performance of the ASMCMC algorithm for observations from the Ricker dynamical system.


---

## Body

## 1. Introduction

It is possible to examine the dynamical characteristics of linear and nonlinear systems using state space modelling because it provides a consistent framework for doing so. The state space modelling process is divided into two stages: the development of a model that describes the underlying system dynamics over time, known as the state space model; and (ii) the development of a model that relates observations to state space variables, known as the measurement state model (also known as the measurement state model). The class of issues covered in this paper is the inference of unknown parameters, denoted by the symbolΩ, that regulate the dynamical system under consideration. As an example of a discrete-time stochastic system representing i and ii, the discrete-time stochastic system is expanded to include Ω to solve the parameter inference issue and is given by
(1)mk=Φk−1,Ωmk−1,uk−1,(2)yk=Ψk,Ωmk,vk,for k=1,2,3,⋯,T with m0~p0 represents the known prior distribution on the initial state space variable m0; in (1) and (2), k represents a generic time step and T denotes the final time, Φk−1,Ω· represents the transfer function in the state space model with noise uk, and Ψk,Ω· represents the nonlinear and non-Gaussian measurement model. Both Φk−1,Ω· and Ψk,Ω· are further known only up to the unknown parameter Ω∈Rp which is the object of inference. Additionally, in (1), uk−1 is assumed to be distributed according to FΩ· which is again parameterized by the components of Ω. In other words, Ω represents the collection of all unknown parameters in the dynamical system of (1) and (2), which are to be inferred based on the observed data. The development of Bayesian parameter inference methodology for the above dynamical system is the key objective of this article.Physical, biological, neuroscience, and object tracking systems are all examples of multidisciplinary fields in which the inference problems that govern the characteristics of nonlinear dynamical systems are applied [1–10]. The most commonly used statistical frameworks for finding solutions to parameter inference problems seem to be frequentist (or non-Bayesian) approaches that use maximum likelihood estimation (ML), such as the Expectation Maximization (EM) procedure [11–16]. It is computationally difficult to infer parameters for nonlinear dynamical systems using frequentist approaches. Furthermore, these procedures return just point estimates rather than the whole distribution, indicating parameter uncertainty. It is necessary to run the estimation method a large number of times, using uncertainty estimation techniques such as bootstrapping, to get estimates of uncertainty. This entails an increase in both the computational load and the intensity.The posterior distribution of unknown parameters, on the other hand, is used to offer both a point estimate and a corresponding estimate of uncertainty when estimating unknown parameters using Bayesian approaches. The aim of Bayesian procedures in a dynamical system setup is to obtain the posterior ofΩ, pΩy1:T, given all observations up to the final time T,y1:Y. The posterior pΩy1:T is generally difficult to obtain in closed form, so Bayesian computational algorithms are utilized to approximate it in the Monte Carlo sense. Markov chain Monte Carlo (MCMC) methods are the most frequently used approaches for estimating the posterior distribution when closed forms are not accessible, and they are the most accurate. Although Bayesian computational methods may be used to derive the posterior distribution of an unknown parameter, they are only justified in the limiting sense and can require a significant number of burn-ins before the MCMC algorithm converges [17–19]. In spite of this, certain computational methods are trivial to conceive and execute on computers, and standard packages are already available for many uncomplicated implementations of these procedures [20–22]. Some of the MCMC procedures are the Metropolis-Hastings (MH) algorithm, Gibbs sampler, and particle Markov Chain Monte Carlo (PMCMC) [23–25].Several sequential Monte Carlo (SMC) approaches have been developed to address the constraints of the Markov Chain Monte Carlo (MCMC) algorithm for parameter estimation in dynamical systems. It is the core concept of SMC to utilise important samples to estimate the posterior ofΩ at each k point in time and to propagate the samples sequentially via a suitable kernel. There exists an extensive literature on SMC methods (see, for example, [26–29]). The SMC-based parameter inference in nonlinear dynamical systems was first addressed in [30] where the Liu and West filter was developed. An artificial evolution of parameters for parameter Ω is used in the Liu and West filter and assumes a mixture of normal distribution for the posterior distributions, pΩy1:k, for k=1,2,⋯,T within the mixture distribution. The tuning parameters govern the extent of the control of overdispersion of the mixture components [30, 31]. To minimize the weight degeneracy or particle decay, the main idea is to generate new samples from the posterior by fitting the mixture to the posterior. The Liu and West filter can generally be applied to any dynamical system, which is the main attraction of this procedure. However, due to the artificial evolution of the unknown parameter, the artificial variability is incorporated, which is the main drawback of this algorithm.The aim of Bayesian procedures is to obtain the posterior ofΩ, pΩy1:T, given all observations up to time T,y1:T. The posterior pΩy1:T is again difficult to obtain in close form, so Bayesian computational algorithms are utilized to approximate it in a Monte Carlo sense.Another major class of SMC methods for parameter inference that does not introduce overdispersion in the posterior ofΩ is particle learning algorithms [32, 33]. The original method is attributed to Storvik [34, 35] resulting in Storvik’s filter (similar approaches are also proposed in [36, 37]). Storvik’s filter assumes that the posterior distribution of Ω given m0:k and y0:k depends on a lower-dimensional set of sufficient statistics that can be recursively updated for each k=1,2,⋯,T. This recursion for sufficient statistics is defined by sk+1=Ssk,mk+1,yk+1, leading to the generation of Ω samples according to Ω~pΩm0:k,y1:k=pΩsk for each k=1,2,⋯,T. Unlike the Liu and West filter, in Storvik’s filter, there is no artificial evolution process for Ω and thus it does not suffer from overdispersion [38]. However, the crucial assumption in Storvik’s filter is the availability of sufficient statistics sk as well as the ability of sampling from the posterior pΩsk given the sufficient statistics sk.Subsequent developments in SMC methods for parameter inference have extended the applicability of Storvik’s filter to a variety of more general settings (see, for example, [39, 40]). The extended Liu and West (ELW) filter to estimate parameters and states [40] divides the parameter set to be inferred Ω, into two sets, Ω0 and γ representing parameters without and with sufficient statistics, respectively. For the Ω0 set, the ELW filter uses the Liu and West filter, where an artificial random error is introduced to the static parameter Ω0. The set of parameters γ updated based on Storvik’s filter γ~pγskΩ0. The sufficient statistics skΩ0 are based on the static parameters. The rest of the parameters have the artificial evolution in which overdispersion is used. The overall set of parameters is represented as Ω=Ω0,γ. The ELW filter applies to a wider class of state space models compared to Storvik’s original procedure but suffers from two drawbacks, namely, (i) artificial overdispersion of the final posterior and (ii) the requirement of the existence of the sufficient statistic skΩ0 for γ and the ability of sampling from the posterior pγskΩ0.Development of statistical methods in specific ecological combines the nonlinear and near chaotic behavior of the system response for various applications [41–44]. A detailed comparison of the inference problem for nonlinear ecology and epidemiology is given in [45]. Nonlinearity is an observer in the experimental research [46]. Although the objectives of epidemiologists and ecologists are different, both are concerned about the persistence of specific species. The mathematical explanations of the population dynamics are similar in both studies [45].The aim and scope of this paper are to use the chaotic epidemiological or ecological model and perform parameter inference. There are two objectives. First is to perform the inference for the proposed application even if sufficient statistics are not available. Second is to use an online method to perform the parameter inference. Many researchers [47, 48] have discussed the relationship between statistics and chaos. The primary inference methodology used in this manuscript is developed in [49], which is a sequential MCMC (SMCMC) procedure to obtain the unknown parameter posterior inference in dynamical systems. However, then in [49], the proposed methodology was only applicable when the considered measurement model is linear and additive Gaussian noise. In this work, the measurement model of the Ricker dynamical system incorporates a non-Gaussian distribution, and the associated dynamical system (a special case of (1) and (2)) incorporates nonlinearities via the transfer function Φk−1,Ω·. The appropriate SMCMC algorithm for inferring Ω is developed for the Ricker dynamical system subsequently.The remainder of this paper is organized into the following sections: In Section2, Ricker’s model is discussed. In Section 3, details of the SMCMC procedure are given. The simulation experiments are given in Section 4. In the last section, Section 5, we state our conclusions, and potential future work is discussed.


## 2. Ricker’s Model

Theoretical ecology relies heavily on mathematical models of competition. Several mathematical models have been suggested to date to characterize the growth of contending populations; some of them are detailed in [50–52], including discrete-time models [53, 54]. The wide range of biological factors that influence ecosystem behavior makes it difficult for researchers to come to a consensus on how to simulate the dynamics of competing populations. Numerous instances of competing species and techniques for mathematical modelling are discussed in one of the pioneering publications on interspecific interaction [55]. The scramble competition has been found to fit the Ricker model [56, 57]. From order to chaos, the Ricker model illustrates dynamics [20–58]. It would be fascinating to observe what dynamic modes emerge when two Ricker maps are joined.The Ricker’s model is a classical discrete population model. It gives the expected density or numbers of individual species at each next generation.(3)m^k=rmk−1e−m∧k−1+uk−1,uk~N0,σ2.Ricker’s model is often used to explain the dynamics of two populations that are linked through migrations [50–58]. As a result, nothing is known about how Ricker communities evolved. In the experimental setup, a noisy model is considered in equation (3); here k represents the time evolution variable. r is a parameter to represent the intrinsic growth rate of the population. The process noise is represented by uk−1, which can also be considered environmental noise. Here the process noise has a zero mean, and the covariance parameter is σ2. Suppose the dynamics of the population are modelled with Ricker’s model, but it is not possible to know the exact population density at any time, so it is necessary to have a measurement model.
(4)yk~Poissonϕmk.Hereyk is the measurement parameter of the individual sampled at any time point k, and ϕ is the scale parameter.By transformingmk=logm^k, it is seen that (3) becomes
(5)mk=logr+logmk−1−mk−1+zk,and the measurement model (4) becomes
(6)yk~Poissonϕemk.Using the transformed variables, equations (5) and (6) can be seen to be special forms of the general state space and measurement model equations given by (1) and (2), respectively. We have Φk−1,Ω·=logr+logmk−1−mk−1, and Ψk,Ω· is the Poisson probability density function with mean ϕemk. Thus, the Ricker dynamical system has three underlying parameters that govern the system, namely, σ2, logr, and ϕ. In this paper, the latter two parameters are taken to be unknown, that is, Ω=logr,ϕ and σ2=0.09 is assumed fixed and known. True values of the parameters for our simulation studies are taken to be Ω=logr,ϕ=3.8,0.7, which the same as given in [45].


## 3. An Augmented Sequential Markov Chain Monte Carlo Algorithm

Sequential Markov chain Monte Carlo algorithm working process is described in the subsequent text. Sequential updating characteristic of SMCMC is used in the proposed technique. The core notion is addressed in [59], which is focused on the Monte Carlo sum, but the cumulative filtering steps required for estimating the probability pyky1:k−1,Ω minimise. State variable mk and unknown parameter Ω are supplemented at the time step of k−1 in [49]. Therefore, the likelihood function changes accordingly from pyky1:k−1,Ω to pykmk−1,Ω. It is possible to achieve the analytical expression to avoid the need for cumulative filtering procedures. This is particularly useful when the amount of k is substantial.The ASMCMC is an iterative procedure that starts from the initial time stepk=1, increases sequentially, and finally ends when k=T. Within the k-th step of the ASMCMC procedure, posterior samples are obtained based on an underlying Markov chain Monte Carlo procedure. This Markov chain Monte Carlo procedure will be called the k-th time step Markov chain Monte Carlo procedure. The target of the k-th step MCMC procedure is the posterior, pΩ,mTy1:k, at time step k. Thus, when k=T, we will have obtained samples from the desired posterior pΩ,mTy1:k as well as pΩ,y1:k via marginalization. Details on the implementation of the k-th step Markov chain Monte Carlo procedure are as follows.Assume that after thek−1-th time step MCMC procedure, M samples (i.e., particles) ξk−1j≡Ωk−1j,mk−1j,j=1,2,⋯,M are available from the posterior density pΩ,mk−1y1:k−1. A class of Gaussian mixture model (GMM) is fitted using ξk−1j,j=1,2,⋯,M following the methodology outlined in [48–60]. This results in the estimated density pMixΩ,mk−1y1:k−1 based on GMMs. The methodology developed and implemented in [48–60] ensures that the fitted GMM is close to the true posterior density pΩ,mk−1y1:k−1, that is,
(7)pMixΩ,mk−1y1:k−1≈pΩ,mk−1y1:k−1,even though the exact form of the latter is unavailable.To initialize thek-th time step MCMC procedure, the M samples ξk−1j,j=1,2,⋯,M are taken to form the starting points of Mk-th step MCMC procedures, that is, Ωk−1j≡Ωkj,0 and mk−1j≡mk−1,kj,0 for j=1,2,⋯,M. Some notations are developed here: denoted by Ωgj,k and mk−1,kj,g, respectively, to be the values of Ω and mk−1 at the g-th cycle of the k-th time step MCMC procedure initialized by kik−10 for j=1,2,⋯,M. In other words, the Mk-th time step MCMC procedures initialized based on Ωkj,0 form separate chains based on separate starting values; this entails that the M MCMC chains can be run in parallel for each time step k=1,2,⋯,T.For details of the steps involved within eachk-th step MCMC procedure, we suppress the j notation. For a generic k-th step MCMC procedure (note that there are M of them), initialize ξk0≡Ωk0,mk−1,k0 as above from outputs of the k−1-th chain. At step k of the kth step MCMC procedure, assume Ωkg,mk−1,kg is already available. To transit from g⟶g+1,
(i)
Generate(8)Ωk∗,mk−1,k∗~qΩk,mk−1,kΩkg,mk−1,kg,where q·· in equation (8) is proposal density.
(ii)
Compute the acceptance probability(9)αkΩg,mg,Ω∗,m∗=minA·B,1here(10)A=pykmk−1,k∗,Ωk∗p^MixΩk∗,mk−1,k∗y1:k−1pykmk−1,kg,Ωkgp^MixΩkg,mk−1,kgy1:k−1(11)B=qΩkg,mk−1,kgΩk∗,mk−1,k∗qΩk∗,mk−1,k∗Ωkg,mk−1,kgwhere pykmk−1∗,Ωk is the probability function of yk given mn−1 and Ω as in ((13)).SetΩkk+1,mk−1,kk+1=Ωk∗,mk−1,k∗ with probability αkΩkg,mk−1,kg,Ωk∗,mk−1,k∗; otherwise, set
(12)Ωkk+1,mk−1,kk+1=Ωkg,mk−1,kg.Continue the iteration fromk+1⟶k+2.The aforementioned Markov chain will converge ask⟶∞ to the stationary (and target) distribution determined by the numerator of the expression A in (10), namely, pykmk−1,Ωp^MixΩ,mk−1y1:k−1. Since
(13)pykmk−1,Ωp^MixΩ,mk−1y1:k−1≈pykmk−1,ΩpΩ,mk−1y1:k−1∝pΩ,mk−1y1:k.The target density of thek-th step MCMC procedure is the posterior density of Ω,mk−1 given y1:k, where the first approximate equality is due to ((7)). In our previous work [49], pykmk−1,Ω was available in closed form due to a linear measurement model and additive noise variables distributed as Gaussian in both the state space and measurement models. In the present context,
(14)pykmk−1,Ω=∫mkΨk,Ωykmkpmkmk−1,Ωdmk.Subsequently, after the burn-in periodB, samples of Ωkg∗ only are collected for a large value of g∗≥B, whereas the samples of mk−1g∗ are discarded. MCMC theory and the marginalization property ensure that Ωkg∗ is distributed according to the posterior pΩy1:k. Since there are M parallel MCMC chains, M such samples of Ω, namely, Ωkg∗,jj=1M, are collected in this way. Subsequently, pure filtering steps as in given in ((15)) are performed to get the samples of mk from pmky1:k,Ω for each Ω in the collection Ωkg∗,jj=1M [48]. Hence, the particle pairs Ωkg∗,j,mkjj=1M constitute M samples from the joint posterior density pΩ,mky1:k and can serve as the input for the next n+1-st time step of the ASMCMC Algorithm 1.
(15)wkj∝pykmkj,Ωpmkjmk−1j,Ωqmkjmk−1j,yk,Ωforj=1,2,⋯,M.Algorithm 1: Augmented sequential MCMC procedure.
mTj,ΩTjj=1M=ASMCMCm0j,Ω0jj=1MInitialize withM initial value Ω0j,m0j prior densityp0Ω,m0.Start for k=1:Tmkj,Ωkjj=1M=AMCMCmk−1j,Ωk−1jj=1M,ykDO:M Parallel MCMC chains initialized from M pairs Ωk−1j,mk−1j, j=1,2,⋯,M.For the genericj-th chain:Start withΩ0,m0≡Ωk−1i,mk−1i (prior value).for k=1:B (burn-in chains)Generate samplesΩ∗,m∗~qΩ∗,m∗Ωk−1,mk−1.Compute the acceptance probabilityαk=minA·B,1asin9and10.SetΩg,mg=Ω∗,m∗ withprobabilityαkElseΩg,mg=Ωk−1,mk−1with probability1−αkend forObtainΩB≡Ωkj, from j=1,2,⋯,Mchains.mkj~pmkjy1:k,Ωkj Samples are based(12)Mixture fittingp^MixΩ,mky1:kOutput:Ωkj,mkji=1M and fitted mixture model p^Mix,n≡p^MixΩ,mky1:k.The posterior distributionpΩy1:T is based on the collected samples of ΩTj.end for


## 4. Results and Discussion

In this section, the ASMCMC methodology is used for parameter inference in the Ricker dynamical system. They-observations were generated from the time-discretized Ricker’s model (state space model given by (5) and the measurement model given by (6)) using the initial point mass prior at m0=7. Starting from the initial state values generated from the prior, the state and measurement systems are updated at every fixed time step of k=1,2,⋯,T=10. The true value for Ω is taken as Ω=logr,ϕ=3.8,0.7, which is the same as the choice made in [45]. The burn-in B was set at B=5000. The estimated posterior density curves are obtained based on the final posterior samples of Ω=logr,ϕ at the completion of the ASMCMC algorithm. These density curves are given in Figures 1 and 2 for the parameters logr and ϕ, respectively. In Figure 1, estimated density curves based on final posterior samples of the parameter logr at the completion of the ASMCMC are presented. The vertical solid black line represents the true value of the parameter logr=3.8, whereas the estimated density curves based on the final posterior samples of the parameter ϕ at the completion of the ASMCMC are presented in Figure 3. The vertical solid black line represents the true value of the parameter ϕ=0.7. We note that the true values of the parameters are well within the support of their respective posterior densities, which gives credence to the parameter inference methodology using the ASMCMC procedure. The following Table 1 represents the simulation parameters considered in the experimental setup.Figure 1
The estimated density curves are based on final posterior samples of the parameterlogr at the completion of the ASMCMC.Figure 2
The estimated density curves are based on final posterior samples of the parameterϕ at the completion of the ASMCMC.Figure 3
Boxplot trajectories representing the distribution ofpΩy1:k based on the ASMCMC algorithm. The boxplots are constructed based on the M final iterates logrkj,B.Table 1
Simulation parameters.
VariableValueDescriptionlogr3.8True value of parameterϕ0.7True value of parameterm07Prior of state variableB5000Burn-in periodTo illustrate the convergence of the ASMCMC procedure, boxplot trajectories represent the distribution ofpΩy1:k for k=1,2,⋯,T can be plotted. These boxplots are constructed based on the M final iteration Ωkg∗,j for j=1,2,⋯,M, collected after burn-in. Figure 3 shows that the boxplot trajectories based on posterior samples of logr. Figure 3 indicates that these boxplot trajectories have stabilized long before reaching the final time point T=10. In other words, the posterior plogry1:k changes much when k approaches the final time point T, which can be taken as an indication that the estimate of plogry1:k based on the ASMCMC sampler has converged. A similar boxplot trajectory plot (shown in Figure 4) is obtained for ϕ which indicates that the convergence of its posterior distribution has been achieved.Figure 4
Boxplot trajectories representing the distribution ofpΩy1:k based on the ASMCMC algorithm. The boxplots are constructed based on the M final iterates ϕkj,B.


## 5. Conclusion

The dynamical system based on Ricker’s model is an example of a nonlinear and non-Gaussian system that has applications in ecology and epidemiology. We develop the ASMCMC algorithm for the Ricker dynamical system to perform Bayesian parameter inference. We observe that the posteriors encompass the true parameter values used to simulate observations from the Ricker dynamical system. Our future work will be to investigate the performance of the ASMCMC algorithm whenσ2 is unknown as well and for high-dimensional dynamical systems that appear in ecology and epidemiology.



---

*Source: 4540366-2022-06-08.xml*