De Gunst, M.

Vujacic, I.

DSpace at VU

Simultaneous perturbation gradient
approximation based Metropolis adjusted
Langevin Markov chain Monte Carlo for inference
of ordinary differential equations
Ivan Vujacˇic´∗1 and Mathisca de Gunst1
1VU University Amsterdam
Abstract: The problem of parameter estimation for models defined by a system
of ordinary differential equations (ODEs) is considered. The most efficient way to
explore the parameter space is by using derivative information. Usual approaches
for obtaining the gradient in ODEs setting like solving sensitivity equations and
using finite difference formulas are computationally costly and not scalable to large
scale systems. In this paper we use simultaneous perturbation gradient approxima-
tion (SPGA), originally proposed in stochastic optimization literature, as a substi-
tute for the gradient in Metropolis adjusted Langevin algorithm (MALA). The ob-
tained algorithm, called Simultaneous Perturbation Gradient Approximation based
Metropolis Adjusted Langevin Markov chain Monte Carlo (SPGA MALA), requires
at most three integration of the ODE system per MCMC step, regardless of the
dimension of the system. This fixed computational costs makes SPGA MALA ap-
plicable to large scale systems. On the other hand, its efficiency is comparable to
that of MALA. We demonstrate its performance of via simulations.
Keywords: Metropolis adjusted Langevin Markov Chain Monte Carlo methods,
Simultaneous perturbation gradient approximation, ordinary differential equations,
parameter estimation.
AMS subject classifications: 60J22, 65C40, 62F15.
1 Introduction
Systems of ordinary differential equations (ODEs) are widely used in science and
engineering for the mathematical modelling of various dynamic processes. We
consider the system of the form{
x′(t) = f(x(t), t;θ), t ∈ [0, T ],
x(0) = ξ,
(1)
where x(t) = (x1(t), . . . , xd(t))
> ∈ Rd is a state vector, ξ in Ξ ⊂ Rd is the initial
condition, θ in Θ ⊂ Rp is a parameter and f is a known function. Given the
values of ξ and θ, we denote the solution of (1) by x(t;θ, ξ). Let us assume that
a process is modelled by the system (1) with ξ0 known and θ0 unknown. For
∗Corresponding author: i.vujacic@vu.nl
2 Vujacˇic´, De Gunst– SPGA MALA
simplicity, assume that we have noisy observations yi(tj), j = 1, . . . , n of all the
states xi(t;θ0, ξ0), i = 1, . . . , d at time points tj ∈ [0, T ], j = 1, . . . , n:
yi(tj) = xi(tj ;θ0, ξ0) + εi(tj), i = 1, . . . , d; j = 1, . . . , n,
where εi(tj) ∼ N (0, σ2i ). The problem is to estimate θ0 from the data Y =
(yi(tj))ij . The methodology presented here can also be used if ξ0 is unknown and
some of the states are unobserved.
In this paper, we adopt Bayesian approach to inference. For some prior density
pi of θ the posterior density is
p(θ|Y, ξ0,σ) ∝ pi(θ)
d∏
j=1
N{Yj,·|X(θ, ξ0)j,·, σjIn},
where σ = (σ1, . . . , σd), X(θ, ξ0) = (xi(tj ;θ, ξ0))ij and In is an identity matrix of
order n. For exploring the parameter space there is an advantage in using gradient
information in MCMC and optimization methods [3, 4]; for concrete examples in
ODE estimation setting see [4, 6]. In the problem we consider, the gradient of the
log-likelihood can be obtained by solving sensitivity equations, which are of order
dp or via the finite difference formulas, which require solving the ODE system at
least p times. Both approaches are computationally costly and not scalable to large
scale systems.
In this paper, we avoid huge computational burden by using simultaneous per-
turbation stochastic approximation (SPGA), introduced by Spall [7]. To obtain
SPGA, the system of the form (1) need be solved at most 2 times, regardless of
the dimension of the system. By using SPGA instead of the gradient in Metropo-
lis adjusted Langevin Markov Chain Monte Carlo (MALA) we obtain a method,
which we call SPGA MALA, that can be used for large scale systems. Although
there is some loss in efficiency of SPGA MALA due to using an approximation of
the derivative this is outweighed by huge computational savings achieved.
The rest of the paper is organized as follows. In sections 2 and 3 reviews of
MALA and SPGA are provided, respectively. Section 4 introduces the proposed
method. In Section 5 we compare performance of MALA and SPGA MALA on
simulated data for various models.
2 Metropolis adjusted Langevin Markov chain Monte
Carlo (MALA)
For the probability density p(θ) let L(θ) = log{p(θ)} denote the log-density. The
MALA proposal [4, p.130] is
θ∗ = θk + 2M∇θL(θk)/2 + 
√
Mzk, (2)
where θk is the value at k-th step, z ∼ N (z|0, Ip),  > 0 is the step size and M is
the weight matrix. The proposal density and acceptance probability are
q(θ∗|θk) = N (θ∗|µ(θk, ), 2M),
α = min{1, p(θ∗)q(θk|θ∗)/p(θk)q(θ∗|θk)}, (3)
Proceedings of the 19th EYSM in Prague 2015 3
respectively, where µ(θk, ) = θk + 2M∇θL(θk)/2. The advantage of MALA over
random walk Metropolis algorithm is that it uses the gradient information which
leads to better exploration of the parameter space. Disadvantage is that it requires
selection of the weight matrix M. In [4], a fully automated algorithm is proposed
for this but it cannot be used in our setting because it requires derivatives. For
more details regarding MALA see [2, 4].
3 Simultaneous perturbation gradient approxima-
tion (SPGA)
In order to estimate partial derivatives via finite difference (FD) approximation the
parameter perturbations are performed along each coordinate separately. For ex-
ample, the estimate of the j-th partial derivative of L(θk) via the central difference
formula is
∂L(θk)
∂θj
≈ L(θ
k + hej)− L(θk − hej)
2h
,
where ej is the j-th unit vector and h is sufficiently small. This requires 2p eval-
uations of L. With simultaneous perturbation (SP), introduced by Spall [7], all
elements of θk are randomly perturbed together. The two sided simultaneous per-
turbation gradient approximation (SPGA) is
∇̂θL(θk) = L(θ
k + h∆)− L(θk − h∆)
2h
(∆−11 ,∆
−1
2 , . . . ,∆
−1
p )
>, (4)
where ∆ = (∆1,∆2, . . . ,∆p)
> is usually a random vector of independent Bernoulli
random variables that take values −1 and 1 with probability 0.5, although other
choices are possible. Two sided SPGA requires two evaluations of L regardless
of the dimension p. FD approximation is superior to SP approximation as an
estimator of the gradient. However, Spall [7] showed that when used in stochas-
tic optimization setting they achieve the same level of statistical accuracy for a
given number of iterations in terms of estimation of the optimum of the objective
function. In the next section, we follow the same idea but in the MCMC setting.
4 Simultaneous perturbation gradient approxima-
tion based Metropolis adjusted Langevin Markov
chain Monte Carlo (SPGA MALA)
SPGA MALA proposal is obtained by substituting the gradient in MALA proposal
(2) with its SPGA, defined in (4):
θ∗ = θk + 2M∇̂θL(θk)/2 + 
√
Mzk. (5)
In view of the MALA proposal density in (3), we require that the density of θ∗
given θk and ∆ is q(θ∗|θk,∆) = N (θ∗|µ̂(θk, ,∆), 2M), where µ̂(θk, ,∆) = θk +
2M∇̂θL(θk)/2. Since ∆ can take 2p values with equal probability it follows that
4 Vujacˇic´, De Gunst– SPGA MALA
the density of θ∗ given θk is the mixture density q(θ∗|θk) = 12p
∑
∆ q(θ
∗|θk,∆).
The proposal mechanism (5) with proposal density q and standard accepetance
probability α = min{1, p(θ∗)q(θk|θ∗)/p(θk)q(θ∗|θk)} defines a valid Markov chain;
it is simply Metropolis Hastings (MH) algorithm where the proposal is a mixture
density q. However, evaluating α is intractable for large p. Because of this, instead
of α we use
α∆ = min{1, p(θ∗)q(θk|θ∗,∆)/p(θk)q(θ∗|θk,∆)}.
In other words, instead of using q which involves calculation of each q(θ∗|θk,∆)
for 2p possible values of ∆, we use the acceptance ratio which involves calculation
of q(θ∗|θk,∆) only for the drawn value of ∆.
This algorithm defines a valid Markov chain since it can be viewed as Metropolis-
Hastings-Green(MHG) algorithm [2, p.41]. MHG algorithm allows state-dependent
mixing or random proposals [1], meaning that on each step the proposal distribution
need not be fixed but can belong to a countable family of proposal distributions.
In our case it is a finite family {q(θ∗|θk,∆) : ∆ = (∆1, . . . ,∆p),∆i ∈ {−1, 1}}.
Using the random proposal instead of the mixture proposal comes with a price. As
pointed out in the discussion section of the article [1], the random proposal method
is less efficient because it accepts fewer proposals. This reduces efficiency of SPGA
MALA. The second reason for reduced efficiency of SPGA MALA is that instead
of the gradient its approximation is used. However, it is clear that SPGA MALA
will be much faster than MALA for large scale systems.
5 Numerical results
In this section we compare the desribed algorithm to MALA on simulated data
generated from the following models.
Fitz Hugh Nagumo (FHN) example. Fitz-Hugh Nagumo system [4] models
the behaviour of spike potentials in the giant axon of squid neurons. It has the
form
x′1(t) = θ3{x1(t)− x1(t)3/3 + x2(t)},
x′2(t) = − 1θ3 {x1(t)− θ1 + θ2x2(t)}.
We have used different notation than in [4], namely (x1, x2) for (V,R) and (θ1, θ2, θ3)
for (a, b, c). We set θ = (0.2, 0.2, 3) and ξ = (−1, 1).
α- pinene example. The following model describes the thermal isomerization of
α-pinene [8].
x′1(t) = −(θ1 + θ2)x1(t),
x′2(t) = θ1x1(t),
x′3(t) = θ2x1(t)− (θ3 + θ4)x3(t) + θ5x5(t),
x′4(t) = θ3x3(t),
x′5(t) = θ4x3(t)− θ5x5(t).
The values of the parameters that we used are θ = (0.1, 0.1, 0.3, 0.1, 0.3) and ξ =
(1, 0, 0, 0, 0).
Hockin model. In [5], a model of the extrinsic blood coagulation is developed
and consists of 34 differential equations and 42 rate constants. Due to lack of space
Proceedings of the 19th EYSM in Prague 2015 5
we do not present the model but refer the reader to the aforementioned reference.
We fixed 32 parameters and estimated the remaining 10. The value of the selected
parameter was set to θ = (0.1, 0.4, 0.1, 0.32, 0.2, 1.05, 2.4, 6, 1.8, 8.2).
From each of the models presented above we generated 200 data points on the
interval [0, 20] and added Gaussian-distributed noise with standard deviation equal
to 0.5. In SPGA (see (4)) we set h = 10e−5 while the gradient in MALA is obtained
by solving sensitivity equations. Ideally, the tuning parameters in MALA should
be chosen in such a way that acceptance rate is between 40% and 70%. As it was
pointed out in Section 2, tuning of MALA is an issue. To simplify, for both MALA
and SPGA MALA we set M = Ip,  = 0.0002p
−1/3 in all the simulations. This
selection achieves the desired acceptance rate in FHN model and was used in [4].
For the other two models this is not the case. However, the most important thing
here is to compare the performance of these two methods for the same selection of
tuning parameters. For comparing sampling efficiency we followed approach used
in [4]. A single Markov chain was initialized on the true mode and 5000 posterior
samples were collected. The effective sample size (ESS) for each parameter was
calculated; the minimum of ESS was used to calculate the time per effectively
independent sample. For each method we ran 10 simulations, using the same data
set. The methods were implemented in the interpreted language MATLAB and all
computations were carried out on an Intel Core i5 computer with 1.3 GHz processor
speed and 4 GB of memory. The results of our simulations are presented in Table
5.
The results of FHN example demonstrate loss in efficiency of SPGA MALA with
respect to MALA; see Section 4 for the discussion. In α-pinene example MALA is
still faster even though the sensitivity equations are of order 25. This is because
the original system and the system of sensitivity equations are both linear. The ex-
ample of Hockin model show the advantage of SPGA MALA. Sensitivity equations
are of order 340 and this heavily affects the computatation time of MALA. On the
other hand, the computation time of SPGA MALA is much smaller compared to
that of MALA, making it much better in terms of the relative speed per effectively
independent sample.
Acknowledgements: This research is supported by the Dutch Technology Foun-
dation STW, which is part of the Netherlands Organisation for Scientific Research
(NWO), and which is partly funded by Ministry of Economic Affairs. Mark Giro-
lami and Ben Calderhead are acknowledged for making their MATLAB code used in
[4] freely available. We thank Itai Dattner and Bart Bakker for useful comments.
References
[1] J. Besag, P. Green, D. Higdon, and K. Mengersen. Bayesian computation and
stochastic systems. Statistical science, pages 3–41, 1995.
[2] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng. Handbook of Markov Chain
Monte Carlo. CRC press, 2011.
6 Vujacˇic´, De Gunst– SPGA MALA
Model
Sampling
method
Time
(s)
Mean ESS
(θ)
Total time
/minimum
mean ESS
Relative
speed
d = 2
p = 3
(θ1, θ2, θ3)
FHN
MALA 363.6 145, 30, 109 12.12 3.4
SPGA
MALA
623.2 84 , 15, 48 41.55 1
d = 5
p = 5
(θ1, θ2, θ3, θ4, θ5)
α-Pinene
MALA 63.7 59, 58, 17, 11, 6 10.62 2.54
SPGA
MALA
134.8 187, 96, 6, 23 , 5 26.96 1
d = 34
p = 10
(θ1, θ2, θ3, θ4, θ5,
θ6, θ7, θ8, θ9, θ10)
Hockin
MALA 1.03e+04
5 6 8 7 7
7 6 5 6 8
2060 1
SPGA
MALA
180.5
7 6 8 8 7
6 6 5 4 7
45.13 45.65
Table 1: Summary of results for 10 runs of the model parameter sampling schemes
for different models with 5000 posterior samples.
[3] A. R. Conn, K. Scheinberg, and L. N. Vicente. Introduction to derivative-free
optimization, volume 8. Siam, 2009.
[4] M. Girolami and B. Calderhead. Riemann manifold langevin and hamiltonian
monte carlo methods. Journal of the Royal Statistical Society: Series B (Sta-
tistical Methodology), 73(2):123–214, 2011.
[5] M. F. Hockin, K. C. Jones, S. J. Everse, and K. G. Mann. A model for the
stoichiometric regulation of blood coagulation. Journal of Biological Chemistry,
277(21):18322–18333, 2002.
[6] A. Raue, M. Schilling, J. Bachmann, A. Matteson, M. Schelke, D. Kaschek,
S. Hug, C. Kreutz, B. D. Harms, F. J. Theis, et al. Lessons learned from
quantitative dynamical modeling in systems biology. PloS one, 8(9):e74335,
2013.
[7] J. C. Spall. Multivariate stochastic approximation using a simultaneous per-
turbation gradient approximation. Automatic Control, IEEE Transactions on,
37(3):332–341, 1992.
[8] I. B. Tjoa and L. T. Biegler. Simultaneous solution and optimization strategies
for parameter estimation of differential-algebraic equation systems. Industrial
& Engineering Chemistry Research, 30(2):376–385, 1991.


Simultaneous perturbation gradient approximation based Metropolis adjusted Langevin Markov chain Monte Carlo for inference of ordinary differential equations

de Gunst, M.C.M.

VU Research Portal

VU Research PortalSimultaneous perturbation gradient approximation based Metropolis adjustedLangevin Markov chain Monte Carlo for inference of ordinary differential equationsVujacic, I.; de Gunst, M.C.M.published inProceedings of the 19th European Young Statisticians Meeting2015document versionPublisher's PDF, also known as Version of recordLink to publication in VU Research Portalcitation for published version (APA)Vujacic, I., & de Gunst, M. C. M. (2015). Simultaneous perturbation gradient approximation based Metropolisadjusted Langevin Markov chain Monte Carlo for inference of ordinary differential equations. In S. StanislavNagy (Ed.), Proceedings of the 19th European Young Statisticians Meeting (pp. 160-165). Publishing House ofthe Faculty of Mathematics and Physics, Charles University in Prague.General rightsCopyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright ownersand it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.            • Users may download and print one copy of any publication from the public portal for the purpose of private study or research.            • You may not further distribute the material or use it for any profit-making activity or commercial gain            • You may freely distribute the URL identifying the publication in the public portal ?Take down policyIf you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediatelyand investigate your claim.E-mail address:vuresearchportal.ub@vu.nlDownload date: 27. May. 2021Simultaneous perturbation gradientapproximation based Metropolis adjustedLangevin Markov chain Monte Carlo for inferenceof ordinary differential equationsIvan Vujačić∗1 and Mathisca de Gunst11VU University AmsterdamAbstract: The problem of parameter estimation for models defined by a systemof ordinary differential equations (ODEs) is considered. The most efficient way toexplore the parameter space is by using derivative information. Usual approachesfor obtaining the gradient in ODEs setting like solving sensitivity equations andusing finite difference formulas are computationally costly and not scalable to largescale systems. In this paper we use simultaneous perturbation gradient approxima-tion (SPGA), originally proposed in stochastic optimization literature, as a substi-tute for the gradient in Metropolis adjusted Langevin algorithm (MALA). The ob-tained algorithm, called Simultaneous Perturbation Gradient Approximation basedMetropolis Adjusted Langevin Markov chain Monte Carlo (SPGA MALA), requiresat most three integration of the ODE system per MCMC step, regardless of thedimension of the system. This fixed computational costs makes SPGA MALA ap-plicable to large scale systems. On the other hand, its efficiency is comparable tothat of MALA. We demonstrate its performance via simulations.Keywords: Metropolis adjusted Langevin Markov Chain Monte Carlo methods,Simultaneous perturbation gradient approximation, ordinary differential equations,parameter estimation.AMS subject classifications: 60J22, 65C40, 62F15.1 IntroductionSystems of ordinary differential equations (ODEs) are widely used in science andengineering for the mathematical modelling of various dynamic processes. Weconsider the system of the form{x′(t) = f(x(t), t;θ), t ∈ [0, T ],x(0) = ξ,(1)where x(t) = (x1(t), . . . , xd(t))> ∈ Rd is a state vector, ξ in Ξ ⊂ Rd is the initialcondition, θ in Θ ⊂ Rp is a parameter and f is a known function. Given thevalues of ξ and θ, we denote the solution of (1) by x(t;θ, ξ). Let us assume thata process is modelled by the system (1) with ξ0 known and θ0 unknown. For∗Corresponding author: i.vujacic@vu.nl2 Vujačić, De Gunst– SPGA MALAsimplicity, assume that we have noisy observations yi(tj), j = 1, . . . , n of all thestates xi(t;θ0, ξ0), i = 1, . . . , d at time points tj ∈ [0, T ], j = 1, . . . , n:yi(tj) = xi(tj ;θ0, ξ0) + εi(tj), i = 1, . . . , d; j = 1, . . . , n,where εi(tj) ∼ N (0, σ2i ). The problem is to estimate θ0 from the data Y =(yi(tj))ij . The methodology presented here can also be used if ξ0 is unknown andsome of the states are unobserved.In this paper, we adopt Bayesian approach to inference. For some prior densityπ of θ the posterior density isp(θ|Y, ξ0,σ) ∝ π(θ)d∏j=1N{Yj,·|X(θ, ξ0)j,·, σjIn},where σ = (σ1, . . . , σd), X(θ, ξ0) = (xi(tj ;θ, ξ0))ij and In is an identity matrix oforder n. For exploring the parameter space there is an advantage in using gradientinformation in MCMC and optimization methods [3, 4]; for concrete examples inODE estimation setting see [4, 6]. In the problem we consider, the gradient of thelog-likelihood can be obtained by solving sensitivity equations, which are of orderdp or via the finite difference formulas, which require solving the ODE system atleast p times. Both approaches are computationally costly and not scalable to largescale systems.In this paper, we avoid huge computational burden by using simultaneous per-turbation stochastic approximation (SPGA), introduced by Spall [7]. To obtainSPGA, the system of the form (1) need be solved at most 2 times, regardless ofthe dimension of the system. By using SPGA instead of the gradient in Metropo-lis adjusted Langevin Markov Chain Monte Carlo (MALA) we obtain a method,which we call SPGA MALA, that can be used for large scale systems. Althoughthere is some loss in efficiency of SPGA MALA due to using an approximation ofthe derivative this is outweighed by huge computational savings achieved.The rest of the paper is organized as follows. In sections 2 and 3 reviews ofMALA and SPGA are provided, respectively. Section 4 introduces the proposedmethod. In Section 5 we compare performance of MALA and SPGA MALA onsimulated data for various models.2 Metropolis adjusted Langevin Markov chain MonteCarlo (MALA)For the probability density p(θ) let L(θ) = log{p(θ)} denote the log-density. TheMALA proposal [4, p.130] isθ∗ = θk + ε2M∇θL(θk)/2 + ε√Mzk, (2)where θk is the value at k-th step, z ∼ N (z|0, Ip), ε > 0 is the step size and M isthe weight matrix. The proposal density and acceptance probability areq(θ∗|θk) = N (θ∗|µ(θk, ε), ε2M),α = min{1, p(θ∗)q(θk|θ∗)/p(θk)q(θ∗|θk)}, (3)Proceedings of the 19th EYSM in Prague 2015 3respectively, where µ(θk, ε) = θk + ε2M∇θL(θk)/2. The advantage of MALA overrandom walk Metropolis algorithm is that it uses the gradient information whichleads to better exploration of the parameter space. Disadvantage is that it requiresselection of the weight matrix M. In [4], a fully automated algorithm is proposedfor this but it cannot be used in our setting because it requires derivatives. Formore details regarding MALA see [2, 4].3 Simultaneous perturbation gradient approxima-tion (SPGA)In order to estimate partial derivatives via finite difference (FD) approximation theparameter perturbations are performed along each coordinate separately. For ex-ample, the estimate of the j-th partial derivative of L(θk) via the central differenceformula is∂L(θk)∂θj≈ L(θk + hej)− L(θk − hej)2h,where ej is the j-th unit vector and h is sufficiently small. This requires 2p eval-uations of L. With simultaneous perturbation (SP), introduced by Spall [7], allelements of θk are randomly perturbed together. The two sided simultaneous per-turbation gradient approximation (SPGA) is∇̂θL(θk) =L(θk + h∆)− L(θk − h∆)2h(∆−11 ,∆−12 , . . . ,∆−1p )>, (4)where ∆ = (∆1,∆2, . . . ,∆p)> is usually a random vector of independent Bernoullirandom variables that take values −1 and 1 with probability 0.5, although otherchoices are possible. Two sided SPGA requires two evaluations of L regardlessof the dimension p. FD approximation is superior to SP approximation as anestimator of the gradient. However, Spall [7] showed that when used in stochas-tic optimization setting they achieve the same level of statistical accuracy for agiven number of iterations in terms of estimation of the optimum of the objectivefunction. In the next section, we follow the same idea but in the MCMC setting.4 Simultaneous perturbation gradient approxima-tion based Metropolis adjusted Langevin Markovchain Monte Carlo (SPGA MALA)SPGA MALA proposal is obtained by substituting the gradient in MALA proposal(2) with its SPGA, defined in (4):θ∗ = θk + ε2M∇̂θL(θk)/2 + ε√Mzk. (5)In view of the MALA proposal density in (3), we require that the density of θ∗given θk and ∆ is q(θ∗|θk,∆) = N (θ∗|µ̂(θk, ε,∆), ε2M), where µ̂(θk, ε,∆) = θk +ε2M∇̂θL(θk)/2. Since ∆ can take 2p values with equal probability it follows that4 Vujačić, De Gunst– SPGA MALAthe density of θ∗ given θk is the mixture density q(θ∗|θk) = 12p∑∆ q(θ∗|θk,∆).The proposal mechanism (5) with proposal density q and standard accepetanceprobability α = min{1, p(θ∗)q(θk|θ∗)/p(θk)q(θ∗|θk)} defines a valid Markov chain;it is simply Metropolis Hastings (MH) algorithm where the proposal is a mixturedensity q. However, evaluating α is intractable for large p. Because of this, insteadof α we useα∆ = min{1, p(θ∗)q(θk|θ∗,∆)/p(θk)q(θ∗|θk,∆)}.In other words, instead of using q which involves calculation of each q(θ∗|θk,∆)for 2p possible values of ∆, we use the acceptance ratio which involves calculationof q(θ∗|θk,∆) only for the drawn value of ∆.This algorithm defines a valid Markov chain since it can be viewed as Metropolis-Hastings-Green(MHG) algorithm [2, p.41]. MHG algorithm allows state-dependentmixing or random proposals [1], meaning that on each step the proposal distributionneed not be fixed but can belong to a countable family of proposal distributions.In our case it is a finite family {q(θ∗|θk,∆) : ∆ = (∆1, . . . ,∆p),∆i ∈ {−1, 1}}.Using the random proposal instead of the mixture proposal comes with a price. Aspointed out in the discussion section of the article [1], the random proposal methodis less efficient because it accepts fewer proposals. This reduces efficiency of SPGAMALA. The second reason for reduced efficiency of SPGA MALA is that insteadof the gradient its approximation is used. However, it is clear that SPGA MALAwill be much faster than MALA for large scale systems.5 Numerical resultsIn this section we compare the desribed algorithm to MALA on simulated datagenerated from the following models.Fitz Hugh Nagumo (FHN) example. Fitz-Hugh Nagumo system [4] modelsthe behaviour of spike potentials in the giant axon of squid neurons. It has theformx′1(t) = θ3{x1(t)− x1(t)3/3 + x2(t)},x′2(t) = − 1θ3 {x1(t)− θ1 + θ2x2(t)}.We have used different notation than in [4], namely (x1, x2) for (V,R) and (θ1, θ2, θ3)for (a, b, c). We set θ = (0.2, 0.2, 3) and ξ = (−1, 1).α- pinene example. The following model describes the thermal isomerization ofα-pinene [8].x′1(t) = −(θ1 + θ2)x1(t),x′2(t) = θ1x1(t),x′3(t) = θ2x1(t)− (θ3 + θ4)x3(t) + θ5x5(t),x′4(t) = θ3x3(t),x′5(t) = θ4x3(t)− θ5x5(t).The values of the parameters that we used are θ = (0.1, 0.1, 0.3, 0.1, 0.3) and ξ =(1, 0, 0, 0, 0).Hockin model. In [5], a model of the extrinsic blood coagulation is developedand consists of 34 differential equations and 42 rate constants. Due to lack of spaceProceedings of the 19th EYSM in Prague 2015 5we do not present the model but refer the reader to the aforementioned reference.We fixed 32 parameters and estimated the remaining 10. The value of the selectedparameter was set to θ = (0.1, 0.4, 0.1, 0.32, 0.2, 1.05, 2.4, 6, 1.8, 8.2).From each of the models presented above we generated 200 data points on theinterval [0, 20] and added Gaussian-distributed noise with standard deviation equalto 0.5. In SPGA (see (4)) we set h = 10e−5 while the gradient in MALA is obtainedby solving sensitivity equations. Ideally, the tuning parameters in MALA shouldbe chosen in such a way that acceptance rate is between 40% and 70%. As it waspointed out in Section 2, tuning of MALA is an issue. To simplify, for both MALAand SPGA MALA we set M = Ip, ε = 0.0002p−1/3 in all the simulations. Thisselection achieves the desired acceptance rate in FHN model and was used in [4].For the other two models this is not the case. However, the most important thinghere is to compare the performance of these two methods for the same selection oftuning parameters. For comparing sampling efficiency we followed approach usedin [4]. A single Markov chain was initialized on the true mode and 5000 posteriorsamples were collected. The effective sample size (ESS) for each parameter wascalculated; the minimum of ESS was used to calculate the time per effectivelyindependent sample. For each method we ran 10 simulations, using the same dataset. The methods were implemented in the interpreted language MATLAB and allcomputations were carried out on an Intel Core i5 computer with 1.3 GHz processorspeed and 4 GB of memory. The results of our simulations are presented in Table5.The results of FHN example demonstrate loss in efficiency of SPGA MALA withrespect to MALA; see Section 4 for the discussion. In α-pinene example MALA isstill faster even though the sensitivity equations are of order 25. This is becausethe original system and the system of sensitivity equations are both linear. The ex-ample of Hockin model show the advantage of SPGA MALA. Sensitivity equationsare of order 340 and this heavily affects the computatation time of MALA. On theother hand, the computation time of SPGA MALA is much smaller compared tothat of MALA, making it much better in terms of the relative speed per effectivelyindependent sample.Acknowledgements: This research is supported by the Dutch Technology Foun-dation STW, which is part of the Netherlands Organisation for Scientific Research(NWO), and which is partly funded by Ministry of Economic Affairs. Mark Giro-lami and Ben Calderhead are acknowledged for making their MATLAB code used in[4] freely available. We thank Itai Dattner and Bart Bakker for useful comments.References[1] J. Besag, P. Green, D. Higdon, and K. Mengersen. Bayesian computation andstochastic systems. Statistical science, pages 3–41, 1995.[2] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng. Handbook of Markov ChainMonte Carlo. CRC press, 2011.6 Vujačić, De Gunst– SPGA MALAModelSamplingmethodTime(s)Mean ESS(θ)Total time/minimummean ESSRelativespeedd = 2p = 3(θ1, θ2, θ3)FHNMALA 363.6 145, 30, 109 12.12 3.4SPGAMALA623.2 84 , 15, 48 41.55 1d = 5p = 5(θ1, θ2, θ3, θ4, θ5)α-PineneMALA 63.7 59, 58, 17, 11, 6 10.62 2.54SPGAMALA134.8 187, 96, 6, 23 , 5 26.96 1d = 34p = 10(θ1, θ2, θ3, θ4, θ5,θ6, θ7, θ8, θ9, θ10)HockinMALA 1.03e+045 6 8 7 77 6 5 6 82060 1SPGAMALA180.57 6 8 8 76 6 5 4 745.13 45.65Table 1: Summary of results for 10 runs of the model parameter sampling schemesfor different models with 5000 posterior samples.[3] A. R. Conn, K. Scheinberg, and L. N. Vicente. Introduction to derivative-freeoptimization, volume 8. Siam, 2009.[4] M. Girolami and B. Calderhead. Riemann manifold langevin and hamiltonianmonte carlo methods. Journal of the Royal Statistical Society: Series B (Sta-tistical Methodology), 73(2):123–214, 2011.[5] M. F. Hockin, K. C. Jones, S. J. Everse, and K. G. Mann. A model for thestoichiometric regulation of blood coagulation. Journal of Biological Chemistry,277(21):18322–18333, 2002.[6] A. Raue, M. Schilling, J. Bachmann, A. Matteson, M. Schelke, D. Kaschek,S. Hug, C. Kreutz, B. D. Harms, F. J. Theis, et al. Lessons learned fromquantitative dynamical modeling in systems biology. PloS one, 8(9):e74335,2013.[7] J. C. Spall. Multivariate stochastic approximation using a simultaneous per-turbation gradient approximation. Automatic Control, IEEE Transactions on,37(3):332–341, 1992.[8] I. B. Tjoa and L. T. Biegler. Simultaneous solution and optimization strategiesfor parameter estimation of differential-algebraic equation systems. Industrial& Engineering Chemistry Research, 30(2):376–385, 1991.

English

https://research.vu.nl/files/1395170/IvanVujacic paper eysm19.pdf

Simultaneous perturbation gradient approximation based Metropolis adjusted Langevin Markov chain Monte Carlo for inference of ordinary differential equations

Abstract

Similar works

Full text

Available Versions

DSpace at VU

VU Research Portal