We propose a novel method for parameter inference that builds on the current research in gradient matching surrogate likelihood spaces. Adopting a three phase technique, we demonstrate that it is possible to obtain parameter estimates of limited bias whilst still adopting the paradigm of the computationally cheap surrogate approximation

Husmeier, Dirk

Lazarus, Alan

Papamarkou, Theodore

English

Enlighten

Inference in complex systems usingmulti-phase MCMC sampling with gradientmatching burn-inAlan Lazarus1, Dirk Husmeier1, Theodore Papamarkou11 School of Mathematics and Statistics, University of Glasgow, UKE-mail for correspondence: a.lazarus.1@research.gla.ac.ukAbstract: We propose a novel method for parameter inference that builds onthe current research in gradient matching surrogate likelihood spaces. Adoptinga three phase technique, we demonstrate that it is possible to obtain parameterestimates of limited bias whilst still adopting the paradigm of the computationallycheap surrogate approximation.Keywords: Parameter inference; Delayed Rejection Adaptive Metropolis; Sur-rogate likelihood; Markov Chain Monte Carlo; Gradient matching1 IntroductionStatistical inference in nonlinear differential equations (DE) is challenging. Thelog-likelihood landscape is typically multimodal and every parameter adaptation,e.g. in an MCMC simulation, requires a computationally expensive numerical in-tegration of the DEs. Using numerical methods to solve the equations resultsin prohibitive computational cost; particularly when one adopts a Bayesian ap-proach in sampling parameters from a posterior distribution. Alternatively, onecan try to reduce this computational complexity by obtaining an interpolantto the data from which one can obtain a comparative objective function thatmatches the gradients of the interpolant and the DEs. By sampling on this cheaprepresentative likelihood surface, bias is introduced to the modelling problem.Current research focuses on reducing this bias by introducing a regularising feed-back mechanism from the DEs back to the interpolation scheme (e.g. Niu et al.2016). The idea is to make the interpolant maximally consistent with the DEs.Although this paradigm has proved to improve performance over na¨ıve gradientmatching, the feedback loop fails to fully eradicate bias in the final estimate.This paper was published as a part of the proceedings of the 32nd Interna-tional Workshop on Statistical Modelling (IWSM), Johann Bernoulli Institute,Rijksuniversiteit Groningen, Netherlands, 3–7 July 2017. The copyright remainswith the author(s). Permission to reproduce or extract any parts of this abstractshould be requested from the author(s).Lazarus et al. 53For this reason, a natural progression would be to sample from the true likeli-hood space whilst reducing computational complexity in the discarded burninsteps. Assuming this hypothesis, we postulate the use of a surrogate likelihood inthe burnin phase alone. Through an example possessing multimodal likelihood,we will show the ability of the algorithm to avoid any local entrapment whilstobtaining accurate parameter estimates.2 MethodConsider time-dependent observations y(t) = x(t) + —where x(t) denotes thesignal and  independent additive zero mean Gaussian noise with variance pa-rameter σ2—whose signals are governed by a system of differential equations:dx(t)dt= f(x(t),θ) (1)dependent on some (partially) unknown parameters θ. Assuming Gaussian noise,we place a GP prior on the latent variable xx(t) ∼ GP(0, k(t, t′)), (2)leading us to a Gaussian distribution, p(xi|φi) = N (xi|0,Ki) for an arbitraryset of time points T = {t1, ..., tn} with entries of Ki given by evaluating kernelfunction k at each element of T ×T (Rasmussen and Williams, 2006). Under ourassumption of Gaussian noise, we consider the joint distribution, p(y,x|φ, σ).Marginalising over latent variables x provides a zero mean distribution for theobservations:y ∼ N (0,K + σ2I) (3)(see Dondelinger et al. (2013) for details). Considering the joint distribution be-tween our signal and observed values, we may implement an elementary trans-formation of a Gaussian distribution to obtain the posterior distribution for oursignal with mean given by:µ(τ) = k(τ, T )(K + σ2I)−1y, (4)where k(τ, T ) denotes evaluation of the kernel function at τ over T . Subsequently,firstly estimating the hyperparameters via ML, we adopt the mean of the posterioras a representation of our signal x. This allows us to proceed under the suppositionthat we have a fixed interpolant for the true signal. Given that the derivative ofeq. 4,∂µ(τ)∂τ= k′(τ, T )(K + σ2I)−1y, (5)is the mean of the Gaussian distributed DE derivative (see section 7.5 of Van-hatalo et al., 2015), we may consider:f(xˆ(t),θ) ∼ N(dxˆ(t)dt, γ2I), (6)54 MCMC parameter inference with surrogate burninwhere γ2 is a fictitious noise term (assumed equal for both gradients) and xˆ isgiven by eq. 4. Contrary to the work done by Dondelinger et al. (2013), fixingthe GP hyperparmeters φ and the interpolant xˆ allows us to abandon the Gibbssampling routine at this stage of the algorithm, further reducing the overall com-putational burden. The corresponding negative log-likelihood term provides agradient matching objective function:pi(θ) = n logγ2 +12γ2∣∣∣∣∣∣∣∣dxˆ(t)dt − f(xˆ(t),θ)∣∣∣∣∣∣∣∣2 (7)which gives a representative computationally tractable surface as a surrogate forthe log-likelihood. This involves the term γ2 representing the mismatch betweenthe gradient obtained from the differential equation and that from explicit differ-entiation of the GP posterior mean. This parameter will be sampled throughoutthe surrogate burnin phase of the algorithm. The proposed sampling scheme in-volves three phases. In the initial burnin phase, samples are drawn from thesurrogate distribution in eq. 7 using a Delayed Rejection Adaptive Metropolis1(DRAM) scheme (Haario et al. (2006)). Assuming a degree of similarity betweenthe surrogate and true likelihood surfaces, this drives the sampler towards theglobal minimum of the true likelihood function until a PSRF2 value of 1.1 hasbeen achieved. From here, we initialise a corrective phase in the true likelihoodspace, correcting for any bias introduced by the inconsistencies between the sur-rogate and true likelihood spaces. Sampling with DRAM, this phase is concludedupon obtainment of a PSRF equal to 1.1. The proceeding sampling phase repli-cates this corrective phase with sampling steps recorded until we achieve a PSRFvalue of 1.05. The stepwise decrease in target PSRF values allows time for theadaptive component of AM to learn the topology and adjust the covariance ac-cordingly. We adopt an uninformative Inv-Gamma(0.001,0.001) prior for σ2 andγ2 and a Ga(4, 0.5) prior for the parameters of the DE. All parameters are sam-pled on the log scale to account for the positivity constraint.3 ResultsWe assess performance on the following DE model of circadian oscillation3:dp1dt=k136 + k2p2− k3, dp2dt= k4p1 − k5 (8)This is a notoriously challenging problem due to the extreme multimodality ofthe likelihood. Following Girolami et al. (2010), we focus on the inference of twoparameters (k3 and k4), setting the other parameters and initial conditions to1Obtained using the adaptive Metropolis component of DRAM with modM-CMC function in the FME package in R.2Obtained at intervals of 20 steps using the gelman.diag function from thecoda package in R.3We used the same differential equations as in Girolami et al 2010. The actualGoodwin oscillator is of a slightly different form, where the terms k3 and k5 arereplaced by k3p1 and k5p2, respectively.Lazarus et al. 55the same fixed values as in Girolami et al. (2010). Five sets of initial parame-ter values for k3 and k4 were obtained using a Sobol sequence over the domain[0, 5]2. Figure ?? shows the chain moving through the k3-k4 parameter domain.Comparing with the traditional method, we observe the ability of the proposedmethod to evade the various local minima. PSRF values of 17.1, 10.4 and 12.3were obtained for k3, k4 and σ2 simulations respectively after 10000 steps us-ing the traditional DRAM method in true likelihood space. Comparatively, theproposed method required 1690 steps in surrogate space, 1690 in the correctivephase and 1010 in the sampling phase to achieve a PSRF of 1.05. The numberof numerical integration steps required are given in Table ?? for each of the tenDRAM chains. In Figure ??, boxplots are given that provide the distribution of 500  1000  1000  1500  2000  2500  3000  3500  4000  4500  5000  5500  6000  6500  7000  7500  8000  8500  9000  9500  10000  10500  11000  11500  12000  12500  13000  13500  14000  14500  15000  15500  16000  16500  17000  17500  18000  18500  19000  19500  20000  20500  21000  21500  22000  22500  23000  23500  24000  24500  25000  25500  26000  26500  27000  27500  28000  28500  29000  29500  30000  30500  31000 1 2 3 4 512345k3k 4●●●●●●●●●●●x −120964828  −99652917  −76565013  −48149132  −46373139  −39269169  −32165199  −25573634  −22591038  −17718481  −15798989  −14411049  −13407007  −11162677  −9006940  −8180082  −7500877  −7176040  −5433731  −5079364  −4990772  −4724996  −3720954  −3661893  −3425647  −3041749  −3012218  −2864565  −2273952  −2126299  −2096768  −1978646  −1890054  −1653808  −1447094  −1417563  −1269910  −1240379  −1151787  −1092726  −1004134  −974603.5  −945072.9  −886011.6  −767889  −708827.7  −679297  −620235.7  −590705.1  −561174.4  −531643.8  −502113.1  −472582.5  −443051.8  −413521.2  −383990.5  −354459.9  −324929.2  −295398.6  −295398.6  −265867.9  −265867.9  −236337.3  −236337.3  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −177276  −177276  −147745.3  −118214.7  −88684.05 1 2 3 4 512345k3k 4●●●●●●●●●●●x −120964828  −99652917  −76565013  −48149132  −46373139  −39269169  −32165199  −25573634  −22591038  −17718481  −15798989  −14411049  −13407007  −11162677  −9006940  −8180082  −7500877  −7176040  −5433731  −5079364  −4990772  −4724996  −3720954  −3661893  −3425647  −3041749  −3012218  −2864565  −2273952  −2126299  −2096768  −1978646  −1890054  −1653808  −1447094  −1417563  −1269910  −1240379  −1151787  −1092726  −1004134  −974603.5  −945072.9  −886011.6  −767889  −708827.7  −679297  −620235.7  −590705.1  −561174.4  −531643.8  −502113.1  −472582.5  −443051.8  −413521.2  −383990.5  −354459.9  −324929.2  −295398.6  −295398.6  −265867.9  −265867.9  −236337.3  −236337.3  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −206806.6  −177276  −177276  −147745.3  −118214.7  −88684.05 1 2 3 4 512345k3k 4●●xx●●●●●●●●●FIGURE 1: Ten chains generated with the proposed method shown in theparameter domain of the negative surrogate log-likelihood space (topleft)and negative log-likelihood space (topright). The bottom plot shows simu-lations generated using the traditional method. The true parameter valueis given by a point at (2,1). The black crosses denote the final point of eachchain.TABLE 1: Number of numerical integration steps (N) for the traditionalmethod. The number required in the proposed scheme is 2700.Chain Index 1 2 3 4 5 6 7 8 9 10N 25583 29778 29835 28452 29708 29672 29768 29530 29859 29847bias in our sampled parameter estimates for each of the five chains. Figure ??56 MCMC parameter inference with surrogate burninprovides RMS deviation in function space obtained using eq. 9,RMSfunction =√1n∣∣∣∣x− xˆ∣∣∣∣2 (9)where x denotes the true signal and xˆ denotes the numerical solution of the DEfor one parameter sample from the sample phase of the multiphase approach andthe post burnin period of the traditional method. This provides a measure of thepredictive accuracy of the MCMC samples.●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●−0.01−0.00500.0050.01k3 k3 k4 k4Bias●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●01234k3 k4BiasFIGURE 2: Boxplots showing the distribution of bias for both methods (left)where red boxes give the bias of the standard DRAM samples (withoutoutliers) and the black boxes give the bias of the proposed method. Theplot on the right gives the bias in both parameters for the DRAM methodwith outliers included.●●●●●●●●●●●●●Novel Method Standard Method0102030RMS Value ●●●●●●●●●●●●●Novel Method Standard Method00.10.20.30.4RMS ValueFIGURE 3: Functional RMS comparison between the proposed method(black) and DRAM (red). On the left, we include outliers (from DRAM)and, on the right, these are removed to enable better scalability of theplots for comparison. The red dotted line denotes a functional RMS equalto zero.4 ConclusionOur work considers the sampling in DE parameter inference as a computationallyefficient three-phase scheme that achieves low levels of bias in sampled parameterestimates (Figure ??). Achieving a PSRF of 1.05, we observe the ability of theLazarus et al. 57algorithm to converge in the parameter space of the circadian oscillator systemof equations; a model for which the standard DRAM procedure fails to replicatethis success (bottom of Figure ??). Considering the results in function space, weobserve the superior performance of the proposed method compared with the tra-ditional DRAM method, showing that the performance improvement is witnessedin both domains of study. These features, along with the vast improvement incomputational efficiency, demonstrate improved parameter inference comparedwith the traditional method.Acknowledgments: This project was supported by a grant from the Engi-neering and Physical Sciences Research Council (EPSRC) of the UK, grant ref-erence number EP/L020319/1.ReferencesDondelinger, F. et al. (2013). ODE parameter inference using adaptive gradientmatching with Gaussian processes. Proceedings of Machine Learning Re-search, Volume 31, pp 216-228.Girolami, M. et al. (2010). System identification and model ranking: the Bayesianperspective. Learning and Inference in Computational Systems Biology.MIT Press, pp. 201-230.Haario, H. et al. (2006). DRAM: Efficient adaptive MCMC. Statistics and Com-puting, Volume 16, Issue 4, pp 339-354.Niu, M. et al. (2016). Fast Parameter Inference in Nonlinear Dynamical Systemsusing Iterative Gradient Matching. Proceedings of Machine Learning Re-search, Volume 48, pp 1699-1707.Rasmussen, C.E. and Williams, C.K.I. (2006). Gaussian Processes for MachineLearning. MIT Press.Vanhatalo, J. et al. (2015). Bayesian Modeling with Gaussian Processes usingthe GPstuff Toolbox. MIT Press.

Inference in Complex Systems Using Multi-Phase MCMC Sampling With Gradient Matching Burn-in

http://eprints.gla.ac.uk/144287/7/144287.pdf

Inference in Complex Systems Using Multi-Phase MCMC Sampling With Gradient Matching Burn-in

Abstract

Similar works

Full text

Available Versions

Enlighten