Revised Implicit Equal-Weights Particle Filter Jacob Skauvold 1 Joint - - PowerPoint PPT Presentation

revised implicit equal weights particle filter
SMART_READER_LITE
LIVE PREVIEW

Revised Implicit Equal-Weights Particle Filter Jacob Skauvold 1 Joint - - PowerPoint PPT Presentation

Revised Implicit Equal-Weights Particle Filter Jacob Skauvold 1 Joint work with Jo Eidsvik 1 , Peter Jan van Leeuwen 2 and Javier Amezcua 2 1 Dept. of Mathematical Sciences, NTNU, Norway 2 Data Assimilation Research Centre, University of Reading, UK


slide-1
SLIDE 1

Revised Implicit Equal-Weights Particle Filter

Jacob Skauvold1

Joint work with Jo Eidsvik1, Peter Jan van Leeuwen2 and Javier Amezcua2

  • 1Dept. of Mathematical Sciences, NTNU, Norway

2Data Assimilation Research Centre, University of Reading, UK

Bergen, 28 May 2018

slide-2
SLIDE 2

Outline

◮ particle filter, sample degeneracy ◮ equal weights, implicit sampling ◮ implicit equal-weights particle filter

slide-3
SLIDE 3

Particle filter

Probability density p(x) represented by weighted ensemble {xi, wi}Ne

i=1

ˆ p(x) =

  • i

wiδ(x − xi) E[g(X)] =

  • g(x)p(x)dx ≈
  • g(x)ˆ

p(x)dx =

  • i

wig(xi) p(x|y) = p(y|x)p(x) p(y) ≈

  • i

w new

i

δ(x − xi) w new

i

= wi · p(y|xi) p(y)

slide-4
SLIDE 4

Importance Sampling (IS) filter

slide-5
SLIDE 5

Sequential Importance Resampling (SIR) filter

A.k.a. the bootstrap filter (Gordon et al., 1993)

slide-6
SLIDE 6

Importance sampling and optimal proposal density

◮ Draw samples from proposal distribution q(x) ◮ Correct weights for difference between q and p

w corrected

i

= p(xi) q(xi)wi Optimal proposal density

◮ Suppose target is filtering distribution at time tn: p(xn|y1:n) ◮ Then choosing q(xn) = p(xn|xn−1 i

, yn) minimizes Var(w n

i )

slide-7
SLIDE 7

Implicit sampling

Implicit Particle Filter (IPF) (Chorin and Tu, 2009)

◮ Want samples from p(x|y) ◮ ξ ∼ g(ξ) ◮ ψ maps mode of g(ξ) to mode of p(x|y)

ξ x ψ

◮ G(ξ) = − log g(ξ), F(x) = − log[p(x|xprev)p(y|x)] ◮ To find x given ξ, solve F(x) − ϕF = G(ξ) − ϕG ◮ w(xn) = p(xn|xn−1)p(y n|xn−1) g(ξ)

  • ∂xn

∂ξ

  • ∝ e−ϕF +ϕG
  • ∂xn

∂ξ

slide-8
SLIDE 8

Equal Weights

Force all the weights to be equal by construction w1 = w2 = . . . = wNe = wtarget

◮ Transformation ψ : ξ → xi involves parameter αi ◮ Weight wi is a function of αi ◮ Choose αi so that wi = wtarget

slide-9
SLIDE 9

Implicit equal-weights particle filter (IEWPF)

Zhu, van Leeuwen and Amezcua (2016)

ˆ xn

i

xn

i

α1/2

i

P1/2ξi xn−1

i

M(xn−1

i

)

◮ xa i = arg maxx p(x|xn−1 i

, y n)

◮ ξ ∼ q(ξ) ◮ Weight of particle i:

wi = w prev

i

· p(xn

i |xn−1 i

, y n)p(y n|xn−1

i

) q(ξ)

  • ∂x

∂ξ

  • = wtarget

◮ Solve for αi to determine xn i for i = 1, . . . , Ne

slide-10
SLIDE 10

Transformation from ξ to x

0.97 0.975 0.98 0.985 0.99 0.995 1 0.97 0.975 0.98 0.985 0.99 0.995 1

g = ξTξ, b = αg

slide-11
SLIDE 11

Gauss-linear test case

250 300 350 400 Time 4 6 8 10 12 x1

Truth Observations IEWPF ensemble

xn = xn−1 + ηn−1 yn = xn

truth + ǫn

Nx = 1000 Ne = 25 ηn−1 ∼ N(0, Q), ǫn ∼ N(0, R), x0 ∼ N(0, B)

slide-12
SLIDE 12

Gauss-linear test case: Ensemble variance over time

20 40 60 80 100 120

Time step

0.1 0.2 0.3 0.4 0.5 0.6

Variance

KF IEWPF

slide-13
SLIDE 13

Gauss-linear test case: Ensemble variance final distribution

slide-14
SLIDE 14

Two-stage IEWPF

Two-stage update scheme: xn

i = xa i + β1/2P1/2ηi + α1/2 i

P1/2ξi ξT

i ηi = 0

xa

i

xn+1

i

α1/2

i

P1/2ξi xn

i

M(xn

i )

β1/2P1/2ηi x∗

i

slide-15
SLIDE 15

Two-stage IEWPF: Ensemble variance over time

20 40 60 80 100 120

Time step

0.1 0.2 0.3 0.4 0.5 0.6

Variance

β = 0.05 β = 0.25 β = 0.30 β = 0.50 KF

slide-16
SLIDE 16

Two-stage Ensemble Variance Final Distribution

slide-17
SLIDE 17

Single-stage IEWPF rank distribution

x(1) ≤ x(2) ≤ . . . ≤ x(r) ≤ xtrue ≤ x(r+1) ≤ . . . ≤ x(Ne) (should be more or less uniform)

slide-18
SLIDE 18

Two-stage IEWPF rank distribution

10 20

Rank

50 100

Frequency β = 0.05

10 20

Rank

50 100

Frequency β = 0.25

10 20

Rank

50 100

Frequency β = 0.30

10 20

Rank

50 100

Frequency β = 0.50

slide-19
SLIDE 19

Non-linear test case

Lorenz96 model with Nx = 40, Ny = 20, Ne = 100 dxi dt = −xi−2xi−1 + xi−1xi+1 − xi + F, i = 1, . . . , Nx.

20 40 60 80 100

Time steps

  • 10
  • 5

5 10 15

x1

xn = M(xn−1) + ηn−1, ηn−1 ∼ N(0, Q) ym = Hxm + ǫn, ǫn ∼ N(0, R)

slide-20
SLIDE 20

Non-linear test case

  • 20

20

x1

Ensemble Truth

  • 20

20

x2

Ensemble Truth 50 100 150 200 250 300 350 400

Time step

50

Variance

x1 x2

slide-21
SLIDE 21

Non-linear test case: Coverage probability

0.5 1 1.5

β

0.4 0.5 0.6 0.7 0.8 0.9

Average coverage probability Coverage prob. of 80 pct prediction interval

slide-22
SLIDE 22

Non-linear test case: Calibration

slide-23
SLIDE 23

Non-linear test case: Rank distribution

slide-24
SLIDE 24

Conclusion

◮ IEWPF ensures equal weights, prevents ensemble degeneracy, but

underestimates variance

◮ Two-stage scheme is able to achieve correct variance, but adds a

tuning parameter Properties under study

◮ Choice of target weight affects quality of estimates ◮ Setting target weight too large means some particles must get lower

weights

◮ Setting target weight low enough for all weights to be equal induces

a bias