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 Bergen, 28 May 2018

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

Particle filter Probability density p ( x ) represented by weighted ensemble { x i , w i } N e i =1 � p ( x ) = ˆ w i δ ( x − x i ) i � � � E [ g ( X )] = g ( x ) p ( x ) dx ≈ g ( x )ˆ p ( x ) dx = w i g ( x i ) i p ( x | y ) = p ( y | x ) p ( x ) � w new ≈ δ ( x − x i ) i p ( y ) i = w i · p ( y | x i ) w new i p ( y )

Importance Sampling (IS) filter

Sequential Importance Resampling (SIR) filter A.k.a. the bootstrap filter (Gordon et al., 1993)

Importance sampling and optimal proposal density ◮ Draw samples from proposal distribution q ( x ) ◮ Correct weights for difference between q and p = p ( x i ) w corrected q ( x i ) w i i Optimal proposal density ◮ Suppose target is filtering distribution at time t n : p ( x n | y 1: n ) ◮ Then choosing q ( x n ) = p ( x n | x n − 1 , y n ) minimizes Var( w n i ) i

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 | x prev ) p ( y | x )] ◮ To find x given ξ , solve F ( x ) − ϕ F = G ( ξ ) − ϕ G ◮ w ( x n ) = p ( x n | x n − 1 ) p ( y n | x n − 1 ) � � � � � ∂ x n � ∂ x n � ∝ e − ϕ F + ϕ G � � � � g ( ξ ) ∂ξ ∂ξ �

Equal Weights Force all the weights to be equal by construction w 1 = w 2 = . . . = w N e = w target ◮ Transformation ψ : ξ �→ x i involves parameter α i ◮ Weight w i is a function of α i ◮ Choose α i so that w i = w target

Implicit equal-weights particle filter (IEWPF) Zhu, van Leeuwen and Amezcua (2016) x n − 1 x n i i α 1 / 2 P 1 / 2 ξ i i M ( x n − 1 x n ) ˆ i i i = arg max x p ( x | x n − 1 ◮ x a , y n ) i ◮ ξ ∼ q ( ξ ) ◮ Weight of particle i : i | x n − 1 , y n ) p ( y n | x n − 1 · p ( x n ) � ∂ x � w i = w prev i i � � � = w target i � � q ( ξ ) ∂ξ � ◮ Solve for α i to determine x n i for i = 1 , . . . , N e

Transformation from ξ to x 1 0.995 0.99 0.985 0.98 0.975 0.97 0.97 0.975 0.98 0.985 0.99 0.995 1 g = ξ T ξ, b = α g

Gauss-linear test case x n = x n − 1 + η n − 1 12 y n = x n truth + ǫ n 10 x 1 8 N x = 1000 Truth 6 Observations IEWPF ensemble N e = 25 4 250 300 350 400 Time η n − 1 ∼ N (0 , Q ) , ǫ n ∼ N (0 , R ) , x 0 ∼ N (0 , B )

Gauss-linear test case: Ensemble variance over time 0.6 KF IEWPF 0.5 0.4 Variance 0.3 0.2 0.1 0 0 20 40 60 80 100 120 Time step

Gauss-linear test case: Ensemble variance final distribution

Two-stage IEWPF Two-stage update scheme: i + β 1 / 2 P 1 / 2 η i + α 1 / 2 i = x a P 1 / 2 ξ i x n i ξ T i η i = 0 x n +1 α 1 / 2 i P 1 / 2 ξ i i x ∗ x n i i β 1 / 2 P 1 / 2 η i M ( x n x a i ) i

Two-stage IEWPF: Ensemble variance over time 0.6 β = 0 . 05 β = 0 . 25 β = 0 . 30 β = 0 . 50 0.5 KF 0.4 Variance 0.3 0.2 0.1 0 0 20 40 60 80 100 120 Time step

Two-stage Ensemble Variance Final Distribution

Single-stage IEWPF rank distribution x (1) ≤ x (2) ≤ . . . ≤ x ( r ) ≤ x true ≤ x ( r +1) ≤ . . . ≤ x ( N e ) (should be more or less uniform)

Two-stage IEWPF rank distribution β = 0 . 05 β = 0 . 25 100 100 Frequency Frequency 50 50 0 0 0 10 20 0 10 20 Rank Rank β = 0 . 30 β = 0 . 50 100 100 Frequency Frequency 50 50 0 0 0 10 20 0 10 20 Rank Rank

Non-linear test case Lorenz96 model with N x = 40, N y = 20, N e = 100 d x i d t = − x i − 2 x i − 1 + x i − 1 x i +1 − x i + F , i = 1 , . . . , N x . 15 10 5 x 1 0 -5 -10 0 20 40 60 80 100 Time steps x n = M ( x n − 1 ) + η n − 1 , η n − 1 ∼ N (0 , Q ) y m = Hx m + ǫ n , ǫ n ∼ N (0 , R )

Non-linear test case 20 Ensemble Truth x 1 0 -20 20 Ensemble x 2 Truth 0 -20 Variance 50 x 1 x 2 0 0 50 100 150 200 250 300 350 400 Time step

Non-linear test case: Coverage probability Coverage prob. of 80 pct prediction interval 0.9 Average coverage probability 0.8 0.7 0.6 0.5 0.4 0 0.5 1 1.5 β

Non-linear test case: Calibration

Non-linear test case: Rank distribution

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