Bell, Andrew

Ignatieva, Anastasia

Worton, Bruce

English

In this paper we consider statistical modelling of volcanic earthquake data. In particular, we investigate the use of Bayesian analysis with Markov Chain Monte Carlo (MCMC) to estimate the parameters of point process models, and make inferences on the models, applied to data collected from the Tungurahua volcano in Ecuador

Edinburgh Research Explorer

     Edinburgh Research Explorer                                      Comparison of MCMC approaches with an application to volcanoearthquake processesCitation for published version:Ignatieva, A, Bell, A & Worton, B 2018, Comparison of MCMC approaches with an application to volcanoearthquake processes. in Proceedings of the 33rd International Workshop on Statistical Modelling.University of Bristol, 33rd International Workshop on Statistical Modelling, Bristol, United Kingdom,16/07/18. <https://digitalcommons.usf.edu/cgi/viewcontent.cgi?article=1019&context=siv>Link:Link to publication record in Edinburgh Research ExplorerDocument Version:Peer reviewed versionPublished In:Proceedings of the 33rd International Workshop on Statistical ModellingGeneral rightsCopyright for the publications made accessible via the Edinburgh Research Explorer is retained by the author(s)and / or other copyright owners and it is a condition of accessing these publications that users recognise andabide by the legal requirements associated with these rights.Take down policyThe University of Edinburgh has made every reasonable effort to ensure that Edinburgh Research Explorercontent complies with UK legislation. If you believe that the public display of this file breaches copyright pleasecontact openaccess@ed.ac.uk providing details, and we will remove access to the work immediately andinvestigate your claim.Download date: 08. Jun. 2022Comparison of MCMC approaches with anapplication to volcano earthquake processesAnastasia Ignatieva1, Andrew F. Bell2, Bruce J. Worton11 School of Mathematics and Maxwell Institute for Mathematical Sciences, TheUniversity of Edinburgh, Edinburgh, UK2 School of GeoSciences, The University of Edinburgh, Edinburgh, UKE-mail for correspondence: Bruce.Worton@ed.ac.ukAbstract: In this paper we consider statistical modelling of volcanic earthquakedata. In particular, we investigate the use of Bayesian analysis with Markov ChainMonte Carlo (MCMC) to estimate the parameters of point process models, andmake inferences on the models, applied to data collected from the Tungurahuavolcano in Ecuador.Keywords: Bayesian modelling; Eruption forecasting; Point processes.1 IntroductionThis paper aims to use statistical modelling to describe the occurrence ofvolcanic earthquakes. The main approach taken is that of using Bayesiananalysis with Markov Chain Monte Carlo (MCMC) to fit point processmodels to the available data, collected from the Tungurahua volcano inEcuador.2 Dataset and modellingThis dataset was recorded in July 2013, and consists of a series of eventtimes which were picked from a stretch of seismic data to identify theindividual earthquakes. The dataset was examined in a study by Bell et al.(2018).The events started at 6:00 on 13 July, and the eruption occurred at 11:46on 14 July. The event rate grew increasingly up until eruption. Plots ofthe data show that the event rate grows at an increasing rate up to theThis paper was published as a part of the proceedings of the 33rd Inter-national Workshop on Statistical Modelling (IWSM), University of Bristol, UK,16-20 July 2018. The copyright remains with the author(s). Permission to repro-duce or extract any parts of this abstract should be requested from the author(s).2 Modelling volcano earthquake processeseruption, with the inter-spike interval (ISI) duration changing from over10 minutes to below 30 seconds. The ISIs are “quasi-periodic”, being moreregular than would be seen if the events followed a Poisson process, andthus not independent (Bell et al., 2018).Applying a material failure approach to describe the physical processesleading a volcanic system to an eruption, the accelerating rate of earth-quakes is described by a power law relationship (Bell et al., 2018):λ(t) = k(tf − t)−p,where k is a constant (related to the amplitude of the signal), tf is thetime of eruption, and p = 1a−1 is a parameter describing the non-linearityof acceleration. At time tf , the rate becomes instantaneously infinite, rep-resenting the eruption (Bell et al., 2018). In the model, λ(t) is the intensityused in the inhomogeneous gamma (IG; parameter α) point process.Details of the MCMC implementation in PyMC3 were investigated, includ-ing the sampling method used and the initialisation process. Attributes ofthe MCMC chain, such as convergence, were examined. Posterior checkswere performed using simulated data, to sense check whether the modelappears appropriate. The fit of the model was assessed further using sta-tistical methods.The MCMC approaches considered included:• No-U-Turn sampler (Hoffman and Gelman, 2014);• Metropolis;• Slice sampling (Neal, 2003).Alternative models were also investigated, and their fits compared to thatof the given inhomogeneous gamma model: inhomogeneous Poisson (IP),inhomogeneous inverse Gaussian (IIG) and inhomogeneous Weibull (IW)models.3 ResultsFigure 1 gives comparison of the MCMC trace plots for the methods forthe IG model. Figures 2 and 3 show posterior plots using MCMC sampledvalues for the IG model. A Kolmogorov-Smirnov goodness of fit approach(Barbieri et al., 2001; Ogata, 1988) gives an effective method of comparisonof various possible models.Ignatieva et al. 3Figure 2.5: Trace plots, first 2k iterations (NUTS, Metropolis, Slice)2.6 NUTS further analysis2.6.1 ConvergenceTo test convergence, 10 chains were run in parallel, starting at randomly selected initialpoints. It can be seen from the trace plot in Figure 2.6 that all of the chains appear to haveconverged to the same distribution. Moreover, the parameter means of each of the chains(shown in Table 2.1) did not vary by more than 0.1 for α, tf , p, and by no more than 1 fork, from the values presented in Table 2.2, which also indicates convergence.Table 2.1: Mean values of parameters for each of 10 randomly started chainsChain α k tf p1 2.476 340.827 1.422 1.5082 2.474 339.679 1.420 1.5013 2.475 340.564 1.421 1.5014 2.476 338.639 1.418 1.4955 2.476 341.698 1.422 1.5066 2.477 339.795 1.420 1.5007 2.475 340.430 1.421 1.5038 2.474 340.166 1.420 1.5009 2.475 339.651 1.420 1.50010 2.475 341.572 1.423 1.509Gelman-Rubin statisticThe Gelman-Rubin statistic is very close to 1 for all the parameters, indicating that thereis no clear evidence against convergence (Gelman and Rubin, 1992).R̂(α) = 1.00001, R̂(k) = 1.00026, R̂(p) = 1.00021, R̂(tf ) = 1.0002217FIGURE 1. race plots for 2,000 iterations: IG model. No-U-Turn sampler (left),Metropolis (middle), Slice (right).Figure 2.8: Autocorrelation plot, 20k iterations, excluding burn-inFigure 2.9: Posterior plots using MCMC sampled values19FIGURE 2. Posterior plots using MCMC sampled values: IG model.4 Modelling volcano earthquake processesbe explained by these parameters being correlated with each other, as can be noted fromFigure 2.9 which shows scatter plots of the parameter values from the samples drawn fromthe posterior (generated from the code in Appendix A.1 provided by Dr. A. Bell (2017)).Lower (higher) values of k tend to occur with lower (higher) values of p and tf .For α, there is no such visible autocorrelation with the other parameters, as the scatterplots of k, tf , and p with α shows a lack of discernible trends.The NUTS sampler works effectively in such situations where some of the parameters arecorrelated, however this significantly slows down other samplers like the Metropolis and Slicesampler, in line with the earlier findings.2.6.3 ThinningThinning is a process by which every k-th iteration of the chain is kept, and the rest arediscarded. This has the effect of reducing autocorrelation within the chain, and reduces MonteCarlo error. However, the discarded values still carry some information about the posteriordistributions. Thus, given that there are no practical restrictions on memory storage in thiscase, thinning is not considered necessary.2.6.4 ResultsThe posterior distributions, with 95% highest posterior density intervals (HPDIs), areshown below in Figure 2.10, with values given in Table 2.2. The resulting HPDIs are reason-ably narrow, and the distributions have regular unimodal shapes.Figure 2.10: Posteriors (excluding burn-in), 20k iterations NUTSTable 2.2: IG model parametersParameter Mean Std. dev. 95% HPDIα 2.476 0.117 [2.252, 2.710]k 339.893 33.207 [285.102, 408.778]tf 1.420 0.060 [1.310, 1.539]p 1.501 0.175 [1.175, 1.852]20FIGURE 3. Posteriors with HPDs: IG model.4 ConclusionsAn IG model was found to produce satisfactory results. It was demonstratedthat the MCMC chain appears to converge to the correct stationary distri-bution, providing reasonable posterior estimates. From review of simulateddata, and Q-Q and K-S plots, it was found that the IG model fits the July2013 data very well. A small number of outliers (around 5% of the data)was noted, and found to correspond to spikes with long preceding ISIs.Some lack of fit was also found in the middle quantiles of the K-S plot,however this only slightly breached the 95% error bounds.Acknowledgments: We are particularly grateful to The Instituto Geosicoof the Escuela Politecnica Nacional (IGEPN) of Ecuador for all their hardmonitoring work and providing the data.ReferencesBarbieri, R., Quirk, M.C., Frank, L.M., Wilson, M.A., and Brown, E.N.(2001). Construction and analysis of non-Poisson stimulus-responsemodels of neural spiking activity. Journal of Neuroscience Methods,105, 25 – 37.Bell, A.F., Naylor, M., Hernandez, S., Main, I.G., Gaunt, H.E., Mothes, P.,Ignatieva et al. 5and Ruiz, M. (2018). Volcanic eruption forecasts from acceleratingrates of drumbeat long-period earthquakes. Geophysical Research Let-ters, 45, 1339 – 1348.Hoffman, M.D. and Gelman, A. (2014). The No-U-Turn sampler: adaptivelysetting path lengths in Hamiltonian Monte Carlo. Journal of MachineLearning Research, 15, 1593 – 1623.Neal, R.M. (2003). Slice sampling. Annals of Statistics, 31, 705 – 741.Ogata, Y. (1988). Statistical models for earthquake occurrences and resid-ual analysis for point processes. Journal of the American StatisticalAssociation, 83, 9 – 27.

Comparison of MCMC approaches with an application to volcano earthquake processes

https://www.research.ed.ac.uk/files/63885694/IWSM97.pdf

Comparison of MCMC approaches with an application to volcano earthquake processes

Abstract

Similar works

Full text

Available Versions

Edinburgh Research Explorer