Variational State and Parameter Estimation

1 Introduction

State-space models are a widely utilised and flexible class of models considered in many scientific and engineering applications [15]. Generally, identification of the parameters of state-space models remains a challenging task [16, 19]. In this paper, nonlinear state-space models of the form


$x_{k + 1}$	$\sim p (x_{k + 1} ∣ x_{k}, θ),$	(1a)
$y_{k}$	$\sim p (y_{k} ∣ x_{k}, θ),$	(1b)

are considered, where $x_{k} \in R^{n_{x}}$ and $y_{k} \in R^{n_{y}}$ are the state and measurement at time step $k$ , respectively, and $θ \in R^{n_{θ}}$ are model parameters. For clarity the presence of an input $u_{k} \in R^{n_{u}}$ is not shown, it is, however, allowed.

In this paper, given a sequence of measurements, $y ≜ y_{1 : T} ≜ {y_{1}, \dots, y_{T}}$ , the estimation of the state and model parameters is considered. Using a Bayesian approach, this amounts to obtaining the distribution $p (x, θ ∣ y)$ , given by

(2)

where $x ≜ x_{1 : T + 1} ≜ {x_{1}, \dots, x_{T + 1}}$ , and $p (θ)$ is a prior distribution over the parameters.

From $p (x, θ ∣ y)$ , distributions of both parameters and states can be obtained via marginalisation. For example, the marginal distribution $p (θ ∣ y)$ is found via

p (θ ∣ y) = \int p (x, θ ∣ y) d x .

(3)

The difficulty with both (2) and (3) is that they are generally intractable. This is due to both the arbitrary form of the distributions considered and the nonlinear integrals involved. As such, these distributions must be approximated. Many differing approximations exist, one class of approaches is to utilise particle-based Monte Carlo methods, an example of this is Hamiltonian Monte Carlo (HMC) [18, 8] methods.

In this paper, an alternative approach, known as variational inference (VI) [10, 2], to approximate $p (x, θ ∣ y)$ is followed. The contribution of this paper is to use VI to approximate $p (x, θ ∣ y)$ . Using a Gaussian assumed density, and a suitable approximation for the integrals, this provides a deterministic method of approximating $p (x, θ ∣ y)$ that can be efficiently performed using gradient-based optimisation with exact first- and second-order derivatives.

2 Variational Inference

Variational inference is a widely used method to approximate, potentially intractable, posterior distributions with a parametric density of an assumed form. The numeric values of this parametrisation are obtained by maximising a likelihood lower bound [2]. In this section, the use of VI in approximating $p (x, θ ∣ y)$ is examined.

As $p (x, θ ∣ y)$ is intractable, it is approximated using an assumed density, parameterised by $β$ , and denoted as $q_{β} (x, θ)$ . The quality of this approximation is measured by

KL [q_{β} (x, θ) ∣∣ p (x, θ ∣ y)],

(4)

the Kullback-Leiber (KL) [13] divergence of $p (x, θ ∣ y)$ from $q_{β} (x, θ)$ . As it is desired that $q_{β} (x, θ)$ is close to $p (x, θ ∣ y)$ the aim is to find a $β$ , denoted $β^{⋆}$ , which minimises this KL divergence and is given by

β^{⋆} = arg min β KL [q_{β} (x, θ) ∣∣ p (x, θ ∣ y)] .

(5)

However, as $p (x, θ ∣ y)$ is intractable $KL [q_{β} (x) ∣∣ p (x, θ ∣ y)]$ cannot be directly minimised or even evaluated. Instead, the concept of VI is used to, equivalently, maximise a likelihood lower bound.

Using conditional probability, this lower bound can be found by expressing the log-likelihood,

log p (y)

, as

log p (y)

= log p (x, θ, y) - log p (x, θ ∣ y) .

(6)

Through addition and subtraction of $log q_{β} (x, θ)$ to the right-hand side of (6) this leads to

log p (y)

= log \frac{p (x, θ, y)}{q_{β} (x, θ)} + log \frac{q_{β} (x, θ)}{p (x, θ ∣ y)} .

(7)

As $log p (y)$ is independent of $x$ and $θ$ , the log-likelihood can alternatively be given by

log p (y) = \int q_{β} (x, θ) log p (y) d τ,

(8)

where $τ = [x^{T}, θ^{T}]^{T}$ .

Substituting (7) into the right-hand side of (8) arrives at

	$log p (y) =$	$\int q_{β} (x, θ) log \frac{p (x, θ, y)}{q_{β} (x, θ)} d τ$
		$+ \int q_{β} (x, θ) log \frac{q_{β} (x, θ)}{p (x, θ ∣ y)} d τ,$		(9)

which will be expressed as

log p (y) = L (β) + KL [q_{β} (x, θ) ∣∣ p (x, θ ∣ y)],

(10)

where

L (β) = \int q_{β} (x, θ) log \frac{p (x, θ, y)}{q_{β} (x, θ)} d τ .

(11)

From (10), and as any KL divergence is non-negative, $L (β)$ is a lower bound to the log-likelihood $log p (y)$ .

Furthermore, from (10) we have

KL [q_{β} (x, θ) ∣∣ p (x, θ ∣ y)] = log p (y) - L (β),

(12)

and therefore


$β^{⋆}$	$= arg min β KL [q_{β} (x, θ) ∣∣ p (x, θ ∣ y)]$	(13a)
	$= arg min β log p (y) - L (β)$	(13b)
	$= arg max β L (β) .$	(13c)

The important property of this is $L (β)$ does not involve the intractable distribution, $p (x, θ ∣ y)$ , rather only the assumed distribution $q_{β} (x, θ)$ which is selected in a tractable form.

The desired approach is therefore to find a value for $β$ , that maximises $L (β)$ via

β^{⋆}

= arg max β L (β),

(14)

for an assumed density form. This provides an approximation for $p (x, θ ∣ y)$ given by $q_{β} (x, θ)$ as desired.

However, as given, this approach is computationally intractable for any large quantity of time steps. This is because $q_{β} (x, θ)$ is a full dense distribution, the dimension of which continues to grow. Fortunately, by expressing $L (β)$ as

$L (β)$	$= \int q_{β} (x, θ) log \frac{p (x, θ, y)}{q_{β} (x, θ)} d τ$
	$= \int q_{β} (x, θ) log p (x, θ, y) d τ$
	$- \int q_{β} (x, θ) log q_{β} (x, θ) d τ,$	(15)

and utilising the Markovian nature of state-space models (see Appendix A for details) to give


	$= \int q_{β} (x_{1}, θ) log p (x_{1}, θ) d [x_{1}, θ]$
	$+ T \sum k = 1 \int q_{β} (x_{k : k + 1}, θ) log p (x_{k + 1}, y_{k} ∣ x_{k}, θ) d [x_{k : k + 1}, θ]$
	$= I_{1} (β) + I_{23} (β),$		(16)

and

	$\int q_{β} (x, θ) log q_{β} (x, θ) d τ$
	$= T \sum k = 1 \int q_{β} (x_{k + 1}, x_{k}, θ) log q_{β} (x_{k + 1}, x_{k}, θ) d [x_{k : k + 1}, θ]$
	$- T \sum k = 2 \int q_{β} (x_{k}, θ) log q_{β} (x_{k}, θ) d [x_{k}, θ]$
	$= I_{4} (β),$		(17)

we have

(18)

without introducing any approximations.

Therefore, calculation of $L (β)$ does not require the full distribution, rather only each $q_{β} (θ, x_{k}, x_{k + 1})$ and as such remains tractable to calculate for large quantities of time steps.

The focus of this paper correspondingly shifts from the, intractable, aim of obtaining $q_{β} (x, θ)$ , to the tractable aim of obtaining each $q_{β} (θ, x_{k}, x_{k + 1})$ .

3 Assumed Gaussian Distribution and Tractable Approximation

In this section, the parametric assumed density chosen and the subsequent formation of a standard form optimisation problem is examined. When selecting the parametric form of $q_{β} (θ, x_{k}, x_{k + 1})$ there are two conflicting goals. Firstly, the parametrisation chosen should be sufficiently flexible to well represent $p (θ, x_{k}, x_{k + 1})$ . Secondly, the parametrisation should be chosen to allow the optimisation to be efficiently performed.

In this paper, we have chosen to use a multivariate Gaussian parameterised as

q_{β_{k}} (θ, x_{k}, x_{k + 1}) = N ⎛ ⎜ ⎝ ⎡ ⎢ ⎣ \begin{matrix} θ x_{k} x_{k + 1} \end{matrix} ⎤ ⎥ ⎦; ⎡ ⎢ ⎣ \begin{matrix} μ_{k} {¯ μ}_{k} {~ μ}_{k} \end{matrix} ⎤ ⎥ ⎦, P_{k}^{\frac{T}{2}} P_{k}^{\frac{1}{2}} ⎞ ⎟ ⎠,

(19)

where $μ_{k} \in R^{n_{θ}}$ , ${¯ μ}_{k} \in R^{n_{x}}$ , ${~ μ}_{k} \in R^{n_{x}}$ , $P_{k}^{\frac{1}{2}} \in R^{(n_{θ} + 2 n_{x}) \times (n_{θ} + 2 n_{x})}$ , and

P_{k}^{\frac{1}{2}} = ⎡ ⎢ ⎣ \begin{matrix} A_{k} & B_{k} & C_{k} 0 & D_{k} & E_{k} 0 & 0 & F_{k} \end{matrix} ⎤ ⎥ ⎦,

(20)

where $A_{k}, D_{k}, F_{k}$ are upper triangular. The parameter $β_{k}$ is therefore given as

β_{k} = {μ_{k}, {¯ μ}_{k}, {~ μ}_{k} A_{k}, B_{k}, C_{k}, D_{k}, E_{k}, F_{k}},

(21)

and related to $β$ via

β = {β_{1}, β_{2}, \dots, β_{T}} .

(22)

Importantly, parameterising each joint distribution using the Cholesky factor,

P_{k}^{\frac{1}{2}}

, allows the components of

I_{23} (β)

corresponding to each time step to be easily approximated using Gaussian quadrature as


	${^I}_{23_{k}} (β)$	$= n_{s} \sum j = 1 w_{j} log p (x_{k + 1}^{j}, y_{k} ∣ x_{k}^{j}, θ_{k}^{j}),$		(23a)

where $w_{j}$ is a weight, and the $n_{s}$ sigma points are denoted as $ξ^{j} \in R^{n_{θ} + 2 n_{x}}$ where $ξ^{j} = {[\begin{matrix} (θ_{k}^{j})^{T}, & (x_{k}^{j})^{T}, & (x_{k + 1}^{j})^{T} \end{matrix}]}^{T}$ is given by linear combinations of the joint mean ${[\begin{matrix} μ_{k}^{T}, {¯ μ}_{k}^{T}, {~ μ}_{k}^{T} \end{matrix}]}^{T}$ and the columns of $P_{k}^{\frac{T}{2}}$ . The sigma points, $ξ^{j}$ , being linear combinations of the elements of $β_{k}$ is important as it significantly simplifies the calculation of both first- and second-order derivatives used in the optimisation.

Furthermore, due to the Gaussian assumption $I_{4} (β)$ , and its derivatives, are available in closed form. With a suitable prior $I_{1} (β)$ is also available in closed form, alternatively it can be similarly approximated using Gaussian quadrature.

This allows $L (β)$ to be approximated by $^L (β)$ as

^L (β)

= I_{1} (β) + {^I}_{23} (β) - I_{4} (β) .

(24)

This parametrisation is, however, over-parameterised as the marginal distributions $q_{β} (θ_{k})$ , $q_{β} (x_{k})$ and $q_{β} (θ, x_{k})$ can be calculated from either $q_{β} (θ_{k + 1}, x_{k - 1}, x_{k})$ or
$q_{β} (θ_{k}, x_{k}, x_{k + 1})$ . Hence, the constraints defined by the feasible set

¯ Ω := {β \in R^{n_{β}} ∣ {¯ c}_{k} (β) = 0, k = 1, \dots, T - 1},

(25)

where $n_{β} = T (n_{θ} + 2 n_{x} + \frac{(n_{θ} + 2 n_{x}) (n_{θ} + 2 n_{x} + 1)}{2})$ and

{¯ c}_{k} (β) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} μ_{k + 1} - μ_{k + 1} {~ μ}_{i + 1} - {¯ μ}_{i} A_{k}^{T} A_{k} - A_{k + 1}^{T} A_{k + 1} A_{k}^{T} C_{k} - A_{k + 1}^{T} B_{k + 1} C_{k}^{T} C_{k} + E_{k}^{T} E_{k} + F_{k}^{T} F_{k} - B_{k + 1}^{T} B_{k + 1} - D_{k + 1}^{T} D_{k + 1} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦,

are required to ensure consistency.

However, as each $A_{k}$ , being a Cholesky factor of a covariance matrix, is non-singular the constraints set can be simplified to

Ω := {β \in R^{n_{β}} ∣ c_{k} (β) = 0, k = 1, \dots, T - 1},

(26)

where

c_{k} (β) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} μ_{k + 1} - μ_{k + 1} {~ μ}_{i + 1} - {¯ μ}_{i} A_{k} - A_{k + 1} C_{k} - B_{k + 1} C_{k}^{T} C_{k} + E_{k}^{T} E_{k} + F_{k}^{T} F_{k} - B_{k + 1}^{T} B_{k + 1} - D_{k + 1}^{T} D_{k + 1} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .

Using this, an approximation of (14) is given by


$^β = arg max β$	$^L (β),$	(27a)
s.t.	$β \in Ω .$	(27b)

This constrained optimisation problem is of a standard form and can be efficiently solved using exact first- and second-order derivatives without more approximations to find a local maximum. The calculation of derivatives required can be performed using automatic differentiation [7] without manual effort. Due to the parametrisation chosen, this can be efficiently performed using standard tools. For the numerical examples provided CasADi [1] has been used.

The result of this optimisation is numerical values for each

q_{{^β}_{k}} (θ, x_{k}, x_{k + 1}) = N ⎛ ⎜ ⎝ ⎡ ⎢ ⎣ \begin{matrix} θ x_{k} x_{k + 1} \end{matrix} ⎤ ⎥ ⎦; ⎡ ⎢ ⎣ \begin{matrix} μ_{k} {¯ μ}_{k} {~ μ}_{k} \end{matrix} ⎤ ⎥ ⎦, P_{k}^{\frac{T}{2}} P_{k}^{\frac{1}{2}} ⎞ ⎟ ⎠,

(28)

which approximates each $p (θ, x_{k}, x_{k + 1} ∣ y)$ . Due to the Gaussian form selected, marginal and conditional distributions of both the parameters and states can be readily extracted as desired.

The proposed approach is summarised in Algorithm 1.

0: Measurements

y_{1 : T}

, prior

p (x_{1}, θ)

, initial estimate of

β

0: Each

q_{{^β}_{k}} (θ, x_{k}, x_{k + 1})

for

k \in 1, \dots, T

- Obtain

^β

from (27) initialised at

β

- Extract each

{^β}_{k}

from

^β

Algorithm 1 State and Parameter Estimation

4 Numerical Examples

In this section, we provide two numerical examples of the proposed approach and compare the results against HMC. The HMC results have been calculating using STAN [3] as detailed in [8].

4.1 Stochastic Volatility

In this example, estimation of $θ = [a, b, log (\sqrt{c})]$ for the following stochastic volatility model:


$x_{k + 1}$	$= a + b x_{k} + \sqrt{c} w_{k},$	(29a)
$y_{k}$	$= \sqrt{e^{x_{k}}} v_{k},$	(29b)

where $w_{k} \sim N (0, 1)$ and $v_{k} \sim N (0, 1)$ is considered using 726 simulated measurements.

For the proposed method fifth-order cubature [9] is used to perform the Gaussian quadrature approximations and the optimisation routine is run until convergence to a locally optimal point requiring $17.6 s$ . The HMC results are obtained using $8 \times 10^{3}$ iterations with half discarded at warm-up requiring $87.0 s$ .

In Figure 1 the marginal distributions of $θ$ using both the proposed approach and HMC are shown. While the marginal distributions for parameters $a$ and $b$ closely match, the VI-based method is overconfident in the distribution of $log (\sqrt{c})$ and cannot closely match the non-Gaussian relationships between the parameters.

Figure 1: Distribution of $θ$ for both the proposed approach (contour plots) and HMC (scatter plots / histograms). The true parameters are shown as vertical lines.

Despite this the state distributions obtained using both methods closely match. Figure 2 highlights this for three differing time steps, which are representative of the results for other time steps.

Figure 2: State distributions for three time steps using the proposed approach (curve) and HMC (histogram) illustrating the high level of similarity.

In this example, the proposed method is demonstrated on a nonlinear and non-Gaussian example and is, generally, representative of the results obtained using HMC at a lower computational cost. This, however, comes with a trade-off, when the underlying distribution is non-Gaussian, the Gaussian assumption of the proposed method may produce overconfident results as shown for the distribution of $log (\sqrt{c})$ .

4.2 Inverted Pendulum

In this section, the proposed method is applied to a rotational inverted pendulum, or Furata pendulum [6]

, using real data collected from a Quanser QUBE—Servo 2. The estimation of a four-dimensional state vector, six parameters

ϕ = [J_{r}, J_{p}, K_{m}, R_{m}, D_{p}, D_{r}]^{T}

, and a seven dimensional covariance matrix

Π

associated with coupled additive Gaussian noise of the process and measurement models is considered using 375 measurements. For further details of the model refer to Appendix B.

For the proposed approach, a third order unscented-transform [11, 23] with $α = 1$ , $κ = 0$ , and $β = 0$ is used to perform the Gaussian quadrature required. A point estimate for $Π$ is also considered, this highlights the flexibility of the proposed approach which allows for distributions of subsets of parameters alongside point estimates of others. In this example, we also exploit the additive form of the state-space model, this allows for a Hessian to be efficiently approximated using only Jacobian’s. Identification problems for systems with additive Gaussian noise commonly use this approach [12, 20]. While the calculations for this example are not fully detailed in Section 3, they represent minor modifications of the resultant optimisation problem.

For the HMC results, $10 \times 10^{3}$ iterations are used to obtain distributions for all states and parameters. As the HMC results calculate a distribution for each element of $Π$ , while the proposed method considers a point estimate, it is expected the results of the proposed method and HMC slightly differ.

In Figure 3 the distribution of $ϕ$ using both the proposed approach and HMC are shown. This illustrates that, despite the approximations introduced, the proposed method has delivered results close to HMC. This has been achieved with a significant computational reduction. The proposed method required a few minutes to calculate while the HMC results required a few hours.

Figure 3: Distribution of inverted pendulum parameters calculated using the proposed method (contour plots) and HMC (scatter plots / histograms).

5 Conclusion

In this paper, a variational based approach to estimating state and parameter distributions is presented for general nonlinear state-space models. The proposed method results in an optimisation problem of a standard form for which, due to the parametrisation selected, derivatives are readily available. While the proposed method assumes a Gaussian density, and is therefore limited compared to more general HMC approaches, the numerical examples illustrate that approximations suitable for many purposes are obtained on both simulated and real data. Compared to state-of-the-art HMC approaches this is achieved at a significantly lower computational cost.

Future research directions for this work include blending the proposed approach with sequential Monte Carlo (SMC) algorithms along the lines of what has been done in e.g. [14, 17]. The motivation for this is to combine the computational efficiency of the proposed method and the theoretical guarantees of SMC to overcome the limitations imposed by the Gaussian assumed densities.

References

[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi: a software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, July 2018.
[2] D. Blei, A. Kucukelbir, and J. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859–877, 2017.
[3] Bob Carpenter, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A Probabilistic Programming Language. Journal of Statistical Software, 76(1), 2017.
[4] S. B. Chitralekha, J. Prakash, H. Raghavan, R. B. Gopaluni, and S. L. Shah.
Comparison of Expectation-Maximization based parameter estimation using Particle Filter, Unscented and Extended Kalman Filtering techniques.
IFAC Proceedings Volumes, 42(10):804–809, 2009.
[5] Jarrad Courts, Adrian Wills, and Thomas B. Schön. Variational Nonlinear State Estimation. arXiv, 2020.
[6] K. Furuta, M. Yamakita, and S. Kobayashi. Swing-up Control of Inverted Pendulum Using Pseudo-State Feedback. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 206(4):263–269, November 1992.
[7] Andreas Griewank and Andrea Walther. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. Society for Industrial and Applied Mathematics, jan 2008.
[8] Johannes Hendriks, Adrian Wills, Brett Ninness, and Johan Dahlin. Practical Bayesian System Identification using Hamiltonian Monte Carlo. arXiv, 2020.
[9] Bin Jia, Ming Xin, and Yang Cheng. High-degree cubature Kalman filter. Automatica, 49(2):510–518, feb 2013.
[10] Michael I. Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul. An Introduction to Variational Methods for Graphical Models. Machine Learning, 37(2):183–233, November 1999.
[11] Simon J. Julier and Jeffrey K. Uhlmann. New extension of the Kalman filter to nonlinear systems. In Ivan Kadar, editor, Signal Processing, Sensor Fusion, and Target Recognition VI. SPIE, July 1997.
[12] Juho Kokkala, Arno Solin, and Simo Särkkä. Sigma-Point Filtering and Smoothing Based Parameter Estimation in Nonlinear Dynamic Systems. Journal of Advances in Information Fusion, 2015.
[13] S. Kullback and R. A. Leibler. On Information and Sufficiency. The Annals of Mathematical Statistics, 22(1):79–86, March 1951.
[14] Fredrik Lindsten, Jouni Helske, and Matti Vihola. Graphical model inference: Sequential Monte Carlo meets deterministic approximations. In Advances in Neural Information Processing Systems 31, pages 8190–8200. 2018.
[15] L. Ljung. System Identification: Theory for the User. Prentice Hall information and system sciences series. Prentice Hall PTR, 1999.
[16] Lennart Ljung. Perspectives on system identification. Annual Reviews in Control, 34(1):1–12, April 2010.
[17] Christian A. Naesseth, Fredrik Lindsten, and David M. Blei. Markovian Score Climbing: Variational Inference with KL(p||q). In Advances in Neural Information Processing Systems 33. 2020.
[18] Radford M. Neal.
Handbook of Markov Chain Monte Carlo
, chapter MCMC using Hamiltonian dynamics, pages 113–162. Chapman and Hall/CRC, 2011.
[19] Brett Ninness. Some System Identification Challenges and Approaches. IFAC Proceedings Volumes, 42(10):1–20, 2009.
[20] Simo Särkkä. Bayesian Filtering and Smoothing. Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2013.
[21] Thomas B. Schön, Adrian Wills, and Brett Ninness. System identification of nonlinear state-space models. Automatica, 47(1):39 – 49, 2011.
[22] Michail D. Vrettas, Yuan Shen, and Dan Cornford. Derivations of Variational Gaussian Process Approximation Framework. Technical report, Aston University, March 2008.
[23] E. A. Wan and R. Van Der Merwe. The unscented Kalman filter for nonlinear estimation. In Proceedings of the IEEE Adaptive Systems for Signal Processing, Communications, and Control Symposium. IEEE, October 2000.

References

[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi: a software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, July 2018.
[2] D. Blei, A. Kucukelbir, and J. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859–877, 2017.
[3] Bob Carpenter, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A Probabilistic Programming Language. Journal of Statistical Software, 76(1), 2017.
[4] S. B. Chitralekha, J. Prakash, H. Raghavan, R. B. Gopaluni, and S. L. Shah.
Comparison of Expectation-Maximization based parameter estimation using Particle Filter, Unscented and Extended Kalman Filtering techniques.
IFAC Proceedings Volumes, 42(10):804–809, 2009.
[5] Jarrad Courts, Adrian Wills, and Thomas B. Schön. Variational Nonlinear State Estimation. arXiv, 2020.
[6] K. Furuta, M. Yamakita, and S. Kobayashi. Swing-up Control of Inverted Pendulum Using Pseudo-State Feedback. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 206(4):263–269, November 1992.
[7] Andreas Griewank and Andrea Walther. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. Society for Industrial and Applied Mathematics, jan 2008.
[8] Johannes Hendriks, Adrian Wills, Brett Ninness, and Johan Dahlin. Practical Bayesian System Identification using Hamiltonian Monte Carlo. arXiv, 2020.
[9] Bin Jia, Ming Xin, and Yang Cheng. High-degree cubature Kalman filter. Automatica, 49(2):510–518, feb 2013.
[10] Michael I. Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul. An Introduction to Variational Methods for Graphical Models. Machine Learning, 37(2):183–233, November 1999.
[11] Simon J. Julier and Jeffrey K. Uhlmann. New extension of the Kalman filter to nonlinear systems. In Ivan Kadar, editor, Signal Processing, Sensor Fusion, and Target Recognition VI. SPIE, July 1997.
[12] Juho Kokkala, Arno Solin, and Simo Särkkä. Sigma-Point Filtering and Smoothing Based Parameter Estimation in Nonlinear Dynamic Systems. Journal of Advances in Information Fusion, 2015.
[13] S. Kullback and R. A. Leibler. On Information and Sufficiency. The Annals of Mathematical Statistics, 22(1):79–86, March 1951.
[14] Fredrik Lindsten, Jouni Helske, and Matti Vihola. Graphical model inference: Sequential Monte Carlo meets deterministic approximations. In Advances in Neural Information Processing Systems 31, pages 8190–8200. 2018.
[15] L. Ljung. System Identification: Theory for the User. Prentice Hall information and system sciences series. Prentice Hall PTR, 1999.
[16] Lennart Ljung. Perspectives on system identification. Annual Reviews in Control, 34(1):1–12, April 2010.
[17] Christian A. Naesseth, Fredrik Lindsten, and David M. Blei. Markovian Score Climbing: Variational Inference with KL(p||q). In Advances in Neural Information Processing Systems 33. 2020.
[18] Radford M. Neal.
Handbook of Markov Chain Monte Carlo
, chapter MCMC using Hamiltonian dynamics, pages 113–162. Chapman and Hall/CRC, 2011.
[19] Brett Ninness. Some System Identification Challenges and Approaches. IFAC Proceedings Volumes, 42(10):1–20, 2009.
[20] Simo Särkkä. Bayesian Filtering and Smoothing. Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2013.
[21] Thomas B. Schön, Adrian Wills, and Brett Ninness. System identification of nonlinear state-space models. Automatica, 47(1):39 – 49, 2011.
[22] Michail D. Vrettas, Yuan Shen, and Dan Cornford. Derivations of Variational Gaussian Process Approximation Framework. Technical report, Aston University, March 2008.
[23] E. A. Wan and R. Van Der Merwe. The unscented Kalman filter for nonlinear estimation. In Proceedings of the IEEE Adaptive Systems for Signal Processing, Communications, and Control Symposium. IEEE, October 2000.

Appendix A Decomposition of $L (β)$

In this section, the decomposition of $L (β)$ required to express $L (β)$ as pairwise joints distributions, rather than as a full distribution is presented. This utilises the Markovian nature of state-space models in a fashion similar to related system identification approaches such as [21, 4, 20, 22, 5].

Firstly, $L (β)$ is expressed in the form of two integrals as

	$L (β)$	$= \int q_{β} (x, θ) log p (x, θ, y) d τ$
		$- \int q_{β} (x, θ) log q_{β} (x, θ) d τ .$

The first integral is decomposed according to


	$= \int q_{β} (x, θ) log p (x, y ∣ θ) d τ + \int q_{β} (x, θ) log p (θ) d τ$
	$= \int q_{β} (x, θ) log p (x_{1} ∣ θ) d τ$
	$+ T \sum k = 1 \int q_{β} (x, θ) log p (x_{k + 1}, y_{k} ∣ x_{k}, θ) d τ$
	$+ \int q_{β} (x, θ) log p (θ) d τ$
	$= \int q_{β} (x, θ) log p (x_{1}, θ) d τ$
	$+ T \sum k = 1 \int q_{β} (x, θ) log p (x_{k + 1}, y_{k} ∣ x_{k}, θ) d τ$
	$= \int q_{β} (x_{1}, θ) log p (x_{1}, θ) d [x_{1}, θ]$
	$+ T \sum k = 1 \int q_{β} (x_{k : k + 1}, θ) log p (x_{k + 1}, y_{k} ∣ x_{k}, θ) d [x_{k : k + 1}, θ]$
	$= I_{1} (β) + I_{23} (β) .$

While the second integral is decomposed according to

	$\int q_{β} (x, θ) log q_{β} (x, θ) d τ$
	$= \int q_{β} (x, θ) log q_{β} (x ∣ θ) d τ + \int q_{β} (x, θ) log q_{β} (θ) d τ$
	$= \int q_{β} (x, θ) log q_{β} (x_{1} ∣ θ) d τ$
	$+ T \sum k = 1 \int q_{β} (x_{x}, θ) log q_{β} (x_{k + 1} ∣ x_{k}, θ) d τ$
	$+ \int q_{β} (x, θ) log q_{β} (θ) d τ$
	$= \int q_{β} (x, θ) log q_{β} (x_{1}, θ) d τ$
	$+ T \sum k = 1 \int q_{β} (x_{x}, θ) log \frac{q_{β} (x_{k + 1}, x_{k}, θ)}{q_{β} (x_{k}, θ)} d τ .$

log \frac{q_{β} (x_{k + 1}, x_{k}, θ)}{q_{β} (x_{k}, θ)} = log q_{β} (x_{k + 1}, x_{k}, θ) - log q_{β} (x_{k}, θ),

this gives

	$\int q_{β} (x, θ) log q_{β} (x, θ) d τ$
	$= T \sum k = 1 \int q_{β} (x, θ) log q_{β} (x_{k + 1}, x_{k}, θ) d τ$
	$- T \sum k = 2 \int q_{β} (x, θ) log q_{β} (x_{k}, θ) d τ,$

and therefore,

	$\int q_{β} (x, θ) log q_{β} (x, θ) d τ$
	$= T \sum k = 1 \int q_{β} (x_{k + 1}, x_{k}, θ) log q_{β} (x_{k + 1}, x_{k}, θ) d [x_{k : k + 1}, θ]$
	$- T \sum k = 2 \int q_{β} (x_{k}, θ) log q_{β} (x_{k}, θ) d [x_{k}, θ]$
	$= I_{4} (β) .$

Together, this allows $L (β)$ to be calculated without requiring the full distribution, $q_{β} (x, θ)$ , but rather only each $q_{β} (θ, x_{k}, x_{k + 1})$ using

Appendix B Inverted Pendulum Model

The state vector used to model the Furata pendulum is $x = {[\begin{matrix} ϑ & α & ˙ ϑ & ˙ α \end{matrix}]}^{T},$ where $ϑ$ and $α$ are the arm and pendulum angles, respectively.

The continuous time dynamics are given by,

˙ x = {[\begin{matrix} ˙ ϑ & ˙ α & ¨ ϑ & ¨ α \end{matrix}]}^{T},

where

[\begin{matrix} ¨ ϑ ¨ α \end{matrix}] = M^{- 1} ([\begin{matrix} τ - D_{r} ˙ ϑ - D_{p} ˙ α - \frac{1}{2} m_{p} l_{p} g sin (α) \end{matrix}] - C),

and

	$C$	$= ⎡ ⎣ \begin{matrix} \frac{1}{2} m_{p} l_{p}^{2} sin (α) cos (α) ˙ ϑ ˙ α - \frac{1}{2} m_{p} l_{p} l_{r} sin (α) {˙ α}^{2} - \frac{1}{4} m_{p} l_{p}^{2} cos (α) sin (α) {˙ ϑ}^{2} \end{matrix} ⎤ ⎦,$
	$τ$	$= \frac{k_{m} (V_{m} - k_{m} ˙ ϑ)}{R_{m}},$
	$M$	$= [\begin{matrix} M_{11} & M_{12} M_{21} & M_{22} \end{matrix}],$

where

	$M_{11}$	$= J_{r} + m_{p} l_{r}^{2} + \frac{1}{4} (m_{p} l_{p}^{2} - m_{p} l_{p}^{2} {cos}^{2} (α)),$
	$M_{12}$	$= M_{21}$
		$= \frac{1}{2} m_{p} l_{p} l_{r} {cos}^{2} (α),$
	$M_{22}$	$= J_{p} + \frac{1}{4} m_{p} l_{p}^{2},$

and $m_{p}$ is the pendulum mass, $l_{r}$ , $l_{p}$ and the rod and pendulum lengths respectively, $J_{r}$ , $J_{p}$ are the rod and pendulum inertias, $R_{m}$ is the motor resistance, $k_{m}$ is the motor constant, $D_{r}$ , $D_{p}$ are damping coefficients for the rod and pendulum respectively, and $V_{m}$ is the applied motor voltage input.

The process model used consists of a two-step Euler discretisation of these continuous time dynamics over an $8 m s$ sampling time subsequently disturbed by noise $v_{k}$ .

Measurements are from encoders on the arm and pendulum angle and current measurement from the motor. The resultant measurement model is

y_{k}

= {[\begin{matrix} ϑ & α & \frac{V_{m} - k_{m} ˙ ϑ}{R_{m}} \end{matrix}]}^{T} + w_{k},

and it is assumed that

[\begin{matrix} v_{k} w_{k} \end{matrix}] \sim N (0, Π) .

Variational State and Parameter Estimation

Authors

Variational Nonlinear System Identification

Constructing a variational family for nonlinear state-space models

Self-learning eigenstates with a quantum processor

A Point Mass Proposal Method for Bayesian State-Space Model Fitting

On-line Bayesian parameter estimation in general non-linear state-space models: A tutorial and new results

The Mixture of Markov Jump Processes: Monte Carlo Method and the EM Estimation

Copula-like Variational Inference

1 Introduction

2 Variational Inference

3 Assumed Gaussian Distribution and Tractable Approximation

4 Numerical Examples

4.1 Stochastic Volatility

4.2 Inverted Pendulum

5 Conclusion

References

References

Appendix A Decomposition of $L (β)$

Appendix B Inverted Pendulum Model

Variational State and Parameter Estimation

Authors

Related Research

Variational Nonlinear System Identification

Constructing a variational family for nonlinear state-space models

Self-learning eigenstates with a quantum processor

A Point Mass Proposal Method for Bayesian State-Space Model Fitting

On-line Bayesian parameter estimation in general non-linear state-space models: A tutorial and new results

The Mixture of Markov Jump Processes: Monte Carlo Method and the EM Estimation

Copula-like Variational Inference

This week in AI

1 Introduction

2 Variational Inference

3 Assumed Gaussian Distribution and Tractable Approximation

4 Numerical Examples

4.1 Stochastic Volatility

4.2 Inverted Pendulum

5 Conclusion

References

References

Appendix A Decomposition of L(β)

Appendix B Inverted Pendulum Model

Appendix A Decomposition of $L (β)$