{'fraction_non_alphanumeric': 0.03437383410023629, 'fraction_numerical': 0.008978982713592838, 'mean_word_length': 4.59738270917444, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Practical Bayesian System Identification using Hamiltonian Monte Carlo

This paper considers Bayesian parameter estimation of dynamic systems using a Markov Chain Monte Carlo (MCMC) approach. The Metroplis-Hastings (MH) algorithm is employed, and the main contribution of the paper is to examine and illustrate the efficacy of a particular proposal density based on energy preserving Hamiltonian dynamics, which results in what is known in the statistics literature as ``Hamiltonian Monte--Carlo'' (HMC). The very significant utility of this approach is that, as will be illustrated, it greatly reduces (almost to the point of elimination) the typically very high correlation in the Metropolis--Hastings chain which has been observed by several authors to restrict the application of the MH approach to only very low dimension model structures. The paper illustrates how the HMC approach may be applied to both significant dimension linear and nonlinear model structures, even when the system order is unknown, and using both simulated and real data.


Introduction
In many important applications, decisions must be made in the presence of noise corrupted data and dynamic system model uncertainty. The impact of these decisions ranges from trivial to catastrophic such as the instability of a control loop. Consequently, the ability to manage uncertainty is vital, and this paper is directed at this problem by illustrating how when dynamic system models are determined from noise corrupted data, the associated estimated system uncertainty can be quantified.
Furthermore, we illustrate how the uncertainty in functions of the estimated dynamic system, such as control system performance measures, e.g. stability margin for a putative control law, can also be quantified.
The topic of estimating dynamic system models from observed data has been studied in earnest for half a century, with the overwhelming majority of solutions being highly effective parametric methods based on optimisation techniques such as least squares (LS) and maximum likelihood (ML) and generalisation such as the prediction error (PE) approach.
These approaches have their origins in the statistics literature. Except for very simple cases, such as linearly parameterised model structures, quantifying the error in the associated parameter estimates that arise from stochastic data corruption relies on the employment of a central limit theorem (CLT). Unfortunately, this gives an error quantification that is approximate, since it relies on assuming the asymptotic in data length N → ∞ CLT quantification is approximately true for finite data length N . Furthermore, it is often very challenging to apply this CLT-based approach to nonlinear dynamic system estimation. This paper is motivated by these difficulties, particularly that of accurate error quantification for nonlinear systems and/or small data length N , and examines the utility of a solution based on a Bayesian methodology. The rationale for this approach is that error bounds based on Bayesian posterior densities are non-asymptotic in data length N , and the ability to inject quite sophisticated prior knowledge can significantly improve the accuracy of the obtained estimates.
Despite these attractive features, the Bayesian approach is not without its challenges. Foremost is that it is an approach that is not based on numerical optimisation methods that underpin LS, ML and PE methods, but instead requires numerical integration methods. In very special cases such as the linearly parametrised linear Gaussian model with so-called "conjugate" priors, the associated integration can be provided in closed form, and in more general cases of very low model dimension the required integrals can be computed numerically using, for example, quadrature methods.
However, in this latter case the number of function evaluations required grows as a parameter-dimension exponential of the number of integration "bins", and hence grows dramatically in both. Furthermore, the range of integration for each parameter to be estimated must be specified, which for the problems considered here is impossible since the parameters themselves are unknown a priori. This renders this numerical integration approach inapplicable except in very special very low parameter dimension situations.
An alternative numerical integration approach studied in this paper is based on constructing a random number generator for which the samples are distributed according to the Bayesian posterior distribution. The utility of this random number generator is that sample-averages of the random number generator realisations can be used to approximate, with arbitrary accuracy for sufficiently large number of samples, the required integrals. For example, it can provide estimates of the conditional mean, which is known to be the minimum mean-squared estimate.
Constructing such a random number generator has received significant research attention and among the possible solutions is the so-called Markov-Chain approach. As the name suggests, this approach constructs a Markov-Chain whose stationary distribution coincides with desired posterior density target distribution. Ensuring that this Markov-Chain does indeed converge to the desired stationary distribution can be achieved using the decades old and ingenious Metropolis-Hastings (MH) algorithm which dates back to the 1950's were it was developed to compute problems in physical chemistry [10,9].
Despite achieving this remarkable outcome, the MH method has a significant drawback in that it requires a certain user-supplied "proposal density". This proposal must be carefully selected to deliver Markov chain realisations that are sufficiently uncorrelated that the approach delivers sufficiently accurate computation of desired integrals without requiring an unreasonably large amount of realisations being generated. Historically, this practical aspect has prohibited the application of the MH algorithm to dynamic system estimation of all but the simplest model structures.
This paper addresses this weakness by employing the so-called "Hamiltonian Monte Carlo" approach, which uses a particular proposal density that, as will be illustrated here, renders the MH-algorithm far less sensitive to user tuning and scales far more successfully to higher dimensions.
Hamiltonian Monte Carlo (HMC) is a relatively recent development originating in the physics community [19] and then taken up by the statistics community [6,12] and now is attracting interest in machine learning [2].
However, despite its great potential, it appears to be largely unrepresented in the control community. This paper is directed at filling this gap by introducing the HMC approach in applications ranging from classical linear transfer function models through to quite general nonlinear state-space models. Furthermore estimation of models of a-priori unknown model order and estimation together with error quantification of functions of the parameter estimates are also illustrated.

Problem Formulation
The problem addressed here is that of using measured data from a dynamic system consisting of inputs u t and outputs y t for t = 1, . . . , N to determine a model for the system.
Typically, such a model will have several unknown parameters which we collect in a vector θ ∈ R n θ and will typically be a function of past inputs and outputs plus also potentially latent variables, such as system states x t .
A very general mathematical description of this model could be written as y t+1 = f (θ, y 1:t , u 1:t+1 , x 1:t , e 1:t+1 ), where we use the shorthand y 1:t to denote the set {y 1 , . . . , y t } etc. , x t is a possible state vector, and e t is an unobserved stochastic sequence where all realisations are independent and identically distributed (i.i.d).
In this paper, the function f can be quite general, but in simple cases could take the form of, for example, a transfer function or linear state-space model.
The aim of this work is to evaluate expectation integrals with respect to a posterior distribution p(θ | y 1:t ) of the form where ϕ : θ → C n is a general function of the parameters, such as the conditional mean, Bode frequency response, or other almost arbitrary functions such as achieved phase margin for a given controller.
Throughout this paper, while the vector θ ∈ R n θ denotes the model parameters the vector η ∈ R nη with n η ≥ n θ will be used to denote the full vector of elements we wish to draw samples from using MCMC. It is important to make this distinction as for some problems η could also include values additional to θ which are latent variables, such as the system state history x 1:N or quantification of any noise sequences such as their variances.

A Markov Chain Monte Carlo Approach
The computational approach studied here relies on the Law of Large Numbers to provide estimates of (2) via so-called Monte-Carlo integration where samples are used to compute a sample average of the integral I in (2) according to Under rather general assumptions on p(η | y 1:N ) and ϕ(·), it can be shown that where the convergence is with probability 1. The "Monte Carlo" epithet is in reference to the famous casino, and a nod to the random aspects of the approach.
Note that ϕ(·) may be chosen as the indicator function which is equivalent to I M being an estimate of the posterior density p(η | y 1:t ), or marginals of it, via a sample histogram of the realisations of all of {η i } or the component of it for which a marginal is of interest.
Importantly, the rate of convergence in (5) depends only on the correlation between the realisations {ϕ(η i )} and not on the dimension n η of η, thus removing the exponential computational load problem associated with quadrature and other deterministic methods mentioned in the introduction.
Furthermore, bounds for integration limits need not be specified a-priori since the realisations naturally arise mainly in regions that significantly affect the evaluation of I and rarely in regions that do not.
However, the ability to use this approach is predicated on being able to draw the samples (3). A remarkably simple and elegant approach to solve this problem is to employ the so-called "Metropolis-Hastings" (MH) algorithm.

The Metropolis-Hastings Algorithm
The Metropolis-Hastings algorithm simulates a Markov chain using two steps at each iteration k. First, a can-
Second, a desired "target" density π(η) from which it is desired to draw realisations is specified, such as π(η) = p(η | y 1:N ) and then the candidate parameter is accepted with probability given by .
The final step is to set η k ← η when the candidate is accepted and η k ← η k−1 when it is rejected. The full procedure is presented in Algorithm 1.
In essence this approach is such that an arbitrary Markov chain that the user can easily simulate, as specified by the "proposal distribution" q(η k | η k−1 ) is then modified to a Markov chain whose stationary density is the required target π(η).

Sufficient Invariance Condition -Detailed Balance
It is important to understand the fundamentals of the design of the MH algorithm. First, observe that it generates a time-homogeneous Markov chain with transition kernel K(η k | η k−1 ) given by where, to avoid unnecessary technicalities that may obscure intuition, we neglect the very rare event that drawing from the proposal q(· | η k−1 ) for η k would deliver exactly the same value η k−1 as in the previous iteration.
Therefore, the probability density π k for the k'th iteration of the chain, is related to the one π k−1 at the previous iteration according to the law of total probability as and therefore, if the chain is to converge to deliver realisations from a stationary distribution π then the latter must satisfy the equation (12) in which case π is referred to as an "invariant" distribution of the chain.
The MH algorithm is specifically designed so that the kernel K(· | ·) it implements has the target π as an invariant distribution by ensuring the chain it simulates satisfies a condition which is sufficient (but not necessary) for invariance, which is termed the "detailed balance" condition: To see the sufficiency of this condition for invariance, note that it implies To see that the very simple design of the MH algorithm ensures detailed balance, we compute Comparing (15) and (16) and noting that the min{·, ·} operation is symmetric establishes that the "detailed balance" condition (13) and hence the invariance condition (12) is ensured.
For more extensive detail on this and related MH convergence analysis the interested reader is referred to [21,13,18].

Constructing proposal distributions
The proposal distribution q(· | ·) is the only tuning aspect of the MH algorithm and determines its performance in terms of how quickly the accuracy of the approximation (4) improves with increasing number of realisations L.
One way of judging this rate is to employ the Central Limit Theorem (CLT) for the estimate (4) which is given by [18]: where where the integrated auto-correlation time (IACT) is introduced and infers the convergence rate approximation On the one hand, this result emphasises the great attraction of this MCMC method for computing integrals since the 1/L convergence rate in (19) is independent of the dimension of η which determines the dimension of any integrals involved.
On the other hand, the dimension of η is still an issue since it is typically observed that as it grows, the complexity of the target π(η) grows, and it becomes increasingly difficult to design a proposal q(· | ·) such that the correlation corr(ϕ(η 0 ), ϕ(η k )) decreases quickly with increasing k.
For example, for a proposal q(· | ·) with large steps designed to explore π(·) well, one or more elements of η are often moved well away from the support of π(·) and draw α(· | ·) very low so that almost all steps are rejected, and hence the chain stays in the same place for very long periods and is highly correlated.
Vice versa, if only tiny steps are proposed by q(· | ·) in order to ensure many acceptances, the chain wanders very slowly, and again in a highly correlated manner akin to a random walk.
Both of these generate a very high, perhaps infinite value for σ 2 ϕ in (19). The literature on MCMC therefore concentrates strongly on the so-called "mixing" of the chain with "strong mixing" meaning a quick decrease in corr(ϕ(η 0 ), ϕ(η k )) with increasing k.
One commonly used approach to address this issue is to use the proposal choice (7) with the choices where Σ post denotes an estimate of the target distribution covariance. The rationale here is that the support of a Gaussian with these choice for mean and covariance approximates the support of the target π.
Typically, this is obtained by starting with a covariance Σ(η k−1 ) = 2 I for some choice of , and then after an initial pilot run {η} k has been obtained, computing Σ π via a sample average estimate from this run.
Often, this involves trial and error choice of to achieve sufficient mixing of the pilot run, which in many cases of non-trivial dimension is not possible.
Another approach for constructing a proposal is to borrow ideas from the non-linear optimisation literature and employ a quasi-Newton optimisation step as a proposal by again using the Gaussian proposal choice (7) with Here, we introduce the gradient of the log-posterior, and the negative inverse Hessian of the log-posterior This strategy takes into account the gradient and curvature of the target distribution and therefore introduces a mode-seeking behaviour, which guides the Markov chain to areas of higher likelihood.
We refer to this type of proposal as the manifold Metropolis adjusted Langevin algorithm (mMALA) as it is known in the statistics literature, see [6] where it has been successfully applied to, for example, econometric problem with up to 20 states and 50+ parameters) [1].
While these proposal choices (20) and (21) are well accepted as improving the performance of the MH method, they fall short of achieving the optimum of delivering a sequence of realisations {η k } that are both uncorrelated, and therefore also involve few rejections. The Hamiltonian Monte Carlo Method (HMC) is designed to achieve these latter goals.

Hamiltonian Monte Carlo
The essential idea of Hamiltonian Monte Carlo (HMC) is to augment the parameter space by adding an additional vector ρ where and then employ Metropolis Hastings to generate samples from a new augmented target density for some user chosen density p(ρ). In this case, the marginal of π(η, ρ) with respect to ρ is the target π(η) of interest by the law of total probability: In HMC, the "augmented" density π is formulated as where with and H(η, ρ) is interpreted as the Hamiltonian of a dynamic system in η(t), ρ(t) that evolves according to Hamilton's equations: so that η can be imagined as a position vector and ρ a momentum vector. In this case V(η) can be considered the potential energy and K(ρ) the kinetic energy of a particle in motion governed by (28).
Set (η , ρ ) ← (ηL, −ρL) as the proposed sample. 5: Sample u uniformly over [0, 1]. 6: if u ≤ min{1, α(η , ρ , η k−1 , ρ k−1 )} given by (31). then 7: end if 11: end for Importantly, Hamiltonian dynamics are energy preserving as can be simply established via Therefore, if we construct a MH algorithm for which a new realisation (η , ρ ) is governed by the dynamics (28) then, if in addition the proposal γ(· | ·) is symmetric the acceptance probability for this new realisation is thus realising the ambition of a high acceptance rate. Realising the ambition of uncorrelatedness between (η , ρ ) and (η, ρ) can be achieved by drawing a random proposal ρ ∼ p(ρ) independant of any previous iterations and then simulating the Hamiltonian system over some user chosen time period t ∈ [0, L] to deliver (η , ρ ) independent of (η, ρ).
This is the HMC algorithm, which for the choice of inde- is formally defined in Algorithm 2 with M referred to as a "mass matrix": There are some important subtleties that require attention. First, the reader may note that the simulation step is a deterministic transformation Φ(η, ρ) → (η , ρ ) and
hence a formal change of variables calculation should be applied when computing the acceptance probability (31). This can be accomplished using the well know change of variables formula [14] which would involve a multiplicative factor in (31) However, an important property of Hamiltonian dynamics is that they are volume preserving, so that this determinant is always equal to one and can be neglected [12].
Furthermore, Hamilton's equations (28) are time reversible in the sense that if they govern the movement of a particle for L seconds, if the momentum of that particle is reversed in sign, then after another L seconds it will return to where it originally was. This explains the negative sign on ρ L in step 4 of Algorithm 2 since it ensures the detailed balance condition (13) holds and that π is an invariant density for the resulting chain [19].
Finally, the reader may also note that between lines 6 and 10 of Algorithm 2 the acceptance ratio is computed and used to accept or reject, even though as established in (32) this acceptance ratio should always be one, and hence no rejections.
The former accounts for the fact that numerical simulation of the Hamiltonian system will necessarily involve some errors, hence exact equality of the acceptance probability to unity cannot be expected.
Simulating the Hamiltonian system from t = 0 to t = L involves computing the integrals ∂η dt (34) Fortunately, there exists a class of so-called "symplectic" integrators designed specifically for the simulation of Hamiltonian systems and to be as close to volume preserving as possible. One deceptively simple one used in HMC applications is the "leap frog" method specified in Algorithm 3 again for the kinetic energy choice (33) which is used at line 3 of the HMC Algorithm 2.

HMC Tuning Choices
Unfortunately, there are some user defined choices that need to be made to implement HMC. Namely, the Hamiltonian simulation duration L, the symplectic integration step size and the kinetic energy mass matrix M which is the covariance of the Normal distribution for new momentum proposals.
While L and can be statically defined [2], a state-ofthe-art implementations uses a so called "No-U-Turn" criteria [7,4] to dynamically adjust the integration time, and the integration step size is dynamically tuned to achieve the desired acceptance rate.
Additionally, there is the choice of mass matrix M which can be chosen by computing an empirical estimate of the target covariance using a warm-up phase of the Markov chain. This choice of mass matrix decorrelates the parameters and gives more uniform energy levels.
However, unless the target distribution is exactly Gaussian no choice of global rotation and scaling will yield completely uniform level sets.

Example
To illustrate the performance advantages of the HMC proposal in a MH algorithm, we consider an example in a two-dimensional parameter space η = {η 1 , η 2 } where the target distribution π(η) has a doughnut-shaped support.
This may seem contrived, however, it is recognised that many common targets in higher dimensional spaces tend to occupy a narrow shell [2].
The result for a random walk proposal with the MH algorithm are shown in Figure 1, where clearly little movement over the first 10 iterations has occurred and most moves have been rejected.
By way of contrast, the moves for the HMC method over the first 10 iterations are shown in Figure 2 and show all moves accepted with much better exploration of this doughnut support target.
A short video showing the evolution of this HMC Markov chain is available at https://youtu.be/cyqh4mUB3bc.

Convergence
To this point, it has only been established that the HMC method implements a Markov chain for which the joint target π(η, ρ) is an invariant density of the chain with marginal over ρ satisfying (24). Being an invariant density it is then possible that the distribution of realisations will converge to π, but to this point there is no guarantee this will actually happen, and indeed it has not been established that π is unique as an invariant density.
Taking this latter question of uniqueness first, and to streamline notation set and define Γ {λ ∈ R n λ : q(λ) > 0} (37) where q(λ ) = q(η , ρ ) is the proposal density achieved by the HMC approach of Algorithm 2 of randomly drawing ρ ∈ N (0, M) and then simulating the Hamiltonian dynamics to deliver a proposed realisation η .
It is clear that the HMC method can only effectively sample the entire support Λ of the target π if Λ ⊆ Γ.
In fact, if this holds π is the unique invariant density.
Lemma 5.1 Suppose that Λ is connected, that (38) holds and that π(λ) is bounded over λ ∈ Λ. Then π(λ) is the unique invariant density of the Markov chain realised by Algorithm 2.
PROOF. The transition kernel K(ξ | λ) for the Markov chain realised by Algorithm 2 is given by Therefore, under the assumption (38) it holds that Define the distribution Π(A) associated with the invariant density π(λ) according to and define the conditional distribution P (A | λ) similarly: Then via (40) This establishes that if the chain converges in a distributional sense, then the target density π and its associated distribution is the only possibility. That this convergence actually occurs can also be established. Finally, we note that while this distributional convergence is a goal in designing Algorithm 2, ultimately its utility is in estimating expected values of function with respect to the target density via sample averages as per (4), (5). This required convergence is established in the following Theorem. R n λ and that ϕ : with probability one.
PROOF. It has already been established in the proof of Lemma 5.2 that the Markov chain realised by Algorithm 2 is Harris recurrent. The result then follows by Theorem 17.1.7 of [11]. 2 Note how central the assumption (38) is to all these convergence results. One of the drawbacks of the HMC method is that it is usually not straightforward to establish that (38) holds in practice since it will critically depend on the Hamiltonian simulation duration L.
If it is too short, then (38) may not hold. If too long then an uneccessary computional burden ensues. This has led some authors to suggest some specifics around the simulation step, such as randomising the simulation duration L [3].

HMC for Bayesian Dynamic System Identification
We now discuss the application of the HMC method to the dynamic system identification problem formulated in section 2 at the beginning of this paper. This involves using HMC with target π(η) = p(η | y 1:N ) the posterior density of the vector η given the observed data y 1:N .
Instrumental to achieving this is the use of Bayes' rule to decompose the posterior as where p(y 1:N | η) is the likelihood of the observed data for a given η, p(η) is a prior for η, and p(y 1:N ) is an ηindependent constant that ensures the area under the posterior density p(η | y 1:N ) is one.
In turn, by repeated applications of Bayes' rule the likelihood can be decomposed as The computation of the "one step ahead prediction density" p(y t | y 1:t−1 , η) then depends on the model structure chosen.
For example, for the general linear transfer function model structure where A(q, θ), · · · , D(q, θ) are polynomials in the shift operator q and parametrised by the vector θ, while {e t } is an i.i.d. sequence of zero mean random variables with variance σ 2 e it is possible to rewrite it as If C(q, θ) and D(q, θ) are constrained to be monic (zero'th order constant term in the polynomials fixed at 1) then [1 − H −1 ]y t is purely a function of y 1:t−1 with no dependence on y t and hence with the definition of the one step ahead predictor then if e t ∼ p e (·) the likelihood becomes (54) Similarly, in the case of linear, possibly time varying state space model structure where x t ∈ R nx is a state vector and A t (θ) · · · D t (θ) are matrices with dimensions conformal to the above and the dimensions of the possibly multivariable (in this case) y t and u t , then if the i.i.d. random variables ν t and e t in (55) are Gaussian distributed then the likelihood is again given by (54) with y t|t−1 computed using the standard Kalman filter recursions and η T = [θ T , vec{Σ} T ].
In the linear but non-Gaussian case, or the non-linear state space case the one-step ahead predictor density can still be formulated but the associated multi-dimensional integral calculation is generally intractable.
A strategy to deal with this problem is to include the full state history Depending on the model structure, these two component may well be tractable to compute. For example, for the possible nonlinear state space structure with w t ∼ p w (·), e t ∼ p e (·) (60) then it is straightforward to evaluate: Note that vice-versa, if smoothed state estimates are also of interest, then we can marginalise over the model parameters to compute expectations with respect to the smoothing density p(x 1:N | y 1:N ).

Choice of Prior
The choice of prior p(η) in (49) can be of great benefit in achieving very effective dynamic system estimation solutions. Indeed, it is one of the great benefits of the Bayesian approach that it provides a principled approach for injecting prior knowledge, to then be combined with evidence from the data via the likelihood p(y 1:N | η) to deliver the final posterior p(η | y 1:N ).
For example, priors can be used to promote smoothness, sparsity, stability and other properties of the estimated model and are closely connected to the technique of "regularisation" in non-Bayesian settings [16,17].
For example, it is very well known that Gaussian prior choice is closely related to imposing a 2-norm η 2 2 regularisation penalty in an optimisation based (e.g. Maximum Likehood) estimation approach. Similarly with L denoting the Laplacian density, the prior choice implies a 1-norm η 2 1 regularisation penalty, and is commonly used to induce sparseness in the estimated η.
It is also possible, and indeed common to formulate a prior in a hierarchical manner where ϑ is referred to as a "hyperparameter", which has its own prior attached to it as indicated above.
An example of this is the "horseshoe prior" [5,17,15] p(η i ) ∼ N (0, β 2 i τ 2 ), β i ∼ C + (0, 1), τ ∼ C + (0, 1), where C + (0, 1) denotes a zero-mean half-Cauchy (over the positive reals) with scale 1. Here, each β i is an independent shrinkage parameter relating to each η i and τ is a global shrinkage parameter. This prior is again designed to promote sparseness in the estimated vector η in that elements will be zero or close to it unless there is evidence from the data y 1:N otherwise. Although more complicated than the Laplace prior (64) also designed to produce sparseness, it has been observed to increase the accuracy of the resulting estimate [5].

Model Averaging
The parameter posterior for hierarchical Bayesian models is still computed using Bayes' rule. However, we are now required to marginalise the posterior over all the possible values of ϑ supported by the data, which corresponds to This is a type of model averaging where all possible choices of ϑ are averaged over and therefore the influence of ϑ is removed from the model. As a result, Bayesian methods are not prone to overfitting as no single value of ϑ is selected.
Although exact marginalisation is intractable for many models of interest, it can easily be performed in the case where we have samples from the joint posterior by discarding the ϑ k values. Sampling from the joint posterior could be performed using the MCMC methods, such as MH or HMC, discussed earlier

Numerical Examples
In this section, we illustrate the use of HMC on a range of examples from system identification. These simulations demonstrate that HMC provides a practical means for sampling from the required posterior distribution; this includes both classical model structures such as ARX and output-error models, and more challenging nonlinear state-space models. Depending on what is appropriate the estimates are assessed by; comparison to the true parameters, comparing the true and estimated frequency response, and/or using the model fit metric [8] given by where y t and y t denote the one-step-head predictor and the observation at time t, respectively.
The implementations are made in MATLAB and Python using the PyStan [20] interface to Stan [4]. The complete source code for these examples is provided on GitHub at https://github.com/jnh277/hmc_lin_sys.

Autoregressive exogenous
This section provides an example of identifying an ARX model from data. It compares the ability of MH, mMALA, and HMC to provide uncorrelated samples from the posterior density η ∼ p(η|y 1:T ) using the model given by where e t ∼ N (0, σ 2 e ) denotes zero-mean additive Gaussian noise. Here, y 1:T and u 1:T denotes the output and input, respectively. Finally, the parameters of the model to identify using input-output data are The ARX model is simulated for T = 1000 data points with order n a = n b = 2 and parameters {−1.5, 0.7, 0, 1, 0.5, 1} using a random binary input. Moreover, we assume that the true model order is known to simplify the inference. For inference, 667 data points are used; with the remainder reserved for validation. The arx command in MATLAB is used to obtain the maximum likelihood estimate.
A Bayesian estimate is obtained using HMC, MH, and mMALA with an L2 (Gaussian) prior on the parameters. Figure 3 shows the trace plots and automatic correlation factor for a 1 as well as the posterior marginal distributions of a 1 and b 1 . The trace is a plot of the state of the Markov chain at each iteration. Note that the chain from MH and mMALA exhibit a much larger auto-correlation than the chain generated by HMC. This is also apparent as the ACF falls off quicker for HMC compared with the other two approaches. Finally, calculating the model fit for each estimate gives 96.09 for Matlab's arx command, 96.08 using HMC, 96.04 using mMALA, and 94.03 using MH.
Recall the discussion in the previous section, a large persistent autocorrelation results in poor mixing and poor performance in the MCMC method. Even on this relatively low order example, the difference in the autocorrelation factor is a strong indicator of the greater efficiency of HMC sampling from the target distribution. This difference becomes even more noticeable as the dimension of the parameter space increases. Even for high dimensional parameter spaces, HMC using the No-U-Turn sampling algorithm is typically able to efficiently provide samples with minimal user tuning. Having made this comparison, the remainder of the examples in this section focus on HMC and details relevant to its use for the model type considered.

Auto-regressive exogenous, continued.
After the previous sanity check, we proceed with a slightly more interesting example where the model orders are unknown. In this setting, the standard approach to maximum likelihood inference would be to include a regularisation term or to use cross-validation to determine a suitable model order. In the Bayesian setting, the standard approach is to use a sparseness prior which regularises the model and turns off unnecessary parameters in an over-parameterised model. Using HMC, this example compares the effect of L1 (exponential) prior, L2 (Gaussian) prior, and the horseshoe prior.
In this example, we use n a = n b = 10 as a guessed maximum model order and compute the maximum likelihood estimate using Matlab's arx function with the popular tuned-correlated (TC) kernel employed to regularise the results. The different impacts that each prior has on the posterior distribution of the parameters are evident in Figure 4. Shown are the marginal distributions for the parameters that are not part of the true model, i.e. for k > 2, which indicates that the horseshoe prior is more effective at reducing the value of these parameters than the L1 or L2 priors. Also shown are the maximum likelihood estimates. As well as promoting sparseness, it would be good if choosing the horseshoe prior still yielded accurate estimates of the parameters. The model fit given using ML is 89.3, and for HMC using L1, L2, and horseshoe priors are respectively 89.2, 89.1, and 89.4. Additionally, the Nyquist diagram of the HMC samples with a horseshoe prior is shown in Figure 5 and compared to the true systems and the ML estimate.

Output error
This section demonstrates the Bayesian estimation of an output error (OE) model using HMC. Consider the model where the measurements are corrupted by additive Gaussian noise e t ∼ N (0, σ e ). The parameters of the model to identify using input-output data With a guessed maximum model order of n b = n f = 10, a Bayesian estimate of the model is given using HMC with a horseshoe prior, achieving a model fit of 99.42, and a maximum likelihood estimate using Matlab's arx command with regularisation achieving a model fit of 99.45. The corresponding Nyquist diagrams are shown in Figure 6.

Identification with measurement outliers
So far, the examples shown have only considered additive Gaussian noise. This assumption, however, is problematic if the measurements include outliers.
The assumption of additive Gaussian noise is not required for Bayesian inference using HMC, and a large variety of noise models can be easily implemented. This section revisits the example given in Section 7.2 and considers noise that has been simulated from a Student's T distribution, e t ∼ T ν, 0, σ 2 e . For low values of ν (the degrees of freedom), the Student's T distribution has fatter tails than a Gaussian distribution and so the measurements can have significant outliers.
In this example, we use n a = n b = 10 as a guessed model order and compare the results from Matlab's arx command with TC regularisation against estimating the model using HMC and a Student's T distribution for the noise. Both the standard deviation and the degrees of freedom of the noise are estimated. Since for ν > 10, the Student's T distribution becomes practically Gaussian, hence using this model is not restrictive. For the parameters, a horseshoe prior is used. This gives the parameters to be estimated using HMC as η = {a 1:na , b 0:n b , σ e , ν}, where ν is constrained within [1, ∞] and is given a Gamma prior, G(2, 0.1).
The Nyquist diagrams corresponding to the estimated models are shown in Figure 7. Despite, Matlab's arx method using a Huber style cost to improve performance when outliers are present it's performance is degraded for this data. In contrast, HMC provides good estimates despite the measurement outliers.
where m p is the pendulum mass, L r , L p and the rod and pendulum lengths, J r , J p are the rod and pendulum inertias, R m and k m are the motor resistance and constant, D p and D r are the pendulum and arm damping constants, respectively.
Measurements of the two angles, θ and α, and the motor current I m made at 8ms intervals are available. These measurements are modelled according to Given that the parameters m p , L p and L r can easily be measured accurately, we focus on estimating the parameters J r , J p , K m , R m , D p , and D r using the input-output data {V m,t , y t } T t=1 . To do this we choose as latent variables x(t) = θ(t) α(t)θ(t)α(t) T . This allows us to write a state space model of the form (55). Additionally, covariance between the process noise, e t , and measurement noise w t was considered giving the joint likelihood as where f and h are given by equations (71) and (72), respectively, combined with an ODE solver. In this case, the Runge-Kutta 4 th order method was used.
HMC was used to infer the joint posterior of x 1:T , J r , J p , K m , R m , D p , D r , and Π from an experimentally collected data set containing 375 input-output pairs. The model parameters were constrained to be positive with a horseshoe prior and an LKJ prior on the covariance matrix. The model parameters were initialised with a value of 1 and the latent variables were initialised at the measurement values, where appropriate. The marginal posterior distributions for the parameters are shown in Figure 8. From this we can see that the most likely value for D p is close to zero as joints use magnetic bearings and whereas D r has a non-zero value introduced by cabling slightly interfering with the smooth movement of the arm rotation.
The estimates of the arm angles and angular velocities are compared to the measured values (and gradients of the measured values) in Figure 9. As HMC provides the full joint posterior of the states and the parameters there is far more information we can look at. Figure 10 gives examples of marginal and joint distributions between parameter-state, parameter-measurement, and  . Posterior distributions of the parameters of the nonlinear inverted pendulum system: km,α at t = 1.8s, and the predicted measurement of Im at t = 1.8s. The joint distribution between each of pair is represented by the scatter plots. Also shown is the measured value for Im indicated by a dashed line.
state-measurement pairs. From this, for instance, we can see that there is correlation between the value for k m and the predicted measurements of I m .

Conclusion
This paper has presented Bayesian inference using Hamiltonian Monte Carlo within the context of system identification. Hamiltonian Monte Carlo is presented as a specific choice of the proposal for Metropolis-Hastings and its advantages over other proposals, such as mMALA, have been identified. In particular, it allows for more efficient sampling from high dimensional parameter spaces with minimal user tuning. This allows the application of HMC for Bayesian system identification to a broad range of model types. These model types include classical models such as ARX and OE, through to nonlinear state-space models. The ability to include the latent variables of a state-space as parameters to be sampled allows joint parameter estimation and state smoothing to be performed using HMC.