close

Enter

Log in using OpenID

embedDownload
Sequential Monte Carlo
Informative Observations
with
Lawrence Murray (CSIRO) and Pierre Del Moral (UNSW)
CSIRO Mathematics, Informatics and Statistics
www.csiro.au
Highly
Problem
• Consider a continuous-time Markov process X(t) ∈ Rd.
• For a sequence of times t0, . . . , tn, write X0:n ≡ {X(t0), . . . , X(tn)}.
• Two scenarios of interest:
– Given x0 and xn, simulate bridges X1:n−1 ∼ p(X1:n−1 | x0, xn, θ)
and estimate the transition density (normalising constant)
p(xn | x0, θ).
– Given yn ∼ p(yn | xn, θ) and X0 ∼ p(X0 | θ), simulate X0:n ∼
p(X0:n | yn, θ) and estimate the marginal likelihood (normalising
constant) p(yn | θ).
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 2 of 30
Motivation
• Observations are highly informative, e.g.
– observation noise is narrow, or
– observations are sparse in time relative to the timescale of
the state process, or
– current θ is poor, e.g. early steps of the Markov chain in
PMCMC.
• Transition density p(xk | xk−1, θ) can be simulated but not evaluated
pointwise.
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 3 of 30
Basic idea
• Introduce a sequence of weighting and (optional) resampling
steps at intermediate times.
• The weighting function is a design choice, will demonstrate a
Gaussian process approach.
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 4 of 30
Importance sampling approach
Firstly note:
p(xn | xn−1)p(X1:n−1 | x0)
p(xn | x0)
∝ p(xn | xn−1)p(X1:n−1 | x0).
p(X1:n−1 | x0, xn) =
Simulate from p(X1:n−1 | x0), weight with p(xn | xn−1) [Pedersen,
1995].
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 5 of 30
Sequential Monte Carlo
Further observe:
p(X1:n−1 | x0) =
p(xn | xn−1)
=
p(xn | x0)
so that
n−1
Y
p(Xk | xk−1)
k=1
n−1
Y
p(xn | xk )
k=1 p(xn | xk−1 )
p(xn | xn−1)p(X1:n−1 | x0)
p(xn | x0)



n−1
n−1
p(x
|
x
)
Y
Y
n
k 

= 
p(Xk | xk−1)
k=1 p(xn | xk−1 ) k=1


n−1
p(x
|
x
)p(X
|
x
)
Y
n
k
k
k−1 

= 
p(x
|
x
)
k=1
n
k−1
Which suggests an SMC implementation [Lin et al., 2010].
p(X1:n−1 | x0, xn) =
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 6 of 30
Algorithm 1
For any bounded function f on (Rd)n−1, we have

E [f (X1:n−1) | x0, xn] = E f (X1:n−1)
n−1
Y
k=1
Gk (X
k−1:k 
) x0
with the weight functions
Gk (Xk−1:k ) :=
p(xn | xk )
.
p(xn | xk−1)
So for k = 1, . . . , n − 1 we can incrementally simulate p(Xk | xk−1),
weight with Gk (Xk−1:k ), and optionally resample.
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 7 of 30
Ornstein–Uhlenbeck process
Consider the Ornstein–Uhlenbeck process satisfing the Ito SDE:
dx = (θ1 − θ2x) dt + θ3 dW,
with parameters θ1 = 0.0187, θ2 = 0.2610 and θ3 = 0.0224 [Aït-Sahalia,
1999]. For step size ∆t, the transition density is
2
p (x(t + ∆t) | x(t)) = φ µ(∆t), σ (∆t) ,
with
θ1
θ1 
+ x(t) −  exp (−θ2∆t)
µ(∆t) =
θ2
θ2
2
θ
σ 2(∆t) = 3 (1 − exp (−2θ2∆t)) .
2θ2


Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 8 of 30
Ornstein–Uhlenbeck process
Bootstrap
Bridge
0.2
0.15
0.15
0.1
0.1
0.05
0.05
0
0
x
0.2
-0.05
-0.05
0
0.2
0.4
0.6
t
0.8
1
0
0.2
0.4
0.6
t
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 9 of 30
0.8
1
Ornstein–Uhlenbeck process
MSE(log z)-1 Mean(t)-1
ESS(z) Mean(t
-5
-2
10
10
10-6
Bridge
10-3
10-7
10-4
10-8
-9
10
-5
-9
10
-8
10
-7
10
Bootstrap
-6
10
-5
10
10
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 10 of 30
-5
10
-4
10
1
Bootstrap
Ornstein–Uhlenbeck process
MSE(log z)-1 Mean(t)-1
ESS(z) Mean(t)-1
CAR(z) Mean
-2
-5
10
10
10-6
10-3
10-7
10-4
10-8
-5
-8
10
-7
10
Bootstrap
-6
10
-5
10
10
-9
-5
10
-4
-3
10
10
-2
10
Bootstrap
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 11 of 30
10
-9
10
-8
10
-7
10
Bootstrap
Ornstein–Uhlenbeck process
ESS(z) Mean(t)
-1
CAR(z) Mean(t)
-2
-1
-5
0
10
-6
10
-3
0
10-7
0-4
0-5 -5
10
10-8
-4
-3
10
10
10-9 -9
10
-2
10
10-8
Bootstrap
10-7
Bootstrap
10-6
10-5
Z
1 X
CAR(z) =
2 ci − 1 ,
Z i=1
where ci is the sum of the ith smallest elements of z.


Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 12 of 30
Federal Funds Rate
x
Bootstrap
Bridge
0.12
0.12
0.115
0.115
0.11
0.11
0.105
0.105
0.1
0.1
0.095
0.095
0.09
0.09
0.085
0.085
0.08
0.08
0
0.5
1
1.5
2
t
2.5
3
3.5
4
0
0.5
1
1.5
2
t
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 13 of 30
2.5
3
3.5
4
Federal Funds Rate
MSE(log z)-1 Mean(t)-1
ESS(z) Mean
-6
-2
10
10
10-8
10-3
Bridge
-10
10
10-4
10-12
-5
10
-14
10
10-16 -16
10
-14
10
-12
-10
10
10
Bootstrap
-8
10
-6
10
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 14 of 30
10-6 -6
10
10-5
10-4
Bootstrap
Federal Funds Rate
MSE(log z)-1 Mean(t)-1
4
ESS(z) Mean(t)-1
CAR(z) Mean
-2
-6
10
10
10-3
10-7
10-4
10-8
-5
-9
10
-12
-10
10
10
Bootstrap
-8
10
-6
10
10-6 -6
10
10
-5
10
-4
10
Bootstrap
-3
10
-2
10
10-10 -10
10
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 15 of 30
10-9
10-8
Bootstrap
Federal Funds Rate
ESS(z) Mean(t)-1
CAR(z) Mean(t)-1
-6
10
10-7
10-8
-9
10
-5
10
-4
10
Bootstrap
-3
10
-2
10
10-10 -10
10
10-9
10-8
Bootstrap
10-7
10-6
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 16 of 30
Algorithm 2
For any bounded function f on (Rd)n−1, we have
n | xn−1 ) Qn−1
x
f (x1:n−1) p(x
H
(X
)
0
k
k−1:k
q(xn | xn−1 ) k=1
n | xn−1 ) Qn−1
E p(x
H
(X
)
k−1:k x0
q(xn | xn−1 ) k=1 k
E [f (X1:n−1) | x0, xn] =
E
with the weight functions
Hk (Xk−1:k ) :=
q(xn | xk )
q(xn | xk−1)
and chosen positive functions q(xn | xk ) for k = 1, . . . , n − 1, with
q(xn | x0) = 1.
So for k = 1, . . . , n − 1 we can incrementally simulate p(Xk | xk−1),
weight with Hk (Xk−1:k ), and optionally resample, then finally weight
with p(xn | xn−1)/q(xn | xn−1).
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 17 of 30
.
SIR
Stochastic SIR (susceptible/infectious/recovered) model.
State satisfies the Ito SDEs
dS = −βSI dt + σ1S dW1
dI = (βSI − νI) dt − σ1S dW1 + σ2I dW2
dR = νI dt − σ2I dW2.
with parameters β, ν, σ1 and σ2.
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 18 of 30
SIR
Convert to Stratonovich SDEs [Kloeden and Platen, 1992, p. 157]...


1
2
dS = −βSI − σ1 S  dt + σ1S ◦ dW1
2


1
1
2
2
dI = βSI − νI + σ1 S − σ2 I  dt − σ1S ◦ dW1 + σ2I ◦ dW2
2
2


1
dR = νI + σ22I  dt − σ2I ◦ dW2.
2
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 19 of 30
SIR
...and then to ODEs with a discrete-time noise innovation [Wilkie,
2004]:
dS
1
∆W1
= −βSI − σ12S + σ1S
dt
2
∆t
dI
1
1
∆W1
∆W2
= βSI − νI + σ12S − σ22I − σ1S
+ σ2 I
dt
2
2
∆t
∆t
dR
1
∆W2
= νI + σ22I − σ2I
,
dt
2
∆t
where each noise term ∆Wn ∼ N (0, ∆t) is an increment of the
Wiener process over a time step of size ∆t.
This permits the use of a numerical integrator that produces an
estimate of the discretisation error.
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 20 of 30
SIR
• An embedded Runge–Kutta scheme gives both order p and order
p−1 solutions via the same stage computations, with the difference
between the two an estimate of the error. The RK4(3)5[2R+]C
method [Carpenter and Kennedy, 1994] is used here.
• An adaptive scheme uses this estimate of the error to adjust
step sizes so as to enforce some error tolerance. Here, the error
tolerance for each state variable X is prescribed as
δabs + δrel|X(t)|.
with δabs = 10−2 (the absolute tolerance) and δrel = 10−5 (the
relative tolerance).
• It then makes sense to introduce a comparable -ball for observations:
Y (t) ∼ U (X(t) − δabs − δrel|X(t)|, X(t) + δabs + δrel|X(t)|) .
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 21 of 30
Data
• Observations of Russian influenza at a boarding school [Anonymous,
1978]
• 763 students
• I observed daily for 14 days
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 22 of 30
Weight functions
Fit a Gaussian process to the time series of each observed variable,
then construct weight functions based on these:
Y (t) ∼ GP (µ(t), C(∆t)) ,
with mean function µ(t) and covariance function C(∆t).
For this example, take a mean function µ(t) = 0 and squared exponential
covariance function parameterised by α and β:
1
C(∆t) = α exp − (∆t)2 .
2β


Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 23 of 30
Weight functions
The functions q(xn | xk ) are then constructed by conditioning the
Gaussian process on the current state:
q(xn | xk ) := φ(µk , σk2 )
with
C(tn − tk )
xk
α
C(tn − tk )2
2
.
σk = α −
α
Note that q is obtained by conditioning on the current state only,
and not past states. This is a computational concession so that the
algorithmic complexity remains linear in the number of particles N .
µk =
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 24 of 30
SIR
0.005
1
0.0045
0.8
0.004
0.6
0.003
ν
β
0.0035
0.4
0.0025
0.002
0.2
0.0015
0.001
0
0
20000
40000
60000
80000
100000
0.04
0
20000
40000
60000
80000
100000
0
20000
40000
60000
80000
100000
1
0.035
0.8
0.03
0.6
σ2
σ1
0.025
0.02
0.4
0.015
0.01
0.2
0.005
0
0
0
20000
40000
60000
80000
100000
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 25 of 30
SIR
700
600
S
500
400
300
200
100
0
2
4
6
8
10
12
0
2
4
6
8
10
12
400
350
300
I
250
200
150
100
50
700
Lawrence 600
Murray Sequential Monte Carlo with Highly Informative Observations: Slide 26 of 30
500
50
SIR
0
2
4
6
0
2
4
6
8
10
12
8
10
12
700
600
R
500
400
300
200
100
0
t
1400
4
160
3.5
140
1000
3
120
800
2.5
100
2
80
1.5
60
1
40
0.5
20
1200
600
400
200
0
0.0018
0.002
0.0022
0.0024
0.0026
0.0028
0.003
0.0032
0.0034
β
0
5
4
3
2
1
0
0.3 0.4 0.5 0.6 0.7 0.8
ν
0
0.0050.010.0150.020.025
σ1
0.20.250.30.350.40.450.50.55
σ2
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 27 of 30
Summary
• Works in a multivariate setting.
• Supports partial observation of the state vector xn.
• Straightforward to extend with an observation model, and should
be particularly useful when observations are informative.
• Supports higher-order discretisations of the forward state process
than Euler–Maruyama.
• Requires only that the forward state process can be simulated and
not that its transition density p(xk | xk−1, θ) can be evaluated
(with some caveats).
• Does not require simulation of the backward state process.
• Will be available, along with some examples, in the next version
of LibBi (out soon, www.libbi.org).
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 28 of 30
References
Anonymous. Influenza in a boarding school. British Medical Journal, 1:587, March 1978. URL http:
//www.ncbi.nlm.nih.gov/pmc/articles/PMC1603269/?page=2.
Y. Aït-Sahalia. Transition densities for interest rate and other nonlinear diffusions. The Journal of Finance,
54(4):1361–1395, 1999. ISSN 1540-6261. doi: 10.1111/0022-1082.00149.
M. H. Carpenter and C. A. Kennedy. Fourth-order 2N-storage Runge-Kutta schemes. Technical Report
Technical Memorandum 109112, National Aeronautics and Space Administration, June 1994.
P. E. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations. Springer, 1992.
M. Lin, R. Chen, and P. Mykland. On generating Monte Carlo samples of continuous diffusion bridges. Journal
of the American Statistical Association, 105(490):820–838, 2010. doi: 10.1198/jasa.2010.tm09057.
A. Pedersen. A new approach to maximum likelihood estimation for stochastic differential equations based on
discrete observations. Scandinavian Journal of Statistics, 22(1):55–71, 1995.
J. Wilkie. Numerical methods for stochastic differential equations. Physical Review E, 70, 2004.
Lawrence Murray Sequential Monte Carlo with Highly Informative Observations: Slide 29 of 30
CSIRO
Lawrence Murray
t +61 8 9333 6480
e lawrence.murray@csiro.au
w http://www.csiro.au
CSIRO Mathematics, Informatics and Statistics
www.csiro.au
1/--pages
Report inappropriate content