Green, PL

Maskell, S

It is now well established that, through the use of sequential Monte Carlo methods, it is possible to track the time-varying state of mechanical systems 'online', using a continous stream of measurements. The best known of these algorithms is the particle filter-a numerical algorithm that can be applied to a large variety of nonlinear problems and which, in recent years, has been used to aid the condition monitoring of many mechanical systems. In this paper, a Sequential Monte Carlo method is used to estimate the parameters of a model from a continuous stream of measurements, with the aim of establishing how one's parameter estimates converge as more data is analysed. Crucially, for reasons described in this paper, this is a situation where a particle filter is unsuitable. The issue is instead resolved using a variant of a Sequential Monte Carlo sampler. It is shown how the algorithm can be used to identify the parameters of a model from large data sets and, within the context of structural dynamics, it is compared with the performance of a similar, Markov chain Monte Carlo method

University of Liverpool Repository

Parameter estimation from big data using a sequential montecarlo samplerP.L. Green, S.MaskellInstitute for Risk and Uncertainty, University of Liverpool, Liverpool, UK, L69 3GQe-mail: p.l.green@liverpool.ac.ukAbstractIt is now well established that, through the use of sequential Monte Carlo methods, it is possible to trackthe time-varying state of mechanical systems ‘online’, using a continous stream of measurements. The bestknown of these algorithms is the particle filter - a numerical algorithm that can be applied to a large varietyof nonlinear problems and which, in recent years, has been used to aid the condition monitoring of manymechanical systems. In this paper, a Sequential Monte Carlo method is used to estimate the parameters ofa model from a continuous stream of measurements, with the aim of establishing how one’s parameter es-timates converge as more data is analysed. Crucially, for reasons described in this paper, this is a situationwhere a particle filter is unsuitable. The issue is instead resolved using a variant of a Sequential Monte Carlosampler. It is shown how the algorithm can be used to identify the parameters of a model from large datasets and, within the context of structural dynamics, it is compared with the performance of a similar, Markovchain Monte Carlo method.1 IntroductionConsider a scenario where one wishes to infer probabilistic estimates of a model’s parameters, θ, from a setof measurement data, z1:n = {z1, ...,zn}. In this context zt represents a vector of measurements taken attime t and so z1:n = {z1, ...,zn} represents the set of all measurements taken up to time n. These could, forexample, represent time histories of acceleration measurements taken from a dynamic structure. Adoptinga Bayesian framework, one can realise probabilistic estimates of θ, conditional on z1:n, from the posteriorparameter distribution:p(θ |z1:n) ∝ p(z1:n|θ)p(θ). (1)At this stage one will often generate samples from p(θ |z1:n) using, for example, Markov chain Monte Carlo(MCMC) methods. When n is large however, this can be very expensive, and so a subset of this training datamust be used. How this subset should be selected is often unclear, and it can be difficult to establish whena sufficient amount of information has been extracted from the ‘full’ dataset. An information-theoretic solu-tion to this issue was proposed in [1] which, while useful in some scenarios, does rely on one approachingthe problem with a reasonable estimate of the most probable vector of model parameters. MCMC methodsdesigned to address this scenario were proposed in [2, 3]. These allow one to track how one’s estimates of θconverge as an increasingly large set of data is analysed. Unfortunately, the computational cost of applyingthese methods in situations where large amounts of data are required can still be prohibitive. This is alsoan issue which cannot be tackled using a particle filter, as this would lead to the influence on one’s initialmeasurements decreasing as more data is analysed (see [4] for more details).The current paper proposes a method which, the authors believe, can tackle this problem at reduced compu-tational cost and in a way which is very well suited to parallel processing. It therefore tackles ‘Big Data’ inthe sense that it aims to facilitate the analysis of data sets which are too large for current methodologies. Theproposed algorithm is essentially a variant of a Sequential Monte Carlo (SMC) sampler (which, somewhatconfusingly, is part of the family of SMC methods).2 Prerequisites2.1 Importance samplingImportance sampling forms the foundation of the method proposed here and, as such, it is given a brief in-troduction in this section. For more information, nice introductions to importance sampling can be found inthe books [5, 6] as well as numerous tutorials.Here, π(θ) is used to represent the target distribution which, in the context of this paper, is the posteriorparameter distribution in equation (1). Consider the situation where one wishes to estimate the expectedvalue of some function, f(θ), (where θ ∼ π(θ)). With importance sampling, one first generates the set ofsamples, {θ1, ...,θN} from a ‘target distribution’, q(θ) (which is relatively easy to sample from). One thencalculates a set of ‘importance weights’ according towi =π(θi)q(θi), i = 1, ..., N. (2)It can then be shown thatE[f(θ)] ≈N∑i=1f(θi)w˜i (3)wherew˜i =wi∑j wj, i = 1, ..., N (4)are referred to as ‘normalised importance weights’. If one then generates a new set of samples, {θ¯1, ..., θ¯N},by sampling with replacement from the original set while ensuring thatPr(θ¯ = θi) = w˜i (5)then it can also be shown that the new set of samples will be an approximate discrete representation of π(θ).The weights of the new samples are, as a result, approximately equal. This method, known as ‘resampling’,helps to remove samples with low weights and produce replicas of those with large weights. Resampling isoften performed when one’s estimates are dominated by relatively few samples (in other words, when themajority of samples have low weights). This is often referred to as ‘degeneracy’ in particle filter literature.Based on the work in [7], resampling occurs if the ‘effective sample size’, Nf , drops below a threshold (N/2in the current paper), where the effective sample size is defined asNeff =1∑i(w˜i)2. (6)It is important to note that the resampling procedure can be readily parallelised (see [8] for more details)and that this suitability for modern, fine grained computer architectures is a property of the algorithm that isproposed in the current paper.2.2 Previous workIn [9] it was suggested that importance sampling could help to address the issues raised at the beginning ofthe current paper. In the situation where one has obtained measurements up to time t, the target distributionis p(θ |z1:t). Using the prior, p(θ), as a proposal density, one’s importance weights at time t can then becalculated according towt =p(z1:t|θ)p(θ)p(θ)= p(z1:t|θ). (7)Assuming that the probability of witnessing separate measurements is independent, such that p(z1:t |θ) =∏ti=1 p(zi |θ), allows these importance weights to be calculated sequentially bywt = p(zt|θ)wt−1 (8)which prevents one from having to analyse the entire dataset every time a new measurement arrives. Unfor-tunately however, as more data arrives, the posterior tends to become very concentrated relative to the priorand many of the samples generated will be associated with very low importance weights. While resamplingcan help to remove these samples, it doesn’t help further exploration of the parameter space in regions ofhigh probability density. In [9] it was suggested that a ‘move step’ could be inserted to address this (possiblyusing MCMC). This step is, however, rather expensive to implement. In the current work it is suggested thatthis move step can be conducted efficiently using a SMC sampler, which is introduced in the next section.3 Sequential Monte Carlo samplersA brief introduction to SMC samplers is given here. A more detailed description is given in [10] while amore applied introduction can be found in [11].Say θk is the state of a system at iteration k and that πk(θk) is defined as the kth target distribution - thecurrent distribution of interest. (The observant reader may notice a change in notation here, where k’s havebeen employed instead of t’s - the reasons for this are explained later in the text.) π(θ1:k) is used to representthe joint distribution over all previous states. With an SMC sampler, rather than targeting πk(θk) directly,one instead targetsπk(θk)k∏k′=2L(θk′−1|θk′) (9)where L(θk′−1|θk′) (known as the ‘L-kernel’) is defined such that∫ {πk(θk)k∏k′=2L(θk′−1|θk′)}dθ1:k−1 = πk(θk). (10)Generating samples from the proposal distribution q(θ1:k) = q(θk |θk−1)q(θ1:k−1), importance weights areassigned according to:wik =πk(θik)∏kk′=2 L(θk′−1|θk′)q(θk |θk−1)q(θ1:k−1). (11)From this one can then go on to show thatwik = wik−1πk(θik)πk−1(θik−1)L(θik−1|θik)q(θik|θik−1)(12)such that, in other words, the importance weights can now be updated in a sequential manner as k increases.At first sight one may consider exchanging the indices k and t such that the target distribution isπt(θt) = p(θ |z1:t) (13)and using the SMC sampler to tackle the problem outlined in Section 1. This would, however, still requireanalysis of the full dataset every time new measurements arrive. In the current paper it is suggested that theSMC sampler can actually be used to facilitate the ‘move step’ of the method proposed in [9]. As will beshown, the advantages of this approach are that• The move step will only be conducted if the effective sample size drops below a pre-defined threshold.• The move step will be suitable for parallelisation.4 Proposed methodologySay one has a model of the formyt = f(xt,θ) (14)where xt represents the input to the model at time t. It is assumed that measurements are made according tozt = h(yt) + ǫ, ǫ ∼ N (0,Σǫ) (15)where the covariance matrix Σǫ is known and h can be a nonlinear function. The proposed methodology forsequential parameter estimation is outlined in algorithm 1.Algorithm 1 Proposed algorithm.Set k = 1, t = 1Sample {θ1t,k, ...,θNt,k} from the prior, p(θ)Find initial weights: wit,k = p(z1 |θit,k), i = 1, ..., Nwhile doNormalise weights: w˜i = wit,k∑j wjt,kEstimate quantities of interest.Find effective sample size: Neff = 1∑i(w˜i)2while Neff < N/2 doResample to get {θ1t,k, ...,θNt,k}Reset weights: wit,k = 1, i = 1, ..., NSet k = k + 1Generate new samples {θ1t,k, ...,θNt,k} from q(θt,k |θt,k−1)Find new weights: wit,k =p(z1:t | θit,k)p(z1:t |θit,k−1)L(θit,k−1 |θit,k)q(θit,k | θit,k−1), i = 1, ..., NNormalise weights: w˜i = wit,k∑j wjt,kFind effective sample size: Neff = 1∑i(w˜i)2end whilet = t+ 1Find new weights: wit,k = wit−1,kp(zt |θt,k)end whileKey points to note are that:• If degeneracy doesn’t occur (such that the effective sample size remains sufficiently large) then theSMC sampler will never be employed, and the algorithm will proceed in a similar manner to thatproposed in [9] (except without ‘move steps’).• The SMC sampler is only employed if degeneracy occurs. The SMC sampler helps to provide a furtherexploration of the parameter space in a way which is suitable for parallel processing.• The SMC sampler will run until the effective sample size has increased above the predefined threshold.• k is used to index the number of times the SMC sampler is run, while t indexes time.One of the advantages of this approach is that, if new data has relatively little effect on the geometry of theposterior, then the samples simply need to be reweighted. It is only when ‘information rich’ data arrives, anda large change in the posterior occurs, that the parameter samples need to be moved (and the entire datasetneeds to be reanalysed). The method therefore embodies the idea of ‘highly informative training data’ out-lined in [1].With regard to the L-kernel, it should be noted that a symmetric proposal distribution is utilised throughoutthe following examples. The L-kernel is chosen to beL(θit,k−1|θit,k) = q(θit,k−1|θit,k) (16)which then allows the weights of the samples generated using the SMC sampler to be calculated fromwit,k =p(z1:t |θit,k)p(z1:t |θit,k−1)L(θit,k−1 |θit,k)q(θit,k |θit,k−1)=p(z1:t |θit,k)p(z1:t |θit,k−1). (17)5 Results5.1 Example 1 - linear systemIn this first example, a linear model is considered, whereyt = θxt (18)and measurements are made according tozt = yt + ǫ, ǫ ∼ N (0, β−1ǫ ) (19)where βǫ is the noise precision. This is a convenient problem as, by choosing a Gaussian prior, p(θ) =N (θ;µ0, β−10 ), closed form expressions for the posterior distribution can be obtained. This therefore allowsone to check that the proposed method can follow the true solution. 100 ‘measurements’ were generated(sampling x from U [0, 1]) and the data was analysed sequentially. The prior moments, noise precision andtrue value of θ wereµ0 = 0, β0 = 1, βǫ = 100, θ = 1 (20)while N = 100 samples were employed. Proposals were drawn according toq(θit,k|θit,k−1) = N (θit,k; θit,k−1, 1). (21)Figure 1 shows that the algorithm performed well, and was able to closely match the true solution. It shouldbe noted that resampling - the reanalysis of the full dataset - was only needed when particularly ‘informative’data points were analysed (in other words, when large changes in the estimated mean and variance occurred).The remainer of the analysis could be performed very cheaply, without the SMC sampling step.0 20 40 60 80 1000.750.80.850.90.9511.051.1tMean0 20 40 60 80 10010−510−410−310−210−1tVariance  True SolutionSMC SamplerResampling EventFigure 1: Applying the proposed methodology to the identification of a linear system.5.2 Example 2 - comparison with MCMCIn this section, the previous example is reanalysed using a MCMC based method. Specifically, a variantof the well-known Transitional MCMC (TMCMC) [12] algorithm is proposed and applied. Here, when theeffective sample size drops below the predefined threshold, new samples are generated in a manner similarto that employed in the TMCMC algorithm. Specifically, Markov chains (which are ergodic and whose sta-tionary distributions are equal to the target), are ‘grown’ from the existing samples. The probability that eachchain will grow is equal to the normalised importance weights of the corresponding existing sample. Thistherefore encourages new samples to be generated in the region where existing samples are highly weighted.The problem described in the previous section was analysed 1000 times using the proposed SMC methodol-ogy as well as this new, TMCMC based methodology.The statistics of the ensemble of results are shown in Figure 2. It is clear that both methodologies are veryclosely matched. The advantage of the SMC approach, however, is that it will always be suitable for parallelprocessing (while the parallelisation of TMCMC depends on the distribution of importance weights). In anextreme case, for example, where one sample has a normalised weight equal to 1 and all others are equalto zero, the TMCMC approach would result in the generation of a single Markov chain which would beimpossible to parralelise.0 20 40 60 80 1000.950.960.970.980.9911.011.021.031.041.05tmean  0 20 40 60 80 10010−510−410−310−210−1tvariance  SMC (mean)SMC (3σ confidence bounds)TMCMC (mean)TMCMC (3σ confidence bounds)Figure 2: Identification of a linear system using the proposed SMC approach and the proposed MCMCapproach. Statistics are taken from an ensemble of 1000 simulations.5.3 Example 3 - simulated two-storey structureAs the final example in this paper, ‘measurement’ data is generated from a simulation of a linear two-storeystructure. The mass of each floor was set equal to 10000 kg. The vector of parameters to be estimated wereθ = (k1, k2, c1, c2)T where the k’s and c’s represent interstorey stiffness and damping terms respectively.Their true values, and details of the prior distributions that were utilised, are given in Table 1.Parameter True value Priork1 30× 106 N (k1; 25× 106, 2, 5 × 1012)k2 20× 106 N (k2; 25× 106, 2, 5 × 1012)c1 30× 103 Γ(c1; 1, 5 × 104)c2 20× 103 Γ(c2; 1, 5 × 104)Table 1: True parameter values and corresponding prior (SI units).To generate training data, the structure was subjected to a band-limited white noise base acceleration, withstandard deviation equal to 0.01 m/s2. The displacement time histories of the two floors were then corruptedwith Gaussian noise, such that the signal to noise ratio was approximately 10. Figure 3 shows how the pro-posed algorithm successfully converges to the true parameter values as more data is analysed. It should benoted that, again, as more data is analysed resampling is needed less frequently. This is because, after aninitial phase, new measurement data is contributing relatively little additional information.0 500 1000 1500 20001.522.533.5x 107k 1 (N/m)0 500 1000 1500 20001.522.533.5x 107k 2 (N/m)  Mean estimateResampling eventTrue value3σ confidence bounds0 500 1000 1500 2000−10123x 105tc 1 (Ns/m)0 500 1000 1500 2000−2−1012x 105c 2 (Ns/m)Figure 3: Parameter estimation of two-storey simulated structure using the proposed SMC methodology.6 ConclusionsIn this paper the authors investigate the parameter identification of systems from data sets which are con-tinuously growing. The aim is to provide an efficient methodology which can be used to track how one’sparameter estimates converge as more data is analysed, thus allowing one to decide when a sufficient quan-tity of data has been analysed. An algorithm is proposed which utilises a Sequential Monte Carlo sampler.The efficiency of the method stems from its ability to quickly analyse data which adds little additional infor-mation, thus allowing it to focus on data which is ‘information rich’. It is demonstrated that the proposedmethod is able to achieve a similar level of repeatability when compared to an alternative, Markov chainMonte Carlo method. The proposed Sequential Monte Carlo method is, however, better suited to parallelimplementation on modern computer architectures.References[1] P.L. Green, E.J. Cross, K. Worden. Bayesian system identification of dynamical systems using highlyinformative training data. Mechanical Systems and Signal Processing, 56:109–122, 2015.[2] P.L. Green. Bayesian system identification of a nonlinear dynamical system using a novel variant ofsimulated annealing. Mechanical Systems and Signal Processing, 52:133–146, 2015.[3] P.L. Green. Bayesian system identification of dynamical systems using large sets of training data: AMCMC solution. Probabilistic Engineering Mechanics, 42:54–63, 2015.[4] D. Crisan, A. Doucet. A survey of convergence results on particle filtering methods for practitioners.Signal Processing, IEEE Transactions on, 50(3):736–746, 2002.[5] C.M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.[6] D.J.C. MacKay. Information theory, inference and learning algorithms. Cambridge university press,2003.[7] A. Kong. A note on importance sampling using standardized weights. University of Chicago, Dept. ofStatistics, Tech. Rep, 348, 1992.[8] S. Maskell, B. Alun-Jones, M. Macleod. A single instruction multiple data particle filter. In NonlinearStatistical Signal Processing Workshop, 2006 IEEE, pages 51–54. IEEE, 2006.[9] N. Chopin. A sequential particle filter method for static models. Biometrika, 89(3):539–552, 2002.[10] P. Del Moral, A. Doucet, A. Jasra. Sequential monte carlo samplers. Journal of the Royal StatisticalSociety: Series B (Statistical Methodology), 68(3):411–436, 2006.[11] S. Maskell. An application of Sequential Monte Carlo samplers: An alternative to particle filtersfor non-linear non-Gaussian sequential inference with zero process noise. In Data Fusion & TargetTracking Conference (DF&TT 2012): Algorithms & Applications, 9th IET, pages 1–8. IET, 2012.[12] J. Ching, Y.C. Chen. Transitional Markov chain Monte Carlo method for Bayesian model updating,model class selection, and model averaging. Journal of engineering mechanics, 133(7):816–832, 2007.

Parameter estimation from big data using a sequential monte carlo sampler

Abstract

Similar works

Full text

Available Versions

University of Liverpool Repository