Spectral Density-Based and Measure-Preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs
SSpectral Density-Based and Measure-Preserving ABCfor partially observed diffusion processes
An illustration on Hamiltonian SDEs
Evelyn Buckwar, Massimiliano Tamborrino, Irene TubikanecInstitute for StochasticsJohannes Kepler University Linz, Austria
Abstract
Approximate Bayesian Computation (ABC) has become one of the major tools of likelihood-free statistical inference in complex mathematical models. Simultaneously, stochastic differentialequations (SDEs) have developed to an established tool for modelling time dependent, real worldphenomena with underlying random effects. When applying ABC to stochastic models, two ma-jor difficulties arise. First, the derivation of effective summary statistics and proper distances isparticularly challenging, since simulations from the stochastic process under the same parameterconfiguration result in different trajectories. Second, exact simulation schemes to generate trajec-tories from the stochastic model are rarely available, requiring the derivation of suitable numericalmethods for the synthetic data generation. To obtain summaries that are less sensitive to theintrinsic stochasticity of the model, we propose to build up the statistical method (e.g., the choiceof the summary statistics) on the underlying structural properties of the model. Here, we focuson the existence of an invariant measure and we map the data to their estimated invariant den-sity and invariant spectral density. Then, to ensure that these model properties are kept in thesynthetic data generation, we adopt measure-preserving numerical splitting schemes. The derivedproperty-based and measure-preserving ABC method is illustrated on the broad class of partiallyobserved Hamiltonian type SDEs, both with simulated data and with real electroencephalography(EEG) data. The proposed ingredients can be incorporated into any type of ABC algorithm anddirectly applied to all SDEs that are characterised by an invariant distribution and for which ameasure-preserving numerical method can be derived.
Keywords
Approximate Bayesian Computation, Likelihood-free inference, Stochastic differential equations,Numerical splitting schemes, Invariant measure, Neural mass models
Acknowledgements
This research was partially supported by the Austrian Science Fund (FWF): W1214-N15, projectDK14.
Over the last decades, SDEs have become an established and powerful tool for modelling timedependent, real world phenomena with underlying random effects. They have been successfullyapplied to a variety of scientific fields, ranging from biology over finance, to physics, chemistry,neuroscience and others. Diffusion processes obtained as solutions of SDEs are typically charac-terised by some underlying structural properties whose investigation and preservation is crucial. a r X i v : . [ s t a t . C O ] J u l xamples are boundary properties, symmetries or the preservation of invariants or qualitative be-haviour such as the ergodicity or the conservation of energy. Here, we focus on a specific structuralproperty, namely the existence of a unique invariant measure. Besides the modelling, it is of pri-mary interest to estimate the underlying model parameters. This is particularly difficult when themultivariate stochastic process is only partially observed through a 1-dimensional function of itscoordinates (the output process), a scenario that we tackle here. Moreover, due to the increasingcomplexity of SDEs, needed to understand and reproduce the real data, the underlying likelihoodis often unknown or intractable. Among several likelihood-free inference approaches, we focus onthe simulation-based ABC method. We refer to Marin et al. (2012) and to the recently publishedbook “Handbook of Approximate Bayesian Computation” for an exhaustive discussion (Sissonet al. 2018).ABC has become one of the major tools for parameter inference in complex mathematicalmodels in the last decade. The method is based on the idea of deriving an approximate posteriordensity targeting the true (unavailable) posterior by running massive simulations from the modelto replace the intractable likelihood. It was first introduced in the context of population genetics;see, e.g., Beaumont et al. (2002). Since then, it has been successfully applied in a wide range offields; see, e.g., Barnes et al. (2012); Blum (2010a); Boys et al. (2008); McKinley et al. (2017);Moores et al. (2015); Toni et al. (2009). Moreover, ABC has also been proposed to infer parame-ters from time series models (see, e.g., Drovandi et al. 2016; Jasra 2015), state space models (see,e.g., Martin et al. 2019; Tancredi 2019) and SDE models (see, e.g., Kypraios et al. 2017; Maybanket al. 2017; Picchini 2014; Picchini and Forman 2016; Picchini and Samson 2018; Sun et al. 2015;Zhu et al. 2016). Several advanced ABC algorithms have been proposed in the literature, suchas, ABC-SMC, ABC-MCMC, sequential-annealing ABC, noisy ABC; see, e.g., Fan and Sisson(2018) and the references therein for a recent review. The idea of the basic acceptance-rejectionalgorithm is to keep a sampled parameter value from the prior as a realisation from the approxi-mate posterior, if the distance between the summary statistics of the synthetic dataset, which isgenerated conditioned on this parameter value, and the summaries of the original reference data issmaller than some tolerance level. The goal of this paper is to illustrate how building up the ABCmethod on the structural properties of the underlying SDE and using a numerical method capableto preserve them in the generation of the data from the model leads to a successful inference evenwhen applying ABC in this basic acceptance-rejection form.The performance of any ABC method depends heavily on the choice of “informative enough”summary statistics, a suitable distance measure and a proper tolerance level (cid:15) . The quality of theapproximation improves as (cid:15) decreases, and it has been shown that, under some conditions, theapproximated ABC posterior converges to the true one when (cid:15) → (cid:15) decreases. A possibility is to use ad-hocthreshold selection procedures; see, e.g., Barber et al. (2015); Blum (2010b); Lintusaari et al.(2017); Prangle et al. (2014); Robert (2016). Here, we fix the tolerance level (cid:15) as a percentile ofthe calculated distances. This is another common practice used, for example, in Beaumont et al.(2002); Biau et al. (2015); Sun et al. (2015); Vo et al. (2015). Instructions for constructing effectivesummaries and distances are rare and they depend on the problem under consideration; see, e.g.,Fearnhead and Prangle (2012) for a semi-automatic linear regression approach, Jiang et al. (2017)for an automatic construction approach based on training deep neural networks and Blum (2010b);Prangle (2018) for two recent reviews. To avoid the information loss caused by using non-sufficientsummary statistics another common procedure is to work with the entire dataset; see, e.g., Jasra(2015); Sun et al. (2015). This requires the application of more sophisticated distances d such asthe Wasserstein metric (Bernton et al. 2019; Muskulus and Verduyn-Lunel 2011) or other distancesdesigned for time series; for an overview see, e.g., Mori et al. (2016) and the references therein.When working with stochastic models, simulations from the stochastic simulator, conditionallyto the same parameter configuration, yield different trajectories. To consider summary statisticsthat are less sensitive to the intrinsic stochasticity of the model (Wood 2010), we choose thembased on the structural property of an underlying invariant measure. The idea is to map the data,i.e., the realisations of the output process, to an object that is invariant for repeated simulationsunder the same parameter setting and that reacts sensitive to small changes in the parameters.2n particular, we map the data to their estimated invariant density and invariant spectral density,taking thus the dependence structure of the dynamical model into account. The distance measurecan then be chosen according to the mapped data.As other simulation-based statistical methods, e.g., MCMC, SMC or machine learning algo-rithms, ABC relies on the ability of simulating data from the model. However, the exact sim-ulation from complex stochastic models is rarely possible, and thus numerical methods need tobe applied. This introduces a new level of approximation into the ABC framework. When thestatistical method is build upon the structural properties of the underlying model, the successfulinference can only be guaranteed when these properties are preserved in the synthetic data gen-erated from the model. However, the issue of deriving a property-preserving numerical methodwhen applying ABC to SDEs is usually seen as not so relevant, and it is usually recommendedto use the Euler-Maruyama scheme or one of the higher order approximation methods describedin Kloeden and Platen (1992); see, e.g., Picchini (2014); Picchini and Forman (2016); Picchiniand Samson (2018); Sun et al. (2015). In general, these standard methods do not preserve theunderlying structural properties of the model; see, e.g., Ableidinger et al. (2017); Malham andWiese (2013); Moro and Schurz (2007); Strømmen Melbø and Higham (2004).Here, we propose to apply structure-preserving numerical splitting schemes within the ABCalgorithm. The idea of these methods is to split the SDE into explicitly solvable subequations andto apply a proper composition of the resulting exact solutions. Standard procedures are, for exam-ple, the Lie-Trotter method and the usually more accurate Strang approach; see, e.g., Leimkuhleret al. (2016). Since the only approximation enters through the composition of the derived explicitsolutions, numerical splitting schemes usually preserve the structural properties of the underlyingSDE and accurately reproduce its qualitative behaviour. Moreover, they usually have the sameorder of convergence as the frequently applied Euler-Maruyama method and are likewise efficient.We refer to Blanes et al. (2009) and Mclachlan and Quispel (2002) for an exhaustive discussion ofsplitting methods for broad classes of ordinary differential equations (ODEs), which partially havealready been carried over to SDEs; see, e.g., Misawa (2001) for a general class of SDEs, Ableidingerand Buckwar (2016) for the stochastic Landau-Lifshitz equations, Br´ehier and Gouden`ege (2019)for the Allen-Cahn equation and Ableidinger et al. (2017) for Hamiltonian type SDEs.The main contribution of this work lies in the combination of the proposed invariant measure-based summary statistics and the measure-preserving numerical splitting schemes within the ABCframework. We demonstrate that a simulation-based inference method, here ABC, can only per-form well if the underlying simulation method preserves the structural properties of the SDE.While the use of preserving splitting schemes within the ABC method yield successful results, ap-plying a general purpose numerical method, such as the Euler-Maruyama discretisation, may resultin seriously wrong inferences. We illustrate the proposed Spectral Density-Based and Measure-Preserving ABC method on the class of stochastic Hamiltonian type equations for which theexistence of an underlying unique invariant distribution and measure-preserving numerical split-ting schemes have been already intensively studied in the literature; see, e.g., Ableidinger et al.(2017); Mattingly et al. (2002); Leimkuhler and Matthews (2015); Milstein and Tretyakov (2004).Hamiltonian type SDEs have been investigated in molecular dynamics, where they are typicallyreferred to as Langevin equations; see, e.g., Leimkuhler and Matthews (2015). Recently, theyhave also received considerable attention in the field of neuroscience as the so-called neural massmodels (Ableidinger et al. 2017).The paper is organised as follows. In Section 2, we recall the acceptance-rejection ABC setting.We introduce the invariant measure-based summary statistics and propose a proper distance. Wethen discuss the importance of considering measure-preserving numerical schemes for the syntheticdata generation when exact simulation methods are not applicable and provide a short introductionto numerical splitting methods. In Section 3, we introduce Hamiltonian type SDEs and recall twosplitting integrators preserving the invariant measure of the model. In Section 4, we validate theproposed method by investigating the stochastic harmonic oscillator, for which exact simulationis possible. In Section 5, we apply the proposed ABC method to the stochastic Jansen and Ritneural mass model (JR-NMM). We refer to Jansen and Rit (1995) for the original version, anODE with a stochastic input function, and to Ableidinger et al. (2017) for its reformulation as a3amiltonian type SDE. This model has been reported to successfully reproduce EEG data. Weillustrate the performance of the proposed ABC method with both simulated and real data. Finalremarks, possible extensions and conclusions are reported in Section 6. Further illustrations ofthe proposed ABC method are available in the provided supplementary material, here reportedas Section 7. A sample code used to generate the main results is available on github. Let (Ω , F , P ) be a complete probability space with the right-continuous and complete filtration F = {F} t ∈ [0 ,T ] . Let θ = ( θ , ..., θ k ), k ∈ N , be a vector of relevant model parameters. We considerthe following n -dimensional, n ∈ N , non-autonomous SDE of Itˆo-type describing the time evolutionof a system of interest dX ( t ) = f ( t, X ( t ); θ ) dt + G ( t, X ( t ); θ ) dW ( t ) X (0) = X , t ∈ [0 , T ] . (1)The initial value X is either deterministic or a R n -valued random variable, measurable withrespect to F . Here, W = ( W ( t )) t ∈ [0 ,T ] is a r -dimensional, r ∈ N , Wiener process with independentand F -adapted components. We further assume that the drift component f : [0 , T ] × R n → R n and the diffusion component G : [0 , T ] × R n → R n × r fulfil the necessary global Lipschitz and lineargrowth conditions, such that the existence and the pathwise uniqueness of an F -adapted strongsolution process X = ( X ( t )) t ∈ [0 ,T ] ∈ R n of (1) is guaranteed; see, e.g., Arnold (1974).We aim to infer the parameter vector θ inherent in the SDE (1), when the n -dimensionalsolution process X is only partially observed through the 1-dimensional and parameter-dependentoutput process Y θ = ( Y θ ( t )) t ∈ [0 ,T ] = g ( X ) , (2)where g : R n → R is a real-valued continuous function of the components of X .Further, we assume a specific underlying structural model property, namely the existence of aunique invariant measure η Y θ on ( R , B ( R )) of the output process Y θ , where B denotes the BorelSigma-algebra. The process has invariant density f Y θ and mean, autocovariance and variancegiven by E [ Y θ ( t )] = η µ ∈ R , Cov[ Y θ ( t ) , Y θ ( s )] := r θ ( t, s ) = r θ ( t − s ) , s ≤ t, Var[ Y θ ( t )] = r θ (0) = η σ ∈ R + . (3)If the solution process X of SDE (1) admits an invariant distribution η X on ( R n , B ( R n )), then theoutput process Y θ inherits this structural property by means of the marginal invariant distributionsof η X . Furthermore, if X (0) ∼ η X , then the process Y θ = ( Y θ ( t )) t ∈ [0 , ∞ ) evolves according to thedistribution η Y θ for all t ≥
0. Our goal is to perform statistical inference for the parametervector θ of the SDE (1), when the solution process X is partially observed through discrete timemeasurements of the output process Y θ given in (2), by benefiting from the (in general unknown)invariant distribution η Y θ satisfying (3). Let y = ( y ( t i )) li =1 , l ∈ N , be the reference data, corresponding to discrete time observations ofthe output process Y θ . Let us denote by π ( θ ) and π ( θ | y ) the prior and the posterior density,respectively. For multivariate complex SDEs, the underlying likelihood is often unknown or in-tractable. The idea of the ABC method is to derive an approximate posterior density for θ byreplacing the unknown likelihood by possibly billions of synthetic dataset simulations generatedfrom the underlying model (1) and mapped to Y θ through (2). The basic acceptance-rejection4BC algorithm consists of three steps: i. Sample a value θ (cid:48) from the prior π ( θ ); ii. Condition-ally on θ (cid:48) , simulate a new artificial dataset from the model (1) and derive the synthetic data y θ (cid:48) = ( y θ (cid:48) ( t i )) mi =0 , t = 0 , t m = T, m ∈ N , from the process Y θ (cid:48) given by (2); iii. Keep the sampledparameter value θ (cid:48) as a realisation from the posterior if the distance d between a vector of sum-mary statistics s = ( s , . . . , s h ) , h ∈ N , of the original and the synthetic data is smaller than somethreshold level (cid:15) ≥
0, i.e., d ( s ( y ) , s ( y θ (cid:48) )) < (cid:15) . Algorithm 1
Acceptance-rejection ABC
Input:
Observed data y Output:
Samples from the posterior π ABC ( θ | y ) Precompute a vector of summary statistics s ( y ) Choose a prior distribution π ( θ ) and a tolerance level (cid:15) for i = 1 : N do Draw θ i = ( θ i , ..., θ ik ) from the prior π ( θ ) Conditionally on θ i , simulate a new realisation y θ i from the output process Y θ Compute the summaries s ( y θ i ) Calculate the distance D i = d ( s ( y ) , s ( y θ i )) If D i < (cid:15) , keep θ i as a sample from the posterior end for When (cid:15) = 0 and s is a vector of sufficient statistics for θ , the acceptance-rejection ABC (sum-marised in Algorithm 1) produces samples from the true posterior π ( θ | y ). Due to the complexityof the underlying SDE (1), we cannot derive non-trivial sufficient statistics s for θ . Moreover,due to the underlying stochasticity of the model, P ( d ( s ( y ) , s ( y θ (cid:48) )) = 0) = 0. Thus, (cid:15) is requiredto be strictly positive. Hence, the acceptance-rejection ABC Algorithm 1 yields samples from anapproximated posterior π ABC ( θ | y ) according to π ( θ | y ) ≈ π ABC ( θ | y ) = π { θ | d ( s ( y ) , s ( y θ )) < (cid:15) } . Besides the tolerance level (cid:15) , which we fix as a percentile of the calculated distances, the qualityof the ABC method depends strongly on the choice of suitable summary statistics combined witha proper distance measure and on the numerical method used to generate the synthetic data fromthe model. In the following, we introduce summaries that are very effective for the class of modelshaving an underlying invariant distribution, we suggest a proper distance based on them and wepropose the use of measure-preserving numerical splitting schemes.
When applying ABC to stochastic models, an important statistical challenge arises. Due to theintrinsic randomness, repeated simulations of the process Y θ under the same parameter vector θ may yield very different trajectories. An illustration is given in Figure 1 (top and middle panels),where we report two trajectories of the output process of the stochastic JR-NMM (25) generatedwith an identical parameter configuration. This model is a specific SDE of type (1), observedthrough Y θ as in (2), and admitting an invariant distribution η Y θ satisfying (3). See Section 5for a description of the model. In the top panel, we visualise the full paths for a time T = 200,while in the middle panel we provide a zoom, showing only the initial part. Proposal 1: To use the property of an invariant measure η Y θ and to map the data y θ to their estimated invariant density ˆ f y θ and invariant spectral density ˆ S y θ . Instead of working with the output process Y θ , we take advantage of the structural modelproperty η Y θ and focus on its invariant density f Y θ and its invariant spectral density S Y θ . Bothare deterministic functions characterized by the underlying parameters θ , and thus invariant for5 Y q ( t ) t Y q ( t ) n S^ y q ( n ) y −5 0 5 10 15 20 . . . . . f^ y q ( y ) Fig. 1.
Two realizations of the output process of the stochastic JR-NMM (25) generated withthe numerical splitting method (17) for an identical choice of θ . The lengths of the time intervalsare T = 200 and T = 3 (to provide a zoom) in the top and middle panel, respectively. The twoinvariant densities and two invariant spectral densities, estimated from the two full datasets shownin the top panel, are reported in the lower panel on the left and right, respectivelyrepeated simulations under the same parameter configuration. The invariant spectral density isobtained from the Fourier transformation of the autocovariance function r θ , and it is given by S Y θ = F{ r θ } ( ω ) = (cid:90) ∞−∞ r θ ( τ ) e − iωτ dτ, (4)for ω ∈ [ − π, π ]. The angular frequency ω relates to the ordinary frequency ν via ω = 2 πν . Sinceboth f Y θ and S Y θ are typically unknown, we estimate them from a dataset y θ . First, we estimatethe invariant density f Y θ with a kernel density estimator, denoted by ˆ f y θ ; see, e.g., Pons (2011).Second, we estimate the invariant spectral density S Y θ (4) with a smoothed periodogram estimator(Cadonna et al. 2017; Quinn et al. 2014), denoted by ˆ S y θ , which is typically evaluated at Fourierfrequencies. Differently from the invariant density, the invariant spectral density does not accountfor the mean E [ Y θ ] but captures the dependence structure of the data coming from the model.We define the invariant measure-based summary statistics s of a dataset y θ as s ( y θ ) := ( ˆ S y θ , ˆ f y θ ) . (5)Figure 1 shows the two estimated invariant densities (left lower panel) and invariant spectraldensities (right lower panel), all derived from the full paths of the output process Y θ (top panel).After performing the data mapping (5), which significantly reduces the randomness in theoutput of the stochastic simulator, the distance d can be chosen among the distance measuresbetween two R -valued functions. Here, we consider the integrated absolute error (IAE) defined byIAE( g , g ) := (cid:90) R (cid:12)(cid:12)(cid:12) g ( x ) − g ( x ) (cid:12)(cid:12)(cid:12) dx ∈ R + . (6)6nother natural possibility could be a distance chosen among the so-called f-divergences (see, e.g.,Sason and Verd´u 2016), or the Wasserstein distance, recently proposed for ABC (Bernton et al.2019). Within the ABC framework (see Step 7 in Algorithm 1), we suggest to use the followingdistance d ( s ( y ) , s ( y θ )) := IAE( ˆ S y , ˆ S y θ ) + w · IAE( ˆ f y , ˆ f y θ ) , (7)returning a weighted sum of the areas between the densities estimated from the original and thesynthetic datasets. Here, w ≥ L times: Draw θ (cid:48) from the prior π ( θ ) Conditionally on θ (cid:48) , simulate two artificial datasets y θ (cid:48) and y θ (cid:48) from the output process Y θ Compute the corresponding summaries (5), i.e., s ( y θ (cid:48) ) = ( ˆ S y θ (cid:48) , ˆ f y θ (cid:48) ) and s ( y θ (cid:48) ) = ( ˆ S y θ (cid:48) , ˆ f y θ (cid:48) ) Determine a value for the weight using (7), i.e., w (cid:48) = IAE( ˆ S y θ (cid:48) , ˆ S y θ (cid:48) )IAE( ˆ f y θ (cid:48) , ˆ f y θ (cid:48) ) Then, we take the median of the resulting L values w (cid:48) . See, e.g., Prangle (2017) for alternativeapproaches for the derivation of weights among summary statistics. Since the densities ˆ f y θ andˆ S y θ are estimated at discrete points, the IAE (6) is approximated applying trapezoidal integration.In Algorithm 1, we assume to observe M ∈ N datasets referring to M realisations of the outputprocess Y θ sampled at l ∈ N discrete points in time, resulting in a matrix y ∈ R M × l of observeddata. The median of the distances (7) computed for each of the M datasets D = median (cid:26)(cid:16) IAE( ˆ S y k , ˆ S y θ ) + w · IAE( ˆ f y k , ˆ f y θ ) (cid:17) Mk =1 (cid:27) (8)is then returned as a global distance in Step 7. Other strategies can be adopted. For example,considering the mean instead yields similar results in all our experiments. One can interpret y asa long-time trajectory (when using simulated observed reference data) or as a long-time recordingof the modelled phenomenon (when using real observed reference data) that is cut into M pieces.Alternatively, y would consist of M independent repeated experiments or simulations, when dealingwith real or simulated data, respectively. As expected, having M >
A crucial aspect of ABC and of all other simulation-based methods is the ability of simulating fromthe model (Step 5 of Algorithm 1). Consider a discretized time grid with the equidistant time step∆ = t i +1 − t i and let ˜ y θ = (˜ y θ ( t i )) mi =1 be a realisation from the process (cid:101) Y θ = ( (cid:101) Y θ ( t i )) mi =1 , obtainedthrough a numerical method, approximating Y θ at the discrete data points, i.e., (cid:101) Y θ ( t i ) ≈ Y θ ( t i ).The lack of exact simulation schemes, i.e., (cid:101) Y θ ( t i ) = Y θ ( t i ), introduces a new level of approximationin the statistical inference. In particular, Algorithm 1 samples from an approximated posteriordensity of the form π ( θ | y ) ≈ π numABC ( θ | y ) := π { θ | d ( s ( y ) , s (˜ y θ )) < (cid:15) } . As a consequence, y θ in Step 5 of Algorithm 1 is replaced by its numerical approximation ˜ y θ .The commonly used Euler-Maruyama scheme yields discretised trajectories of the solutionprocess X of the SDE (1) through (Kloeden and Platen 1992) (cid:101) X ( t i +1 ) = (cid:101) X ( t i ) + f ( t i , (cid:101) X ( t i ); θ )∆ + G ( t i , (cid:101) X ( t i ); θ ) ξ i , (9)7 y h Y q ( y ) h Y~ q num ( y ) h Y~ q e ( y ) D =10 - −0.4 −0.2 0.0 0.2 0.4 y h Y q ( y ) h Y~ q num ( y ) h Y~ q e ( y ) D =3*10 - −0.4 −0.2 0.0 0.2 0.4 y h Y q ( y ) h Y~ q num ( y ) h Y~ q e ( y ) D =4.5*10 - Fig. 2.
Comparison of the true invariant density of the weakly damped stochastic harmonicoscillator (23) (blue solid lines) with the densities estimated using a kernel density estimatorapplied on data y θ generated by the measure-preserving splitting scheme (22) (orange dashedlines) and the Euler-Maruyama method (9) (green dotted lines) with time step ∆ up to time T = 10 . The values of the time steps are ∆ = 10 − (left figure), 3 · − (central figure) and4 . · − (right figure), respectivelywhere ξ i are Gaussian vectors with null mean and variance ∆ I n , where I n denotes the n × n -dimensional identity matrix. As previously discussed, in general, the Euler-Maruyama methoddoes not preserve the underlying invariant distribution η Y θ . Proposal 2: To adopt a numerical method for the synthetic data generation thatpreserves the underlying invariant measure of the model.
We apply numerical splitting schemes within the ABC framework and provide a brief accountof their theory. Let us assume that the drift f and the diffusion G of SDE (1) can be written as f ( t, X ( t ); θ ) = d (cid:88) j =1 f [ j ] ( t, X ( t ); θ ) , G ( t, X ( t ); θ ) = d (cid:88) j =1 G [ j ] ( t, X ( t ); θ ) , d ∈ N . The goal is to decompose f and G in a way such that the resulting d subequations dX ( t ) = f [ j ] ( t, X ( t ); θ ) dt + G [ j ] ( t, X ( t ); θ ) dW ( t ) , for j ∈ { , . . . , d } , can be solved exactly. Note that, the terms G [ j ] can be null, resulting indeterministic equations (ODEs). Let X [ j ] ( t ) = ϕ [ j ] t ( X ) denote the exact solutions (flows) of theabove subequations at time t and starting from X . Once these explicit solutions are derived, aproper composition needs to be applied. Here we use the Strang approach (cid:16) ϕ [1]∆ / ◦ ... ◦ ϕ [ d − / ◦ ϕ [ d ]∆ ◦ ϕ [ d − / ◦ ... ◦ ϕ [1]∆ / (cid:17) ( x ) , x ∈ R n , that provides a numerical solution for the original SDE (1).In Figure 2, we illustrate how the numerical splitting method preserves the underlying invariantmeasure of the weakly damped stochastic harmonic oscillator (23), independently from the choiceof the time step ∆. This is a specific SDE of type (1), observed through Y θ as in (2) and with aknown invariant distribution η Y θ . See Section 3 for the detailed numerical splitting scheme andSection 4 for a description of the model. In contrast, the Euler-Maruyama scheme performs worseas ∆ increases. Each subplot shows a comparison of the true invariant density (blue solid lines)and the corresponding kernel estimate ˆ f y θ based on a path y θ from the model, generated from themeasure-preserving numerical splitting scheme (22) (dashed orange lines) or the Euler-Maruyamaapproach (dotted green lines). The data are generated under T = 10 and different values for thetime step, namely ∆ = 10 − , 3 · − , 4 . · − .8 .4 Notation We apply the summary statistics (5) and the distance (8) in Algorithm 1. We use the notationAlgorithm 1 (i) for the Spectral Density-Based ABC method when the synthetic data are simulatedexactly, Algorithm 1 (ii) for the Spectral Density-Based and Measure-Preserving ABC methodwhen a measure-preserving numerical splitting scheme is applied and Algorithm (1) (iii) when wegenerate the data with the non-preserving Euler-Maruyama scheme.To evaluate the performance of the proposed ABC method, we analyse the marginal posteriordensities, denoted by π ∗ ABC ( θ j | y ), j ∈ { , ..., k } , obtained from the posterior density π ∗ ABC ( θ | y )corresponding to π ABC ( θ | y ), π numABC ( θ | y ) or π e ABC ( θ | y ), depending on whether we obtain it fromAlgorithm 1 (i), (ii) or (iii). Following this notation, we define by ˆ θ ∗ ABC ,j the marginal ABCposterior means. We illustrate the proposed ABC approach on Hamiltonian type SDEs and define the n -dimensional( n = 2 d , d ∈ N ) stochastic process X := ( Q , P ) (cid:48) = ( Q ( t ) , P ( t )) (cid:48) t ∈ [0 ,T ] , consisting of the two d -dimensional components Q = ( X , ..., X d ) (cid:48) and P = ( X d + , ..., X ) (cid:48) , where (cid:48) denotes the transpose. The n -dimensional SDE of Hamiltonian type with initial value X = ( Q , P ) (cid:48) and d -dimensional ( r = d ) Wiener process W describes the time evolution of theprocess X by d (cid:18) Q ( t ) P ( t ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) X ( t ) = (cid:18) ∇ P H ( Q ( t ) , P ( t )) −∇ Q H ( Q ( t ) , P ( t )) − θ P ( t ) + G ( Q ( t ); θ ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) f ( X ( t ); θ ) dt + (cid:18) O d Σ θ (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) G ( θ ) dW ( t ) . (10)We denote with O d the d × d -dimensional zero matrix and with ∇ Q and ∇ P the gradient withrespect to Q and P , respectively. The SDE (10) consists of 4 parts, each representing a specific typeof behaviour. In this configuration, the first is the Hamiltonian part involving H : R d × R d → R +0 given by H ( Q , P ) := 12 ( (cid:107) P (cid:107) R d + (cid:107) Λ θ Q (cid:107) R d ) , where Λ θ = diag[ λ , ..., λ d ] ∈ R d × d is a diagonal matrix. The second is the linear damping part ,described by the matrix Γ θ = diag[ γ , ..., γ d ] ∈ R d × d . The third is the non-linear displacementpart , consisting of the non-linear and globally Lipschitz continuous function G : R d → R d . Thefourth corresponds to the diffusion part , given by Σ θ = diag[ σ , ..., σ d ] ∈ R d × d . Under the requirement of non-degenerate matrices Λ θ , Γ θ and Σ θ , i.e., strictly positive diagonalentries, Hamiltonian type SDEs as in (10) are often ergodic. As a consequence, the distribution ofthe solution process X (and thus of the output process Y θ ) converges exponentially fast towardsa unique invariant measure η X on ( R n , B ( R n )) (and thus η Y θ on ( R , B ( R )); see, e.g., Ableidingeret al. (2017) and the references therein. 9 .2 Measure-Preserving numerical splitting schemes Two splitting approaches for SDE (10) are provided, see Ableidinger et al. (2017). Due to thenon-linear term G , the SDE (10) cannot be solved explicitly. With the purpose of excluding G ,the Hamiltonian type SDE (10) is split into the two subsystems d (cid:18) Q ( t ) P ( t ) (cid:19) = (cid:18) ∇ P H (( t ) , P ( t )) −∇ Q H ( Q ( t ) , P ( t )) − θ P ( t ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) f [1] ( X ( t ); θ ) dt + (cid:18) O d Σ θ (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) G [1] ( θ ) dW ( t ) , (11) d (cid:18) Q ( t ) P ( t ) (cid:19) = (cid:18) d G ( Q ( t ); θ ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) f [2] ( Q ( t ); θ ) dt, (12)where 0 d denotes the d -dimensional zero vector. This results in the linear SDE with addi-tive noise (11) and the non-linear ODE (12) that can be both explicitly solved. Indeed, since ∇ P H ( Q ( t ) , P ( t )) = P ( t ) and ∇ Q H ( Q ( t ) , P ( t )) = Λ θ Q ( t ), Subsystem (11) can be rewritten as dX ( t ) = A · X ( t ) dt + B dW ( t ) , t ≥ , (13)with A = (cid:18) O d I d − Λ θ − θ (cid:19) and B = (cid:18) O d Σ θ (cid:19) . The exact path of System (13) is obtained through X ( t i +1 ) = e A ∆ · X ( t i ) + ξ i , (14)where ξ i are d -dimensional Gaussian vectors with null mean and variance C (∆), where the matrix C ( t ) follows the dynamics of the matrix-valued ODE˙ C ( t ) = AC ( t ) + C ( t ) A (cid:48) + BB (cid:48) , (15)see Arnold (1974). Moreover, since the non-linear term G depends only on the component Q , theexact path of Subsystem (12) is obtained through X ( t i +1 ) = X ( t i ) + (cid:18) d ∆ G ( Q ( t i ); θ ) (cid:19) . (16)We apply the Strang approach given by( ϕ b ∆ / ◦ ϕ a ∆ ◦ ϕ b ∆ / )( x ) , x ∈ R n , (17)where ϕ at and ϕ bt denote the exact solutions (14) and (16) of (11) and (12), respectively. Hence,given X ( t i ), we obtain the next value X ( t i +1 ) by applying the following three steps: X b = X ( t i ) + (cid:18) d ∆2 G ( Q ( t i ); θ ) (cid:19) X a = e A ∆ · X b + ξ i X ( t i +1 ) = X a + (cid:18) d ∆2 G ( Q a ; θ ) (cid:19) The derivation of the two subsystems is not unique. For example, another possibility is tocombine the stochastic term with the non-linear part, yielding the subsystems d (cid:18) Q ( t ) P ( t ) (cid:19) = (cid:18) ∇ P H ( Q ( t ) , P ( t )) −∇ Q H ( Q ( t ) , P ( t )) − θ P ( t ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) f [1] ( X ( t ); θ ) dt, (18) d (cid:18) Q ( t ) P ( t ) (cid:19) = (cid:18) d G ( Q ( t ); θ ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) f [2] ( Q ( t ); θ ) dt + (cid:18) O d Σ θ (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) G [2] ( θ ) dW ( t ) . (19)10he exact path of (18) is given by X ( t i +1 ) = e A ∆ · X ( t i ) , (20)while the exact path of (19) is obtained through X ( t i +1 ) = (cid:18) Q ( t i ) P ( t i ) + ∆ G ( Q ( t i ); θ ) + Σ θ · ξ i (cid:19) , (21)where ξ i are d -dimensional Gaussian vectors with null mean and variance ∆ I d . The Strang ap-proach is now given by ( ϕ c ∆ / ◦ ϕ d ∆ ◦ ϕ c ∆ / )( x ) , x ∈ R n , (22)where ϕ ct and ϕ dt denote the exact solutions (20) and (21) of (18) and (19), respectively. Thus,given X ( t i ), the next value X ( t i +1 ) is obtained via: X c = e A ∆2 · X ( t i ) X d = X c + (cid:18) d ∆ G ( Q c ; θ ) + Σ θ · ξ i (cid:19) X ( t i +1 ) = e A ∆2 · X d The ABC procedure is coded in the computing environment R (R Development Core Team 2011),using the package Rcpp (Eddelbuettel and Fran¸cois 2011), which offers a seamless integration of R and C++ , drastically reducing the computational time of the algorithms. The code is thenparallelised using the R-packages foreach and doParallel , taking advantage of the for-loop inthe algorithm. All simulations are run on the HPC cluster
RADON1 , a high-performing multiplecore cluster located at the Johannes Kepler University Linz. To obtain smoothed periodogram es-timates, we apply the R-function spectrum . It requires the specification of a smoothing parameter span . In all our experiments, we use span = 5 T . In addition, we avoid using a logarithmic scaleby setting the log parameter to “no”. To obtain kernel estimates of the invariant density, we applythe R-function density . Here, we use the default value for the smoothing bandwidth bw and setthe number of points at which the invariant density has to be estimated to n = 10 . The invariantspectral density is estimated at the default values of the spectrum function. A sample code ispublicly available on github at https://github.com/massimilianotamborrino/sdbmpABC . In this section, we illustrate the performance of the proposed ABC approach on a model problem(weakly damped stochastic harmonic oscillator) of Hamiltonian type (10) with vanishing non-linear displacement term G ≡
0. Linear SDEs of this type reduce to (13) and allow for an exactsimulation of sample paths through (14). Therefore, we can apply the Spectral Density-Based ABCAlgorithm 1 (i) under the optimal condition of exact, and thus η Y θ -preserving data generation. Itsperformance is illustrated in Subsection 4.2. To investigate how the numerical error in the syntheticdata generation impinges on the ABC performance, in Subsection 4.3 we compare π ABC ( θ | y ) withthe posterior densities π numABC ( θ | y ) and π eABC ( θ | y ) obtained from Algorithm 1 (ii) and (iii) usingthe measure-preserving numerical splitting scheme (22) and the non-preserving Euler-Maryuamamethod (9), respectively. 11 .1 Weakly damped stochastic harmonic oscillator: The model and itsproperties We investigate the 2-dimensional Hamiltonian type SDE d (cid:16) Q ( t ) P ( t ) (cid:17) = (cid:16) P ( t ) − λ Q ( t ) − γP ( t ) (cid:17) dt + (cid:16) σ (cid:17) dW ( t ) , (23)with strictly positive parameters γ , λ and σ . Depending on the choice of γ and λ , (23) modelsdifferent types of harmonic oscillators, which are common in nature and of great interest in classicalmechanics. Here, we focus on the weakly damped harmonic oscillator, satisfying the condition λ − γ >
0. Our goal is to estimate θ = ( λ, γ, σ ) assuming that the solution process X = ( Q , P ) (cid:48) is partially observed through the first coordinate, i.e., Y θ = Q . An illustration of the performanceof Algorithm 1 (i) for the critically damped case satisfying λ − γ = 0, when only the secondcoordinate is observed, is reported in the supplementary material. The solution process X of SDE(23) is normally distributed according to X ( t ) ∼ η X ( t ) := N (cid:16) e At · E [ X ] , Var[ e At · X ] + C ( t ) (cid:17) , with A and C introduced in (13) and (15), respectively. The invariant distribution η X of thesolution process X is given by η X = lim t →∞ η X ( t ) = N (cid:32)(cid:18) (cid:19) , (cid:32) σ γλ σ γ (cid:33)(cid:33) . Consequently, the structural property η Y θ of the output process Y θ becomes η Y θ = N (cid:18) , σ γλ (cid:19) , (24)and the stationary dependency is captured by the autocovariance function r θ (∆) = σ λ e − γ ∆ (cid:20) γ cos( κ ∆) + 1 κ sin( κ ∆) (cid:21) , where κ = (cid:112) λ − γ . To compare the performances of Algorithm 1 (i)-(iii) on the same data, we consider the same M = 10 observed paths simulated with the exact scheme (14), using a time step ∆ = 10 − over atime interval of length T = 10 . As true parameters for the simulation of the reference data, wechoose θ t = ( λ t , γ t , σ t ) = (20 , , . We use the exact simulation scheme (14) to generate N = 2 · synthetic datasets in [0 , T ]and with the same time step as the observed data. We choose independent uniform priors, inparticular, λ ∼ U (18 , , γ ∼ U (0 . , . , σ ∼ U (1 , . The tolerance level (cid:15) is chosen as the 0 . th percentile of the calculated distances. Hence, we keep10 of all the sampled values for θ . In all the considered examples (see also the supplementarymaterial), the performance of the ABC algorithms for the estimation of the parameters of SDE (23)does not improve when incorporating the information of the invariant densities into the distance(7). This is because the mean of the invariant distribution (24) is zero. Hence, to reduce thecomputational cost, we set w = 0 and base our distance only on the invariant spectral density,estimated by the periodogram. 12igure 3 (top panels) shows the marginal ABC posterior densities π ABC ( θ j | y ) (blue lines) andtheir flat uniform priors π ( θ j ) (red lines). The proposed ABC Algorithm 1 (i) provides marginalposterior densities centred around the true values θ t , represented by the black vertical lines. Theposterior means are given by(ˆ λ ABC , ˆ γ ABC , ˆ σ ABC ) = (20 . , . , . . In the lower panels of Figure 3, we report the pairwise scatterplots of the kept ABC posteriorsamples. Note that, since the kept values of λ are uncorrelated with those of the other parame-ters, the support of the obtained marginal posterior density is approximately the same as whenestimating only θ = λ or θ = ( λ, γ ) (cf. supplementary material). Vice versa, since the kept ABCposterior samples of the parameters γ and σ are correlated, the support of π ABC ( γ | y ) is largerthan that obtained when estimating θ = ( λ, γ ). Despite this correlation, Algorithm 1 (i) allowsfor a successful inference of all the three parameters. . . . . . . l
18 19 20 21 22 p ( l ) p ABC ( l |y ) . . . . . . . g p ( g ) p ABC ( g |y ) . . . . . . . s p ( s ) p ABC ( s |y ) ll lll ll l ll l ll l lll l llll l llll lll lll ll ll lll ll ll l ll ll l lll l llll l ll l llll ll ll ll l lll ll ll llll ll lll llll lll l ll ll l lll ll lll lll lll ll ll lllll l l lll ll ll l l l lll ll lll l llll l lll l ll ll ll lll l llll lll l ll llll l ll llllll lll lll ll l ll l l ll l llll l lll l llll ll l lll l lll lll ll lll llll ll ll l llll llll l lll lll ll ll llll ll ll lll lllll ll l lllll l lll l lll ll lll ll l llll ll lll llll ll ll l ll l llll l ll ll l ll llll ll ll ll ll lll lll l ll ll lll ll l lllll ll lll ll l ll l ll lll ll lll l ll l ll l ll lll l ll ll ll lll l llll llll lll ll lll ll lllllll l l lll ll l l llll llll lllll ll ll lll lll l ll l lll l llll ll lll ll ll l lll l lll lll llll ll ll ll l ll l ll lll lll ll l lll ll lll ll llll ll ll l lll l ll l l l ll lll ll lll ll ll ll llll l llll lll ll lllll ll l lll ll llll l lll llll ll lllll l l lll ll ll ll ll ll lllll lll ll ll ll l ll l ll llll lll ll ll ll ll l ll ll lll ll l l lll l l l lll l ll ll ll ll ll ll ll lll ll l ll llll l l ll ll lll l l lll l lll llll l ll lll lll l ll l l ll l lll l llll lll ll llll l lll ll ll ll l ll lll ll l ll ll l lllll lllll ll ll ll llll l l l ll lll llll l lll l l ll lll lllll l l ll lll l lll ll l llll l ll lll l ll lll ll l l ll ll l ll ll lll l l llll l llll ll l l ll lll l l ll l ll ll l ll ll l ll ll ll ll l ll ll ll l ll lll llllll lll l ll l lll l ll ll lll lll ll llll llll l l llll ll ll ll llll l . . . . . . . . l g ll lll ll l ll l ll l lll l llll l llll lll lll ll ll lll ll ll l ll ll l lll l llll l ll l llll ll ll ll l lll ll ll llll ll lll llll lll l ll ll l lll ll lll lll lll ll ll lllll l l lll ll ll l l l lll ll lll l llll l lll l ll ll ll lll l llll lll l ll llll l ll llllll lll lll ll l ll l l ll l llll l lll l llll ll l lll l lll lll ll lll llll ll ll l llll llll l lll lll ll ll llll ll ll lll lllll ll l lllll l lll l lll ll lll ll l llll ll lll llll ll ll l ll l llll l ll ll l ll llll ll ll ll ll lll lll l ll ll lll ll l lllll ll lll ll l ll l ll lll ll lll l ll l ll l ll lll l ll ll ll lll l llll llll lll ll lll ll lllllll l l lll ll l l llll llll lllll ll ll lll lll l ll l lll l llll ll lll ll ll l lll l lll lll llll ll ll ll l ll l ll lll lll ll l lll ll lll ll llll ll ll l lll l ll l l l ll lll ll lll ll ll ll llll l llll lll ll lllll ll l lll ll llll l lll llll ll lllll l l lll ll ll ll ll ll lllll lll ll ll ll l ll l ll llll lll ll ll ll ll l ll ll lll ll l l lll l l l lll l ll ll ll ll ll ll ll lll ll l ll llll l l ll ll lll l l lll l lll llll l ll lll lll l ll l l ll l lll l llll lll ll llll l lll ll ll ll l ll lll ll l ll ll l lllll lllll ll ll ll llll l l l ll lll llll l lll l l ll lll lllll l l ll lll l lll ll l llll l ll lll l ll lll ll l l ll ll l ll ll lll l l llll l llll ll l l ll lll l l ll l ll ll l ll ll l ll ll ll ll l ll ll ll l ll lll llllll lll l ll l lll l ll ll lll lll ll llll llll l l llll ll ll ll llll l . . . . . . . l s lll ll lll lll llll ll llll ll lll ll ll llll ll l ll ll ll l lllll lll ll llll lll lll ll l ll l l ll l lllll lll l lll ll l llll lll lllll l ll lll l lll lll llll ll lll lll l lll ll l ll l lllllll ll ll lll llll ll lll l ll l l lll ll ll ll ll l lll ll l ll lll l ll l lllll llll lll llll l l lll l lll l lllll ll lll lll lll lll lll ll lll l lllll ll ll ll lll ll ll l llll ll ll l lll lll l lll llll ll ll l ll lllll l ll ll lll ll l ll llll ll ll lll l l lll l l ll l lll lll l ll ll lll l l lllll l ll ll ll lll ll l llllllll ll ll l ll l ll ll llll ll l ll ll l ll lll ll ll ll ll l lllll llll l ll lll ll ll ll ll llll l l ll ll ll lll llll l lll ll ll ll l l lll lll ll ll lll ll lll l ll ll ll ll ll ll l llll ll l llll ll ll lll l l lll ll lll l lll ll ll ll l l lll llll llll lll ll lll ll l lllll ll lll lll lll ll ll l ll l l llll ll l ll lll l llllll ll ll l l lll ll ll l lll lll l lll ll ll llll ll ll ll lll ll lll ll llll ll l llll ll l lll l l ll lll ll l lll ll lll llll l lll l lll l l lll l lll l lll ll lll lll ll l ll l lll ll l lll l lll lll l ll l ll l lll l l ll ll ll ll l ll ll ll ll lll l ll lll ll lllll ll lll lll ll lll l lll lll ll ll ll ll ll ll l ll l ll l ll lll ll ll ll l ll llll l ll l ll l llll ll lll l l l ll llll l ll ll l ll ll lll ll ll lll l l lll llll l ll ll llll l l ll lll ll lll ll l ll ll ll l lll l ll ll lllll lll llll lll l l ll ll llll l lll ll l lll l l llll l ll lll l lll ll llll l l ll l . . . . . . . g s Fig. 3.
Top panels: ABC marginal posterior densities π ABC ( θ j | y ) (blue lines) of θ = ( λ, γ, σ )of the weakly damped stochastic harmonic oscillator (23) and uniform priors (red lines). Theposteriors are obtained from Algorithm 1 (i). The vertical lines represent the true parametervalues. Lower panels: Pairwise scatterplots of the kept ABC posterior samples In Figure 4, we report the approximated marginal posteriors π ABC ( θ j | y ) (blue solid lines) and π numABC ( θ j | y ) (orange dashed lines) obtained with the same priors, (cid:15) , T , w , M and N as before,for different values of the time step ∆. In particular, we choose ∆ = 5 · − (top panels),∆ = 7 . · − (middle panels) and ∆ = 10 − (lower panels). The posteriors obtained fromAlgorithm 1 (ii) successfully targets π ABC ( θ | y ), even for a time step as large as ∆ = 10 − . On thecontrary, Algorithm 1 (iii) is not even applicable. Indeed, the numerical scheme computationallypushes the amplitude of the oscillator towards infinity, resulting in a computer overflow, i.e., (cid:101) Y θ ( t i ) ≈ ∞ . Thus, neither ˆ f ˜ y θ nor ˆ S ˜ y θ can be computed and the density π e ABC ( θ | y ) cannot bederived. 13 . . . . . . l
18 19 20 21 22 p ( l ) p ABC ( l |y ) p ABCnum ( l |y ) . . . . . . . g p ( g ) p ABC ( g |y ) p ABCnum ( g |y ) D =5*10 - . . . . . . . s p ( s ) p ABC ( s |y ) p ABCnum ( s |y ) . . . . . . l
18 19 20 21 22 p ( l ) p ABC ( l |y ) p ABCnum ( l |y ) . . . . . . . g p ( g ) p ABC ( g |y ) p ABCnum ( g |y ) D =7.5*10 - . . . . . . . s p ( s ) p ABC ( s |y ) p ABCnum ( s |y ) . . . . . . l
18 19 20 21 22 p ( l ) p ABC ( l |y ) p ABCnum ( l |y ) . . . . . . . g p ( g ) p ABC ( g |y ) p ABCnum ( g |y ) D =10 - . . . . . . . s p ( s ) p ABC ( s |y ) p ABCnum ( s |y ) Fig. 4.
ABC marginal posterior densities of θ = ( λ, γ, σ ) of the weakly damped stochasticharmonic oscillator (23) obtained from Algorithm 1 (i) with the exact simulation method (14)(blue solid lines) and Algorithm 1 (ii) combined with the splitting scheme (22) (orange dashedlines) for different choices of the time step ∆. In particular, ∆ = 5 · − (top panels), 7 . · − (middle panels) and 10 − (lower panels). The red horizontal lines denote the uniform priors andthe black vertical lines the true parameter valuesAs a further illustration of the poor performance of the Euler-Maruyama scheme, even forsmaller choices of ∆, we now consider the simplest possible scenario where we only estimate oneparameter, namely θ = λ . We set N = 10 , M = 10, (cid:15) = 1 st percentile and we choose a uniformprior λ ∼ U (10 , π e ABC ( λ | y ), we simulate the synthetic data using theEuler-Maruyama method with the time steps ∆ = 10 − , 2 . · − and 3 . · − . Figure 5 showsthe three ABC posterior densities π ABC ( θ | y ) (blue solid lines), π numABC ( θ | y ) (orange dashed lines)and π e ABC ( θ | y ) (green dotted lines) for the different choices of ∆. The horizontal red lines andthe black vertical lines denote the uniform prior and the true parameter value, respectively. In allcases, Algorithm 1 (iii) does not lead to a successful inference. In addition, these results are notstable for the different choices of ∆, and the derived ABC posterior density may not even coverthe true parameter value. 14 . . . . . . l
18 19 20 21 22 D =10 - p ( l ) p ABC ( l |y ) p ABCnum ( l |y ) p ABCe ( l |y ) . . . . . . l
18 19 20 21 22 D =2.5*10 - p ( l ) p ABC ( l |y ) p ABCnum ( l |y ) p ABCe ( l |y ) . . . . . . l
18 19 20 21 22 D =3.5*10 - p ( l ) p ABC ( l |y ) p ABCnum ( l |y ) p ABCe ( l |y ) Fig. 5.
ABC posterior densities of θ = λ of the weakly damped stochastic oscillator (23) obtainedfrom Algorithm 1 (i) using the exact simulation scheme (14) (blue solid lines), (ii) using thesplitting scheme (22) (orange dashed lines) and (iii) using the Euler-Maruyama method (9) (greendotted lines) for different choices of the time step ∆. The horizontal red lines and the verticalblack lines represent the uniform priors and the true parameter values, respectively We now illustrate the performance of Algorithm 1 (ii) by applying it to the stochastic JR-NMM.We rely on the efficient numerical splitting scheme (17) to guarantee measure-preserving syntheticdata generation within the ABC framework. After estimating the parameters from simulated data,we infer them from real EEG data. In the available supplementary material, we illustrate theperformance of Algorithm 1 (ii) also on the non-linear damped stochastic oscillator, an extendedversion of the weakly damped harmonic oscillator discussed in Section 4.
The stochastic JR-NMM describes the electrical activity of an entire population of neurons throughtheir average properties by modelling the interaction of the main pyramidal cells with the sur-rounding excitatory and inhibitory interneurons. The model has been reported to successfullyreproduce EEG data, and is applied in the research of neurological disorders such as epilepsy orschizophrenia (Wendling et al. 2000, 2002). The model is a 6-dimensional SDE of the form d (cid:18) Q ( t ) P ( t ) (cid:19) = (cid:18) P ( t ) − Γ Q ( t ) − P ( t ) + G ( Q ( t ); θ ) (cid:19) dt + (cid:18) θ (cid:19) dW ( t ) , (25)where the 6-dimensional solution process is given by X = ( Q , P ) (cid:48) with the two components Q = ( X , X , X ) (cid:48) and P = ( X , X , X ) (cid:48) . None of the coordinates of X is directly observed.Only the difference between the second and third coordinates can be measured with EEG-recordingtechniques, yielding the output process Y θ = X − X . In (25), the diagonal diffusion matrix is given by Σ θ =diag[ σ , σ , σ ] ∈ R × with coefficients σ i > i = 4 , ,
6. The matrix Γ=diag[ a, a, b ] ∈ R × is also diagonal with coefficients a, b >
0, rep-resenting the time constants of the excitatory and inhibitory postsynaptic potentials, respectively.15he non-linear displacement term is given by G ( Q ; θ ) = Aa [Sigm( X − X )] Aa [ µ + C Sigm( C X )] Bb [ C Sigm( C X )] , where the sigmoid function Sigm: R → [0 , v max ] is defined asSigm( x ) := v max r ( v − x )] , with v max > v ∈ R describingthe value for which 50 % of the maximum firing rate is attained and r > v . The parameters entering in G are µ , A , B and C i , i = 1 , , , ∈ R + . Thecoefficients A and B describe the average excitatory and inhibitory synaptic gain, respectively.The parameters C i are internal connectivity constants, which reduce to only one parameter C , byusing the relations C = C , C = 0 . C , C = 0 . C and C = 0 . C ; see Jansen and Rit (1995). Not all model parameters of the JR-NMM are of biological interest or can be simultaneouslyidentified. For example, the noise coefficients σ and σ were introduced mainly for mathematicalconvenience in Ableidinger et al. (2017). To guarantee the existence of a unique invariant measure η X on ( R , B ( R )), they are required to be strictly positive. However, from a modelling pointof view, only the parameter σ := σ plays a role. Hence, we fix σ = 0 .
01 and σ = 1. Thecoefficients A , B , a , b , v , v max and r have been experimentally recorded; see, e.g., Jansen et al.(1993); Jansen and Rit (1995); van Rotterdam et al. (1982). Thus, we fix them according to thesevalues reported, for example, in Table 1 of Ableidinger et al. (2017). In contrast, the connectivityparameter C , which represents the average number of synapses between the neural subpopulationsand controls to what extent the main population interacts with the interneurons, varies underdifferent physiological constraints. Changing C allows, for example, a transition from α -rhythmicactivity to epileptic spiking behaviour; see, e.g., Ableidinger et al. (2017). Here, we focus on the α -rhythmic activity. Since the parameters σ and µ are new in the SDE version (25), they havenot yet been estimated. They can be interpreted as stochastic and deterministic external inputscoming from neighbouring or more distant cortical columns, respectively. Thus, together with theinternal connectivity parameter C , they are of specific interest. Before inferring θ = ( σ, µ, C ), wetake into account the coefficients A and B to discuss a model-specific issue of identifiability. For the original JR-NMM, it has been shown that different combinations of the parameters A , B and C yield the same type of output, namely the α -rhythmic brain activity. Applying theproposed Spectral Density-Based and Measure-Preserving ABC Algorithm 1 (ii) for the inferenceof θ = ( A, B, C ), with given µ = 220 and σ = 2000, we confirm that the same non identifiabilityarises for the SDE version (25). We choose M = 30 observed paths generated assuming θ t = ( A t , B t , C t ) = (3 . , , , as suggested in the literature (Jansen and Rit 1995). The reference and synthetic data are gen-erated over a time interval of length T = 200 and using a time step ∆ = 2 · − . Within thealgorithm, we generate N = 2 . · synthetic datasets. We choose the weight w in (7) accordingto the procedure introduced in Subsection 2.2 (based on L = 10 iterations) and fix the tolerancelevel (cid:15) = 0 . th percentile. Further, we choose independent uniform prior distributions, namely A ∼ U (1 , , B ∼ U (10 , , C ∼ U (10 , . . . . . . . A p ( A ) p ABCnum ( A|y ) B
10 30 50 70 90 . . . . . p ( B ) p ABCnum ( B|y ) C
10 140 270 400 530 . . . p ( C ) p ABCnum ( C|y ) lll ll lll l llll l ll l l llll ll ll ll lll ll llllll llll lll ll ll l llll lll lllll l lll llll llll l llll l l lll lll ll lll ll l l ll l lll lll lll l ll ll ll ll l lll ll ll l lll ll ll ll ll ll l ll llll l l ll l lll l lll lllllll l lll ll l lll l l ll ll l ll lll ll ll ll l l lll ll lll l ll l ll l lllll l lll ll ll llll ll lll ll lll ll ll l lll l ll ll ll ll lll l ll l llll lll l llll l ll ll ll lll lll l l ll l lllll lll llll l ll ll l ll ll llll ll llll lll lll ll l lll l lll ll lll lll ll l l ll lll l ll ll ll ll l lll lll l l ll l ll llll ll l l l lll l lll ll ll ll l ll l ll lll l ll lll lll l lllll ll l lll ll lll l ll l llll l ll l ll l lll ll ll lllll l llll ll l l l ll ll ll lll lll ll lll l ll ll ll llll lll l lll lll lll l ll lll llll l ll lll l ll lll l ll ll ll ll l ll ll llll llll ll ll ll ll ll l l ll lll lll lll l l ll lll l lll lll ll l lll ll lll lll ll lll l ll ll l llll lll l lll lllll ll l ll ll llll ll l ll lll lll l ll lllll ll lll ll lll llll l ll ll lll l ll llll l lll ll ll ll l lll l ll ll llll l ll l ll ll lll l l ll l llll lll l l l lll lll ll lllll ll ll ll ll llll l ll l llll lll ll lll ll ll ll ll ll l ll ll lll llll lll l ll ll ll l ll l lll l ll ll l l lll lll lll ll lllll l ll ll ll ll l lll ll ll l lll l ll ll lll lll lllll l llll ll llll l lll lllll ll ll lll l lll ll llll l lll lll lllll ll lll l ll ll ll lll ll lll l lllllll lll ll ll lll lll ll lll ll ll ll ll ll ll llll lll l A ll ll B lll ll lll l llll l ll l l llll ll ll ll lll ll llllll llll lll ll ll l llll lll lllll l lll llll llll l llll l l lll lll ll lll ll l l ll l lll lll lll l ll ll ll ll l lll ll ll l lll ll ll ll ll ll l ll llll l l ll l lll l lll lllllll l lll ll l lll l l ll ll l ll lll ll ll ll l l lll ll lll l ll l ll l lllll l lll ll ll llll ll lll ll lll ll ll l lll l ll ll ll ll lll l ll l llll lll l llll l ll ll ll lll lll l l ll l lllll lll llll l ll ll l ll ll llll ll llll lll lll ll l lll l lll ll lll lll ll l l ll lll l ll ll ll ll l lll lll l l ll l ll llll ll l l l lll l lll ll ll ll l ll l ll lll l ll lll lll l lllll ll l lll ll lll l ll l llll l ll l ll l lll ll ll lllll l llll ll l l l ll ll ll lll lll ll lll l ll ll ll llll lll l lll lll lll l ll lll llll l ll lll l ll lll l ll ll ll ll l ll ll llll llll ll ll ll ll ll l l ll lll lll lll l l ll lll l lll lll ll l lll ll lll lll ll lll l ll ll l llll lll l lll lllll ll l ll ll llll ll l ll lll lll l ll lllll ll lll ll lll llll l ll ll lll l ll llll l lll ll ll ll l lll l ll ll llll l ll l ll ll lll l l ll l llll lll l l l lll lll ll lllll ll ll ll ll llll l ll l llll lll ll lll ll ll ll ll ll l ll ll lll llll lll l ll ll ll l ll l lll l ll ll l l lll lll lll ll lllll l ll ll ll ll l lll ll ll l lll l ll ll lll lll lllll l llll ll llll l lll lllll ll ll lll l lll ll llll l lll lll lllll ll lll l ll ll ll lll ll lll l lllllll lll ll ll lll lll ll lll ll ll ll ll ll ll llll lll l A ll ll C lll ll lll l llll l ll l l llll ll ll ll l ll ll lll lll llll ll l ll ll l llll lll lll ll l lll ll ll l lll l llll l l lll lll ll lll ll l l ll l lll lll lll l ll ll ll ll l lll ll ll l lll ll ll l l ll ll l ll llll l l ll l lll l lll lllllll l lll ll l l ll l l ll ll l ll lll ll ll ll l l lll ll lll l ll l ll l lllll l lll ll ll llll ll l ll l l lll ll ll l lll l ll llll ll lll l lll llll lll l l lll l ll ll ll lll lll l l ll l l llll l ll ll ll l ll ll l ll ll llll ll llll lll lll ll l lll l l ll ll lll lll ll l l ll lll l ll ll ll ll l l ll lll ll ll l ll llll ll l l l llll lll ll ll ll l ll l ll lll l ll l ll lll l lllll ll l lll ll lll l ll l llll l ll l ll l lll ll ll lllll l ll ll ll l l l ll ll ll lll lll ll lll l ll ll ll llll lll l ll l lll lll l ll lll llll l ll ll l l ll lll l ll llll ll l ll ll llll llll ll l l ll ll ll l l ll lll lll lll l l ll lll l lll llll l l ll l ll lll lll ll lll l ll ll l ll ll l ll l lll l llll ll l ll ll l lll ll l ll lll lll l ll lll ll ll lll ll lll llll l ll ll lll l ll l lll l lll ll ll ll l lll l ll ll llll l ll l ll ll lll l l lll llll lll l l l lll lll ll lllll ll ll ll l lll ll l ll l l lll lll ll lll ll ll ll ll ll lll ll l ll llll lll l ll l l ll l ll l lll l ll ll l l lll lll l ll ll lllll l ll ll ll ll l lll ll ll l lll l ll ll lll lll l ll ll l llll ll llll l lll lllll ll ll lll l lll ll llll l lll l ll ll lll ll lll l ll ll ll lll ll lll l lllllll ll l ll lllll lll ll lll ll ll ll ll llll llll lll l B ll ll
10 30 50 70 90 C Fig. 6.
Top panels: ABC marginal posterior densities π numABC ( θ j | y ) (blue lines) of θ = ( A, B, C )of the stochastic JR-NMM (25) obtained from Algorithm 1 (ii). The horizontal red lines andthe vertical black lines represent the uniform priors and the true parameter values, respectively.Middle panels: Pairwise scatterplots of the kept ABC posterior samples. Lower panels: Twodifferent views of a 3-dimensional scatterplot of the kept ABC posterior samples within a cuboidformed by the prior. The green dot corresponds to θ t and the red, orange and grey dots representhighlighted samples from the ABC posterior lying on the invariant manifoldFigure 6 (top panels) shows the marginal ABC posterior densities π numABC ( θ j | y ) and the uni-form prior densities π ( θ j ). Clearly, the parameters cannot be inferred simultaneously. The keptABC posterior values of the parameters A , B and C are strongly correlated, as observed in thepairwise scatterplots (middle panels) and in the 3-dimensional scatterplot (two different views,17ower panels). The cuboid covers all possible values for θ drawn from the prior. After running theABC algorithm, the kept values of θ from the ABC posterior form an invariant manifold, in thesense that all the parameter values θ lying on this manifold yield similar paths ˜ y θ of the outputprocess. This is shown in Figure 7, where we report four trajectories that have been simulatedwith the same random numbers but using the parameters θ t (green dot in Figure 6) and threeof the kept ABC posterior samples lying on the invariant manifold (red, orange and grey dots inFigure 6). A segment of T = 10 is split in the top and middle panels. In addition, we visualise thecorresponding estimated invariant densities (bottom left) and invariant spectral densities (bottomright). This explains why the parameters A , B and C are not simultaneously identifiable from theobserved data. Since the internal connectivity parameter C has an important neuronal meaning,in the following we assume A and B to be known and infer θ = ( σ, µ, C ). The estimation of θ = ( σ, µ ) when C is known is reported in the supplementary material. t Y q ( t ) t Y q ( t ) n S^ y q ( n ) y −5 0 5 10 15 20 . . . . . f^ y q ( y ) Fig. 7.
Top and middle panel: Four paths of the output process Y θ = X − X of the stochasticJR-NMM (25) generated under θ t (green lines) and with the three highlighted kept ABC posteriorsamples lying on the invariant manifold of Figure 6 (red, orange and grey lines) using the samerandom numbers. Lower panels: Corresponding estimated invariant densities (left) and estimatedspectral densities (right) θ = ( σ, µ, C )Now, we keep the same ABC setting as before and choose independent uniform priors π ( θ j )according to σ ∼ U (1300 , , µ ∼ U (160 , , C ∼ U (129 , . The reference data are simulated under θ t = ( σ t , µ t , C t ) = (2000 , , . In Figure 8, we report the marginal ABC posterior densities π numABC ( θ j | y ) (blue lines), the uniformprior densities π ( θ j ) (red lines) and the true parameter values θ t (black vertical lines). We obtain18nimodal posterior densities, centred around the true parameter values. The posterior density of σ is slightly broader compared to that obtained when C is known (cf. Figure 19 of the supplementarymaterial). This results from a weak correlation that we detect among the kept ABC posteriorsamples of the parameters σ and C (figures not reported). The posterior means are equal to(ˆ σ numABC , ˆ µ num ABC , ˆ C numABC ) = (1992 . , . , . , and are thus close to θ t . These results suggest an excellent performance of the proposed SpectralDensity-Based and Measure-Preserving ABC Algorithm 1 (ii). . . . . . . s p ( s ) p ABCnum ( s |y ) m
160 180 200 220 240 260 280 . . . . . p ( m ) p ABCnum ( m |y ) . . . . . C
129 131 133 135 137 139 141 p ( C ) p ABCnum ( C|y ) Fig. 8.
ABC marginal posterior densities π numABC ( θ j | y ) (blue lines) of θ = ( σ, µ, C ) of the stochasticJR-NMM (25) obtained from Algorithm 1 (ii). The horizontal red lines and the vertical blacklines represent the uniform priors and the true parameter values, respectivelySimilar satisfactory results are obtained even when adding a fourth parameter, for example,when inferring θ = ( σ, µ, C, b ) (cf. Figure 20 of the supplementary material). When applyingAlgorithm 1 (ii) to real EEG data (cf. Figure 21 of the supplementary material), the marginalposterior for b is centred around the value b = 50, which is that reported in the literature. Dueto the existence of underlying invariant manifolds, identifiability issues, similar to those reportedin Figure 6, arise when adding further or other coefficients, revealing model-specific issues for thestochastic JR-NMM.To illustrate again the importance of the structure-preservation within the ABC method, wenow apply Algorithm 1 (iii) combined with the Euler-Maruyama scheme (9). We use the sameconditions as before, except for a smaller time step ∆ = 10 − used for the generation of theobserved reference data with the Euler-Maruyama method aiming for a realistic data structure.In Figure 9, we report the marginal ABC posterior densities π e ABC ( θ j | y ) (top panels) and theuniform prior densities. In the 3-dimensional scatterplot of Figure 9 (lower panel), the green dotsin the middle of the cuboid represent the kept ABC posterior samples when applying Algorithm1 (ii) (see the previous results reported in Figure 8), which are nicely spread-out around the trueparameter vector θ t (black dot). The red dots correspond to the kept ABC posterior samples from π e ABC ( θ | y ). Hence, Algorithm 1 (iii) based on the Euler-Maruyama scheme provides a posteriorthat is far off from the true parameter vector. 19 . . . . p ( s ) p ABCEM ( s |y ) . . . . . . . m
160 180 200 220 240 260 280 p ( m ) p ABCEM ( m |y ) C
129 131 133 135 137 139 141 p ( C ) p ABCEM ( C|y ) Fig. 9.
Top panels: Marginal ABC posterior densities π e ABC ( θ j | y ) (blue lines) of θ = ( σ, µ, C )of the stochastic JR-NMM (25) obtained from Algorithm 1 (iii) using the non-preservative Euler-Maruyama scheme (9). The horizontal red lines and the vertical black lines represent the uniformpriors and the true parameter values, respectively. Lower panel: 3-dimensional scatterplot of thekept ABC posterior samples using Algorithm 1 (ii) (green dots; see the previous results reportedin Figure 8) and Algorithm 1 (iii) (red dots). The cuboid is formed by the prior. The black dotcorresponds to θ t Finally, we use the Spectral Density-Based and Measure-Preserving ABC Algorithm 1 (ii) toestimate the parameter vector θ = ( σ, µ, C ) of the stochastic JR-NMM from real EEG recordings.We use M = 3 α -rhythmic recordings, rescaled to a realistic range. The EEG data were sampledaccording to a sampling rate of 173.61 Hz, i.e., a time step ∆ of approximately 5 .
76 ms overa time interval of length T = 23 . .Figure 10 shows the first 20 seconds of one of the observed EEG datasets. Here, we simulate N = 5 · synthetic paths from the output process of the stochastic JR-NMM (25) over the sametime interval T as the real data, with a time step ∆ = 2 · − and (cid:15) = 0 . nd percentile. Wechoose independent uniform priors π ( θ j ) according to σ ∼ U (500 , , µ ∼ U (70 , , C ∼ U (120 , . The data are available on: http://ntsa.upf.edu/downloads/andrzejak-rg-et-al-2001-indications-nonlinear-deterministic-and-finite-dimensional .0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.05.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.010.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.015.0 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0 Fig. 10.
Visualisation of the first 20 seconds of one of the used α -rhythmic EEG segmentsrecorded with a sampling rate of 173 .
61 Hz, i.e., ∆ ≈ .
76 msFigure 11 shows the resulting marginal ABC posterior densities π numABC ( θ j | y ) and the uniformprior densities π ( θ j ). All ABC marginal posteriors are unimodal, with means given by(ˆ σ numABC , ˆ µ num ABC , ˆ C numABC ) = (1859 . , . , . . Since µ and σ have not been estimated before, we cannot compare the obtained results with thoseavailable in the literature. The ABC posterior density for C is centred around C = 135 that isthe reference literature value for α -rhythmic EEG data. . . . . s
700 1300 1900 2500 3100 p ( s ) p ABCnum ( s |y ) m
70 130 190 250 310 370 . . . . p ( m ) p ABCnum ( m |y ) . . . . . . C
120 126 132 138 144 150 p ( C ) p ABCnum ( C|y ) Fig. 11.
Marginal ABC posterior densities π numABC ( θ j | y ) (blue lines) of θ = ( σ, µ, C ) of the stochasticJR-NMM (25) fitted on real EEG data using Algorithm 1 (ii). The red lines correspond to theuniform priorsIn Figure 12, we report the first 10 seconds of a trajectory of the output process of the fittedstochastic JR-NMM (25), generated with the numerical splitting scheme (17), choosing ∆ = 2 · − and T = 23 .
6. Note how the path shows a similar oscillatory behaviour as in Figure 10. This isconfirmed by noting the satisfactory matches between the invariant densities (bottom left) and theinvariant spectral densities (bottom right) estimated from the EEG recording (red dashed lines)21nd from the fitted model (blue solid lines). The match is poor only for low frequencies of theinvariant spectral density, even when choosing broader priors. This may result from a lack of fit ofthe JR-NMM or of stationarity in the considered EEG data. A deeper investigation of the modeland of its ability in reproducing real EEG data is currently under investigation, but it is out ofthe scope of this work. t Y q ^ ( t ) t Y q ^ ( t ) n S^ y q ( n ) JR - NMMEEGy −5 0 5 10 15 20 . . . . . f^ y q ( y ) JR - NMMEEG
Fig. 12.
Top and middle panel: Visualisation of 10 s of a sample path of the stochastic JR-NMM(25) generated with the numerical splitting scheme (17) using (ˆ σ numABC , ˆ µ num ABC , ˆ C numABC ), the marginalABC posterior means derived from Algorithm 1 (ii). The chosen time step is ∆ = 2 · − and T = 23 . s . Lower panels: Corresponding estimated invariant density (solid blue line, left) andinvariant spectral density (solid blue line, right) plotted against those estimated from the realEEG dataset shown in Figure 10 (red dashed lines) When performing parameter inference through ABC, crucial and non-trivial tasks are to proposesuitable summary statistics and distances to compare the observed and the synthetic datasets.When the underlying models are stochastic, repeated simulations from the same parameter set-ting yield different outputs, making the comparison between the observed and the synthetic datamore difficult. To derive summary statistics that are less sensitive to the intrinsic randomness ofthe stochastic model, we propose to map the data to their invariant density and invariant spec-tral density, estimated by a kernel density estimator and a smoothed periodogram, respectively.By doing this, different trajectories of the output process are mapped to the same objects onlywhen they are generated from the same underlying parameters, provided that all parameters aresimultaneously identifiable. These transformations are based on the existence of an underlyinginvariant measure for the model, fully characterised by the parameters. A necessary condition ofABC, and of all other simulation-based methods, is the ability to generate data from the model.This is often taken for granted but, in general, it is not the case. Indeed, exact simulation is rarelypossible and property-preserving numerical methods have to be derived.22he combination of the measure-preserving numerical splitting schemes and the use of thespectral density-based distances in the ABC algorithm leads to a successful inference of the pa-rameters, as illustrated on stochastic Hamiltonian type equations. We validated the proposedABC approach on both linear model problems, allowing for an exact simulation of the syntheticdata, and non-linear problems, including an application to real EEG data. Our choice of the cru-cial ingredients (summary statistics and distances based on the underlying invariant distributionand a measure-preserving numerical method) yields excellent results even when applied to ABCin its basic acceptance-rejection form. However, they can be directly applied to more advancedABC algorithms. In contrast, the ABC method based on the Euler-Maruayma scheme drasticallyfails. Its performance may improve for “small enough” time steps. However, there is a trade-offbetween the runtime and the acceptance performance of Algorithm 1 (iii). Indeed, the simula-tion of one trajectory with a time step 10 − requires approximately hundred times more thanthe generation of one trajectory using a time step 10 − . Hence, a runtime of a few hours wouldturn to months. In addition, even for “arbitrary small” time steps, one cannot guarantee thatthe Euler-Maruyama scheme preserves the underlying invariant measure. For these reasons, it iscrucial to base our ABC method on the reliable measure-preserving numerical splitting schemecombined with the invariant measure-based distances. Our results were discussed in the case ofan observable 1-dimensional output process. However, the approach can be directly applied to d -dimensional output processes, d >
1, as long as the underlying SDEs are characterised by aninvariant distribution and a measure-preserving numerical method can be derived. In particular,one can compute the distances in (8) for each of the d components and derive a global distanceby combining them, e.g., via their sum. Moreover, to account for possible dependences betweenthe observed components, one can incorporate the cross-spectral densities which are expected toprovide further information resulting in an improvement of the performance of the method. Aninvestigation in this direction is currently undergoing. Finally, our proposed ABC method may bealso used to investigate invariant manifolds characterised by sets of parameters yielding the sametype of data, as illustrated on the stochastic JR-NMM. This may result in a better understandingof the qualitative behaviour of the underlying model and its ability of reproducing the true featuresof the modelled phenomenon. In this supplementary material, we extend the illustration of the performance of the proposedABC approach by more examples. In particular, we consider two additional SDEs. First, thecritically damped harmonic oscillator (23), fulfilling λ − γ = 0, for which an exact simulationof sample paths is possible, allowing for a validation of Algorithm 1 (i). Second, a non-linearweakly damped stochastic oscillator, for which we need to apply a measure-preserving numericalsplitting scheme, and thus investigate Algorithm 1 (ii). Moreover, we also report the simultaneousinference of the new parameters θ = ( σ, µ ) in the stochastic JR-NMM (25), when the connectivityparameter C is known (while in the main manuscript, we estimate θ = ( σ, µ, C )). Finally, wereport the estimation of θ = ( σ, µ, C, b ), based on both simulated and real EEG data. We denote by Model Problem 1 (MP1) the critically damped harmonic oscillator obtained from(23) with λ − γ = 0 (introduced below), and with Model Problem 2 (MP2) the weakly dampedharmonic oscillator, satisfying λ − γ > - . - . . . . Y q ( t ) Fig. 13.
Two paths of the output process Y θ = P of the critically damped stochastic harmonicoscillator (MP1) for a noise intensity σ = 2 and parameter γ = 1 t Y q ( t ) Fig. 14.
Two paths of the output process Y θ = Q of the weakly damped stochastic harmonicoscillator (MP2) for a noise intensity σ = 2, γ = 1 and a damping force λ = 20 We recall the harmonic oscillator (23), focusing on the critically damped case, i.e., λ − γ = 0.We assume that the 2-dimensional process X = ( Q , P ) (cid:48) is partially observed through the secondcomponent, i.e., Y θ = P . The invariant distribution η X of the process X is given by η X = lim t →∞ η X ( t ) = N (cid:32)(cid:18) (cid:19) , (cid:32) σ γ σ γ (cid:33)(cid:33) . Consequently, the invariant distribution η Y θ of the output process Y θ equals η Y θ = N (cid:18) , σ γ (cid:19) , and the autocovariance function is given by r θ (∆) = σ e − γ ∆ (cid:20) γ − ∆ (cid:21) . At first, we estimate one specific parameter θ of the model problems, keeping the others fixed.For MP1, we set θ = γ and fix σ = 2. For MP2, we focus on θ = λ , fixing γ = 1 and σ = 2. TheABC Algorithm 1 (i) is applied to both model problems with M = 10 observed paths simulatedwith the exact scheme (14) using a time step ∆ = 10 − over a time interval of length T = 10 .In addition, we generate N = 10 synthetic datasets over the same time domain with equal timesteps using the exact simulation scheme (14). We set the tolerance level to (cid:15) = 1 st percentile of24he calculated distances. Furthermore, we choose uniform prior distributions π ( θ ) according to θ = γ ∼ U (0 . ,
5) for MP1 ,t tλ ∼ U (10 ,
30) for MP2 . We use the same parameter setting as in Figure 13 and Figure 14 for the simulation of the observedreference datasets. In particular, the true parameter values are θ t = γ t = 1 for MP1 ,t tλ t = 20 for MP2 . In Figure 15, we report the results obtained from the proposed Spectral Density-Based ABCAlgorithm 1 (i). The left panel (referring to MP1) and the right panel (referring to MP2) showthe ABC posterior densities π ABC ( θ | y ) (blue lines). The horizontal red and vertical black linesdenote the prior densities and the true parameter values, respectively. It is remarkable howthe flat uniform prior densities are updated by means of the observed data resulting in narrowand unimodal posterior densities that are centered around the true parameter values. The ABCposterior means for this and the other scenarios (i.e. the inference of two and three parameters)are reported in Table 1. g p ( g ) p ABC ( g |y ) . . . . . . l
18 19 20 21 22 p ( l ) p ABC ( l |y ) Fig. 15.
ABC posterior densities π ABC ( θ | y ) (blue lines) of MP1 (left panel) and MP2 (right panel)obtained from Algorithm 1 (i). The horizontal red and vertical black lines denote the uniformpriors (not shown according to the full domain) and the true parameter values, respectively We aim for the simultaneous estimation of two parameters, keeping the parameter σ = 2 fixed inMP2. In particular, we consider θ = ( γ, σ ) for MP1 and θ = ( λ, γ ) for MP2. We apply Algorithm1 (i) combined with the exact scheme (14) under the same values for M , ∆ and T as before. Now,we generate N = 5 · synthetic datasets and fix (cid:15) = 0 . nd percentile of the calculated distances,keeping the same number of ABC posterior samples as before. We choose the independent uniformpriors π ( θ j ) according to θ j = γ ∼ U (0 . , .
01) and σ ∼ U (1 ,
3) for MP1 ,t tλ ∼ U (18 ,
22) and γ ∼ U (0 . , .
01) for MP2 . θ t = ( γ t , σ t ) = (1 ,
2) for MP1 ,t t ( λ t , γ t ) = (20 ,
1) for MP2 . The ABC marginal posterior densities π ABC ( θ j | y ) (blue lines) are reported in the left andmiddle panels of Figure 16 for MP1 (top panels) and MP2 (lower panels), while the right panelsof Figure 16 show the scatterplots of the kept ABC posterior samples. Also in this case, theposteriors are unimodal and centered around the true parameter values. Note that, the supportof π ABC ( λ | y ) for MP2 is approximately the same as in Figure 15, suggesting that, in the case ofinferring two parameters, the proposed ABC method is able to identify the same region for λ as inthe case of estimating one parameter. The reason is that the kept ABC posterior samples of λ and γ are not correlated, as it can be observed in the right lower panel of Figure 16. On the contrary,the support of π ABC ( γ | y ) for MP1 is broader than in Figure 15, due to a correlation among thekept ABC posterior samples of γ and σ (cf. right top panel of Figure 16). In spite of this, theABC marginal posterior density resembles that derived when estimating only one parameter (cf.left panel of Figure 15). g p ( g ) p ABC ( g |y ) s p ( s ) p ABC ( s |y ) lll lll ll llll lll ll ll lll lll lllll llll ll llll ll ll l lll llll lll l lll ll lll ll ll ll ll ll l ll lll lll lll ll ll l lll l lll lllll l ll lll ll lll ll l ll l ll l ll ll llllll ll lll ll ll lll lllllllll l ll l ll lll l ll ll ll llllll lll ll ll l lll ll l l ll l l ll l llll ll ll l ll lllll llll ll l l l lll ll lll ll l lll ll l l ll llll l l l ll ll lll l ll llll ll lll l ll l llll l ll ll l ll ll ll l ll ll lll ll ll ll ll l ll llllll ll l lll lll ll llll l l ll l ll l l llll ll ll ll l lll lll lll ll ll ll l ll lll l l ll ll lll ll ll ll l lll ll ll ll l ll l lll l l ll l l lll ll l l lll ll l lll l lll ll lll ll ll ll ll l lll l ll ll ll lll ll l lll lll l ll ll ll ll ll l lll ll ll ll lllll ll l llll llll ll llll llll ll ll ll l l ll lll lll ll lll ll l ll l lll lllll ll lll ll ll ll lll l lll ll l ll ll l ll l l l ll ll l ll llll l ll ll llll ll ll ll l ll ll l ll ll ll ll llll l l ll l ll llll ll ll l ll l ll ll ll ll llll ll ll ll lll ll l lllll l lll ll l lll l ll l ll ll lll ll lllll ll l l ll llll lll l llll l ll lll l l lll l ll l ll lll lll ll lll ll l lllllll lll l ll ll l ll ll l l ll lll lll l l l ll l lll lllllll lll l lll ll ll ll l lll ll lll ll ll lll ll ll ll ll l llll l ll lll l l lll ll l lll ll ll ll l lll l lll ll ll lll llll ll lll l lll l ll l lllll l lll lll lll lll ll ll ll ll ll l llll l lll ll lll lll lllll ll l ll lll ll llll l lllll l llll l lll ll ll ll lll l l ll ll l l ll ll ll l llll ll l l g . . . . s . . . . . . l
18 19 20 21 22 p ( l ) p ABC ( l |y ) g p ( g ) p ABC ( g |y ) llll l llllll lll lll l ll ll ll l lll ll llllll llllll llll ll lllll ll ll lll llll l ll lll ll lll l ll lll lll llll l l ll ll l lll lll llll l ll l ll lll lll lll lll l llll ll lll lllll lll lll llll llll ll llll ll lll lll llll l lll llll ll lll ll lllll ll ll ll ll ll ll lll lll l llll ll l lll l llll ll ll lll lllll ll l lll lll lll lll l llll l l llllll l l ll ll lll ll l lll lll l l l ll lll ll llll l llll lll llll l ll llll ll lll l ll l lll ll ll llll l lll ll lll lll ll ll l ll lll ll ll llll ll lllll l ll ll ll lll l ll ll l ll lll lll ll lll ll l lll l ll llll l lll ll l l llll ll ll ll l ll ll lllll l ll l ll lll ll llll ll ll lll ll ll ll llll l ll ll ll ll ll lll l lll lll lll ll l lll lll l ll ll llllll llllll lll ll ll ll llll lll ll l lllll l l llll lll l l ll ll llll l lll l lll llll ll ll lll lll lll ll l lll ll ll ll ll ll lll ll ll llll ll l l ll l lll l ll ll l ll l ll ll llll llll ll l llll l ll l lll l ll l lll ll ll llll ll lll l llll lllll ll l lll l lllll l ll l lll l l lll ll lllll l lll l lll ll ll lll ll ll ll lll l ll llll ll ll llll l l llll lll ll ll l ll l ll l l ll ll l l ll ll l ll ll ll l llll l ll ll l lll ll ll l ll ll ll lll l lll l ll lll l lll l l lll lll ll ll ll lll lll lll lll lll ll llll ll ll llll lll lll l lll lll lll ll l ll lll lll ll ll ll l lll ll ll ll ll ll l lll ll l l lll ll l lllll lll llll lll ll ll l l l ll lll lll lll lll lllll ll ll ll l l lll . . . . . . l g Fig. 16.
ABC marginal posterior densities π ABC ( θ j | y ) (blue lines) of MP1 (left and middle toppanels) and MP2 (left and middle lower panels) obtained from Algorithm 1 (i). The horizontalred and vertical black lines denote the uniform priors and the true parameter values, respectively.Scatterplots of the kept ABC posterior samples for MP1 and MP2 are reported in the right topand right lower panel, respectively The last goal is the simultaneous inference of all the three parameters θ = ( λ, γ, σ ) of MP2. InFigure 3 (top panels) we report the ABC marginal posterior densities (blue lines) and the priordensities (red lines). In the lower panels, we show the pairwise scatterplots of the kept ABCposterior samples. The kept posterior values of λ turned out to be not correlated with those of Task 3 is already presented in Subsection 4.2 of the main manuscript. For completeness, we recall it here. γ and σ are correlated (cf. right lowerpanel of Figure 3), leading to a support for γ broader than that in Figure 16. The ABC marginalposterior densities shown in Figures 3, 15 and 16, and the results reported in Table 1 highlight thegood performance of the proposed Spectral Density-Based ABC Algorithm 1 (i) under the optimalcondition of exact, and thus measure-preservative data simulation from the underlying model.Table 1: Parameters of interest, true parameter values and ABC posterior means θ θ t ˆ θ ABC
MP1 γ . γ, σ ) (1 ,
2) (0 . , . λ
20 20 . λ, γ ) (20 ,
1) (20 . , . λ, γ, σ ) (20 , ,
2) (20 . , . , . We now consider an extended non-linear version of the previously studied Model Problem 2. Dueto the non-linearity in the model, an exact simulation scheme is not available. Hence, we considerthe measure-preserving numerical splitting scheme (17), and thus investigate the performance ofAlgorithm 1 (ii).
We consider a stochastic oscillator that incorporates a high-amplitude sine wave represented bythe non-linear displacement term G ( Q ) = − sin( Q ). In particular, we study the 2-dimensionalSDE d (cid:18) Q ( t ) P ( t ) (cid:19) = (cid:18) P ( t ) − λ Q ( t ) − γP ( t ) + G ( Q ( t )) (cid:19) dt + (cid:18) σ (cid:19) dW ( t ) (26)with the strictly positive parameters θ = ( λ, γ, σ ). The condition λ − γ > X = ( Q , P ) (cid:48) is partially observed through the firstcoordinate, i.e., Y θ = Q . Figure 17 shows two realisations of the output process generated withthe same choice of parameters. t - . - . . . Y q ( t ) Fig. 17.
Two paths of the output process Y θ = Q of the non-linear stochastic oscillator (26) for θ = ( λ, γ, σ ) = (20 , ,
2) 27 .2.2 Parameter inference from simulated data
We assume to observe M = 30 paths of the output process simulated with the measure-preservingnumerical scheme (17) over a time interval of length T = 10 using a time step ∆ = 10 − and thesame true parameter values as in Figure 17, i.e., θ t = ( λ t , γ t , σ t ) = (20 , , . We then use the same T and ∆ to generate N = 2 · synthetic datasets within ABC. We furtherchoose the tolerance level (cid:15) = 0 . th percentile of the calculated distances, set w = 0 in (7) anduse independent uniform prior distributions λ ∼ U (18 , , γ ∼ U (0 . , .
01) and σ ∼ U (1 , . Figure 18 shows the ABC marginal posterior densities π numABC ( θ j | y ). They are unimodal, narrowand centered around the true parameter values. The ABC posterior means are given by(ˆ λ numABC , ˆ γ num ABC , ˆ σ num ABC ) = (20 . , . , . . In spite of the presence of the non-linear term G , the inference via Algorithm 1 (ii) yields resultssimilar to those obtained for MP2 when applying Algorithm 1 (i) under the exact data generation. . . . . l
18 19 20 21 22 p ( l ) p ABCnum ( l |y ) g p ( g ) p ABCnum ( g |y ) s p ( s ) p ABCnum ( s |y ) Fig. 18.
ABC marginal posterior densities π numABC ( θ j | y ) (blue lines) of θ = ( λ, γ, σ ) of the non-linear weakly damped stochastic oscillator (26) obtained from Algorithm 1 (ii). The horizontalred and vertical black lines denote the uniform priors and the true parameter values, respectively θ =( σ, µ ) of the stochastic JR-NMM We now estimate θ = ( σ, µ ) of the stochastic JR-NMM (25) (see Section 5 of the main manuscript).These are new parameters introduced by Ableidinger et al. (2017) in the SDE reformulation of theoriginal JR-NMM (Jansen and Rit 1995). Differently from the other parameters, these parametershave not yet been estimated in the literature. Here, we fix C = 135 and apply Algorithm 1 (ii)with M = 30, N = 5 · , ∆ = 2 · − and T = 200. We fix (cid:15) = 0 . nd percentile of the calculateddistances and choose uniform prior distributions according to σ ∼ U (1300 , µ ∼ U (160 , . The true parameter values used to generate the observed data are given by θ t = ( σ t , µ t ) = (2000 , . Figure 19 shows the ABC marginal posterior densities π numABC ( θ j | y ) (left and middle top panels)obtained by applying the Spectral Density-Based and Measure-Preserving ABC Algorithm 1 (ii).28he posteriors are centered around the true parameter values, leading to marginal ABC posteriormeans given by (ˆ σ numABC , ˆ µ numABC ) = (1985 . , . . From the scatterplot of the kept ABC posterior samples of σ and µ (right top panel), we concludethat they are not correlated. The successful performance of the proposed ABC approach is alsovisible by looking at the contour plot of the ABC posterior density (lower panel). Indeed, theproposed algorithm is able to detect a plain region of posterior values for θ around θ t . . . . . . s p ( s ) p ABCnum ( s |y ) . . . . . . . m
160 180 200 220 240 260 280 p ( m ) p ABCnum ( m |y ) ll l ll lll llll ll ll ll lll ll ll l ll llll l lll ll ll ll lll ll lll l ll lll l lll ll ll ll l ll lll llll lll lll lll ll lllll lll ll ll l ll llll lll lll ll ll lll l l llll ll l ll lll llll lllll ll lll ll lll llll llll lllll ll ll l ll llll l ll ll lll l lll l l llll ll llll lll lll ll lll ll ll l ll ll llll l ll l lll l ll ll ll lll l ll ll llll l llll ll lll ll ll l ll lll lll ll lll lll l lll ll ll ll lll l lll lll llll ll ll ll l ll ll l ll l ll ll lll lll l lll ll ll llll lll ll llllllll llll ll lll ll lll ll ll llll ll l ll ll lll l ll llll ll lllll l ll l lll l lll ll ll llll lll ll ll ll l ll l ll ll l ll ll ll lll l ll l llll l ll ll lll ll l llll lll l ll l l ll ll ll l ll l lll lll lll lll l l ll ll l l llllll ll lllll llll llll l llll l ll ll ll llll llll ll l ll l llll ll ll l lllll lll ll ll ll lll ll l ll lll l ll ll ll l lll l ll l lll ll lll llll lll llll ll llll ll ll ll lll ll l l ll ll l ll l ll ll ll lll l lll ll l llll l lll lll ll llll ll l lll ll ll ll ll lll lll ll l ll l ll lll llll ll ll l lll ll lll llll lll lll lll lll l l ll llll ll lll l ll ll l lll lll lll ll ll ll ll ll l ll ll lllll ll l ll l llll ll llll l l ll llllll lll l lll l l llll ll l l ll lll ll llll ll ll llll ll ll l llll ll llll llll ll l lll lll l lllll l lll ll l llll ll lll l ll ll l lll ll l lll ll lll l l lll ll ll lllll l ll llll ll lllll ll llll l llll ll llll llll l lll ll llll ll ll l l lll l ll l l s m sig_PD − − − − − − − − . sm Fig. 19.
ABC marginal posterior densities π numABC ( θ j | y ) (blue lines, left and middle top panels)of θ = ( σ, µ ) of the stochastic JR-NMM (25) obtained from Algorithm 1 (ii). The horizontal redand vertical black lines denote the uniform priors and the true parameter values, respectively.Scatterplot of the kept ABC posterior samples (right top panel) and contour plot of the ABCposterior density (lower panel) θ = ( σ, µ, C, b ) of thestochastic JR-NMM We now demonstrate that we obtain satisfactory results even when inferring the four parameters θ = ( σ, µ, C, b ) of the stochastic JR-NMM (25). Since the parameters of main interest are σ , µ and C , in the main manuscript (see Section 5) we did not take into account the well-reportedcoefficient b , which takes the value b = 50 in the literature; see, e.g., Jansen and Rit (1995) andthe references therein. 29 .4.1 Inference from simulated data We start with inferring θ = ( σ, µ, C, b ) from simulated data and apply Algorithm 1 (ii) for M = 30, N = 5 · , ∆ = 2 · − and T = 200. We fix (cid:15) = 0 . th percentile and use the following uniformpriors σ ∼ U (1300 , , µ ∼ U (160 , ,C ∼ U (129 , b ∼ U (44 , . The reference data is generated under θ t = ( σ t , µ t , C t , b t ) = (2000 , , , . s . . . . p ( s ) p ABCnum ( s |y ) m
160 180 200 220 240 260 280 . . . . . p ( m ) p ABCnum ( m |y ) . . . . . C
129 131 133 135 137 139 141 p ( C ) p ABCnum ( C|y ) . . . . . . . . b
44 46 48 50 52 54 56 p ( b ) p ABCnum ( b|y ) Fig. 20.
ABC marginal posterior densities π numABC ( θ j | y ) (blue lines) of θ = ( σ, µ, C, b ) of thestochastic JR-NMM (25) obtained from Algorithm (1) (ii). The horizontal red lines and thevertical black lines represent the uniform priors and the true parameter values, respectivelyIn Figure 20, we report the marginal ABC posterior densities π numABC ( θ j | y ), which are againcentered around the true parameter values. The marginal posterior means are given by(ˆ σ numABC , ˆ µ numABC , ˆ C numABC , ˆ b numABC ) = (1992 . , . , . , . . .4.2 Inference from real EEG data Finally, we infer θ = ( σ, µ, C, b ) from real EEG data. Algorithm 1 (ii) is applied under the sameconditions as in Subsection 5.3 of the main manuscript, except fixing (cid:15) = 0 . nd percentile ofcalculated distances and choosing the uniform priors according to σ ∼ U (500 , , µ ∼ U (70 , ,C ∼ U (120 , b ∼ U (40 , . . . . . s
700 1300 1900 2500 3100 p ( s ) p ABCnum ( s |y ) m
70 130 190 250 310 370 . . . . p ( m ) p ABCnum ( m |y ) . . . . . C
120 126 132 138 144 150 p ( C ) p ABCnum ( C|y ) . . . . . b
40 45 50 55 60 p ( b ) p ABCnum ( b|y ) Fig. 21.
Marginal ABC posterior densities π numABC ( θ j | y ) (blue lines) of θ = ( σ, µ, C, b ) of thestochastic JR-NMM (25) fitted on real EEG data using Algorithm (1) (ii). The red lines correspondto the uniform priorsFigure 21 shows the unimodal marginal ABC posterior densities π numABC ( θ j | y ), yielding posteriormeans given by (ˆ σ numABC , ˆ µ numABC , ˆ C numABC , ˆ b numABC )= (1902 . , . , . , . . Focusing on the coefficient b , the corresponding marginal posterior density is centered around b = 50, which is the value reported in the literature.31 eferences Ableidinger, M., Buckwar, E.: Splitting Integrators for the Stochastic Landau–Lifshitz Equation.SIAM J. Sci. Comput. 38, A1788–A1806 (2016)Ableidinger, M., Buckwar, E., Hinterleitner, H.: A Stochastic Version of the Jansen and Rit NeuralMass Model: Analysis and Numerics. J. Math. Neurosci. 7(8) (2017)Andrzejak, R. G., Lehnertz, K., Mormann, F., Rieke, C., David, P., Elger, C. E.: Indicationsof nonlinear deterministic and finite-dimensional structures in time series of brain electricalactivity: Dependence on recording region and brain state. Phys. Rev. E 64, 061907 (2001)Arnold, L.: Stochastic differential equations: theory and applications. Wiley, New York (1974)Barber, S., Voss, J., Webster, M.: The rate of convergence for approximate Bayesian computation.Electron. J. Stat. 9(1), 80–105 (2015)Barnes, C., Filippi, S., Stumpf, M., Thorne, T.: Considerate approaches to constructing summarystatistics for ABC model selection. Stat. Comput. 22(6), 1181–1197 (2012)Beaumont, M. A., Zhang, W., Balding, D. J.: Approximate Bayesian Computation in PopulationGenetics. Genetics 162(4), 2025–2035 (2002)Bernton, E., Jacob, P. E., Gerber, M., Robert, C. P.: Approximate Bayesian computation withthe Wasserstein distance. J. Roy. Stat. Soc. B (2019)Biau, G., C´erou, F., Guyader, A.: New Insights Into Approximate Bayesian Computation. Ann.I. H. Poincare B 51(1), 376–403 (2015)Blanes, S., Casas, F., Murua, A.: Splitting and composition methods in the numerical integrationof differential equations. Bol. Soc. Esp. Mat. Apl. 45 (2009)Blum, M. G. B.: Approximate Bayesian Computation: A Nonparametric Perspective. J. Am. Stat.Assoc. 105(491), 1178–1187 (2010a)Blum, M. G. B.: Choosing the Summary Statistics and the Acceptance Rate in ApproximateBayesian Computation. In: Lechevallier, Y., Saporta, G. (eds) Proceedings of COMPSTAT2010: Physica-Verlag HD, Heidelberg: pp 47–56 (2010b)Boys, R. J., Wilkinson, D. J., Kirkwood, T. B. L.: Bayesian inference for a discretely observedstochastic kinetic model. Stat. Comput. 18, 125–135 (2008)Br´ehier, C. E., Gouden`ege, L.: Analysis of Some Splitting Schemes for the Stochastic Allen-CahnEquation. Discrete Cont. Dyn.-B 24, 4169–4190 (2019)Cadonna, A., Kottas, A., Prado, R.: Bayesian mixture modeling for spectral density estimation.Stat. Prob. Lett. 125, 189–195 (2017)Drovandi, C. C., Pettitt, A. N., McCutchan, R.: Exact and Approximate Bayesian Inferencefor Low Integer-Valued Time Series Models with Intractable Likelihoods. Bayesian Anal. 11,325–352 (2016)Eddelbuettel, D., Fran¸cois, R.: Rcpp: Seamless R and C++ Integration. J. Stat. Soft. 40(8), 1–18(2011)Fan, Y., Sisson, S. A. (2018) ABC samplers. In: Sisson, S. A., Fan, Y., Beaumont, M. (eds)Handbook of Approximate Bayesian Computation: CRC Press, Taylor & Francis Group: chap 4,pp 87–123 32earnhead, P., Prangle, D.: Constructing summary statistics for approximate Bayesian computa-tion: semi-automatic approximate Bayesian computation. J. Roy. Stat. Soc. B 74(3), 419–474(2012)Jansen, B. H., Rit, V. G.: Electroencephalogram and visual evoked potential generation in amathematical model of coupled cortical columns. Biol. Cybern. 73(4), 357–366 (1995)Jansen, B. H., Zouridakis, G., Brandt, M. E.: A neurophysiologically-based mathematical modelof flash visual evoked potentials. Biol. Cybern. 68, 275–283 (1993)Jasra, A.: Approximate Bayesian Computation for a Class of Time Series Models. Int. Stat. Rev.83(3), 405–435 (2015)Jiang, B., Wu, T.-y., Zheng, C., Wong, W. H.: Learning summary statistics for ApproximateBayesian Computation via deep neural network. Stat. Sinica 27(4), 1595–1618 (2017)Kloeden, P. E., Platen, E.: Numerical Solution of Stochastic Differential Equations. Springer,Berlin (1992)Kypraios, T., Neal, P., Prangle, D.: A tutorial introduction to Bayesian inference for stochasticepidemic models using Approximate Bayesian Computation. Math. Biosci. 287, 42–53 (2017)Leimkuhler, B., Matthews, C.: Molecular dynamics: with deterministic and stochastic numericalmethods. Springer International Publ., Cham (2015)Leimkuhler, B., Matthews, C., Stoltz, G.: The computation of averages from equilibrium andnonequilibrium Langevin molecular dynamics. IMA. J. Numer. Anal. 36(1), 16–79 (2016)Lintusaari, J., Gutmann, M., Dutta, R., Kaski, S., Corander, J.: Fundamentals and RecentDevelopments in Approximate Bayesian Computation. Syst. Biol. 66(1), e66–e82 (2017)Malham, S. J., Wiese, A.: Chi-square simulation of the CIR process and the Heston model. Int.J. Theor. Appl. Finance 16(3) (2013)Marin, J.-M., Pudlo, P., Robert, C. P., Ryder, R.: Approximate Bayesian computational methods.Stat. Comput. 22(6), 1167–1180 (2012)Martin, G. M., McCabe, B. P. M., Frazier, D. T., M., W., Robert, C. P.: Auxiliary Likelihood-Based Approximate Bayesian Computation in State Space Models. J. Comput. Graph. Stat. 0(0), 1–31 (2019)Mattingly, J. C., Stuart, A. M., Higham, D. J.: Ergodicity for SDEs and approximations: locallyLipschitz vector fields and degenerate noise. Stoch. Proc. Appl. 101(2), 185–232 (2002)Maybank, P., Bojak, I., Everitt, R.: Fast approximate Bayesian inference for stable differentialequation models (2017) https://arxiv.org/abs/1706.00689
McKinley, T. J., Vernon, I., Andrianakis, I., McCreesh, N., Oakley, J., Nsubuga, R., Goldstein, M.,White, R.: Approximate Bayesian Computation and Simulation-Based Inference for ComplexStochastic Epidemic Models. Stat. Sci. 33(1), 4–18 (2017)Mclachlan, R., Quispel, G.: Splitting methods. Acta Numer. 11, 341–434 (2002)Milstein, G. N., Tretyakov, M. V.: Stochastic numerics for mathematical physics. Scientific com-putation: Springer, Berlin (2004)Misawa, T.: A Lie Algebraic Approach to Numerical Integration of Stochastic Differential Equa-tions. SIAM J. Sci. Comput. 23(3), 866–890 (2001)Moores, M. T., Drovandi, C. C., Mengersen, K., Robert, C. P.: Pre-processing for approximateBayesian computation in image analysis. Stat. Comput. 25, 23–33 (2015)33ori, U., Mediburu, A., Lozano, J. A.: Distance measures for time series in R: The TSdist package.R. J. 8, 455–463 (2016)Moro, E., Schurz, H.: Boundary Preserving Semianalytic Numerical Algorithms for StochasticDifferential Equations. SIAM J. Sci. Comput. 29, 1525–1549 (2007)Muskulus, M., Verduyn-Lunel, S.: Wasserstein distances in the analysis of time series and dynam-ical systems. Physica. D. 240(1), 45–58 (2011)Picchini, U.: Inference for SDE models via Approximate Bayesian Computation. J. Comput.Graph. Stat. 23(4), 1080–1100 (2014)Picchini, U., Forman, J. L.: Accelerating inference for diffusions observed with measurement errorand large sample sizes using Approximate Bayesian Computation. J. Stat. Comput. Simul. 86(1), 195–213 (2016)Picchini, U., Samson, A.: Coupling stochastic EM and approximate Bayesian computation forparameter inference in state-space models. Comput. Stat. 33(1), 179–212 (2018)Pons, O.: Functional Estimation for Density, Regression Models and Processes. World ScientificPublishing, Singapore (2011)Prangle, D.: Adapting the abc distance function. Bayesian Anal. 12(1), 289–309 (2017)Prangle, D.: Summary Statistics in Approximate Bayesian Computation. In: Handbook of Ap-proximate Bayesian Computation: Chapman & Hall: pp 125–152 (2018)Prangle, D., Blum, M. G. B., Popovic, G., Sisson, S. A.: Diagnostic tools for approximate Bayesiancomputation using the coverage property. Aust. NZ. J. Stat. 56(4), 309–329 (2014)Quinn, B., Clarkson, I., Mckilliam, R.: On the periodogram estimators of periods from interleavedsparse, noisy timing data. In: 2014 IEEE Stat. Signal Processing Workshop: pp 232–235 (2014)R Development Core Team: R: A Language and Environment for Statistical Computing. R Foun-dation for Statistical Computing, Vienna, Austria (2011)Robert, C. P.: Approximate Bayesian Computation: A Survey on Recent Results. In: Cools,R., Nuyens, D. (eds) Monte Carlo and Quasi-Monte Carlo Methods: Springer InternationalPublishing, Cham: pp 185–205 (2016)van Rotterdam, A., Lopes da Silva, F., van den Ende, J., Viergever, M. A., Hermans, A.: A modelof the spatial-temporal characteristics of the alpha rhythm. Bull. Math. Biol. 44, 283–305 (1982)Sason, I., Verd´u, S.: f -Divergence Inequalities. IEEE T. Inform. Theory 62(11), 5973–6006 (2016)Sisson, S. A., Fan, Y., Beaumont, M.: Handbook of Approximate Bayesian Computation. Chap-man & Hall/CRC Handbooks of Modern Statistical Methods: CRC Press, Taylor & FrancisGroup (2018)Strømmen Melbø, A. H., Higham, D. J.: Numerical simulation of a linear stochastic oscillatorwith additive noise. Appl. Numer. Math. 51(1), 89–99 (2004)Sun, L., Lee, C., Hoeting, J. A.: Parameter inference and model selection in deterministic andstochastic dynamical models via approximate Bayesian computation: modeling a wildlife epi-demic. Environmetrics 26(7), 451–462 (2015)Tancredi, A.: Approximate Bayesian inference for discretely observed continuous-time multi-statemodels. Biometrics (2019) https://doi.org/10.1111/biom.13019https://doi.org/10.1111/biom.13019