Short Course State Space Models, Generalized Dynamic Systems and - PowerPoint PPT Presentation

Short Course State Space Models, Generalized Dynamic Systems and Sequential Monte Carlo Methods, and their applications in Engineering, Bioinformatics and Finance Rong Chen Rutgers University Peking University 1

Part Three: Advanced Sequential Monte Carlo 3.1 Mixture Kalman Filter 3.1.1 Conditional Dynamic Linear Models 3.1.2 Mixture Kalman Filters 3.1.3 Partial Conditional Dynamic Linear Models 3.1.4 Extend Mixture Kalman Filters 3.1.5 Future Directions 3.2 Constrained SMC 3.3 Parametrer Estimation in SMC 2

3.1.1 Conditional Dynamic Linear Models Indicator Λ t : (unobserved) If Λ t = λ then x t = H λ x t − 1 + W λ w t y t = G λ x t + V λ v t where w t ∼ N (0 , I ) and v t ∼ N (0 , I ) and independent. Given the trajectory of the indicator { Λ 1 , . . . Λ t } , the system is linear and Gaussian. 3

Example: Tracking a target in clutter Introducing an indicator I t taking values in { 0 , 1 , . . . , m t } . I t = 0 true signal missing. I t = i , then y ( i ) is the true signal. t x t = Hx t − 1 + Ww t y ( i ) = Gx t + V v t if I t = i t y ( j ) ∼ Unif (∆) if I t � = j t and P ( I t = 0) = p d and P ( I t = i ) = (1 − p d ) /m t Given the trajectory of the indicator { I 1 , . . . I t } , the system is linear and Gaussian. 4

Example: Tracking a target with non-Gaussian innovations. x t = Hx t − 1 + Ww t y t = Gx t + V v t where w t ∼ t k 1 , v t ∼ t k 2 . � χ 2 Note that t k = N (0 , 1) / k /k Introducing indicators Λ t = (Λ t 1 , Λ t 2 ) . √ k 1 x t = Hx t − 1 + √ λ 1 Ww t if Λ t 1 = λ 1 √ k 2 y t = Gx t + √ λ 2 V v t if Λ t 2 = λ 2 with v t ∼ N (0 , I ) , w t ∼ N (0 , I ) and Λ t 1 ∼ χ 2 k 1 , Λ t 2 ∼ χ 2 k 2 . Given the trajectory of the indicator { Λ 1 , . . . Λ t } , the system is linear and Gaussian. 5

Example: Tracking a target with random (Gaussian) acceleration plus maneuvering x t = Hx t − 1 + Fs I t u t + Ww t y t = Gx t + V v t where u t , w t and v t are all N (0 , I ) independent. I t maneuvering status: I t = 0 , no maneuvering, s 0 = 0 I t = 1 , slow maneuvering, s 1 I t = 2 , fast maneuvering, s 2 With known transition matrix P = P ( I t +1 | I t ) . 6

3.1.2 Mixture Kalman Filter: Let y t = ( y 1 , . . . , y t ) and Λ t = (Λ 1 , . . . , Λ t ) . Note that � p ( x t | y t )= p ( x t | Λ t , y t ) dF ( Λ t | y t ) and p ( x t | Λ t , y t ) ∼ N ( µ t ( Λ t ) , σ 2 t ( Λ t )) where KF t ( Λ t ) ≡ ( µ t ( Λ t ) , σ 2 t ( Λ t )) can be obtained from Kalman filter. p ( x t | y 1 , . . . y t ) is a mixture Gaussian distribution. 7

(Sequential) Monte Carlo Filter: a discrete sample with weight { ( x (1) t , w (1) t ) , . . . , ( x ( m ) , w ( m ) ) } = ⇒ p ( x t | y 1 , . . . , y t ) t t Mixture Kalman Filter: a discrete sample with weight { ( λ (1) t , w (1) t ) , . . . , ( λ ( m ) , w ( m ) ) } = ⇒ p ( Λ | y 1 , . . . , y t ) t t and a random mixture of Normal distributions m w ( i ) t N ( µ t ( λ ( i ) t ( λ ( i ) � t ) , σ 2 t )) = ⇒ p ( x t | y 1 , . . . , y t ) . i =1 Hence m � w ( i ) f ( x ) φ ( x ; µ t ( λ ( i ) t ( λ ( i ) � t ) , σ 2 E ( f ( x t ) | y 1 , . . . , y t ) ≈ t )) dx t i =1 8

Benefit: improved efficiency V ar [ f ( x t ) | y t ] ≥ V ar [ E ( f ( x t ) | Λ t , y t ) | y t ] Example: X ∼ N (Λ , σ 2 1 ) and Λ ∼ N (0 , σ 2 2 ) . Estimate µ = E ( X ) (1) directly sample from X ∼ N (0 , σ 2 1 + σ 2 2 ) , �� m = σ 2 1 + σ 2 i =1 X i � 2 V ar (ˆ µ ) = V ar m m (2) sample Λ ∼ N (0 , σ 2 2 ) . � m � m i =1 E ( X | Λ i ) i =1 Λ i µ = ˆ = m m �� m = σ 2 i =1 Λ i � 2 V ar (ˆ µ ) = V ar m m 9

Algorithm: At time t , we have a sample ( λ ( i ) t , KF ( i ) t , w ( i ) t ) For t + 1 , (1) : generate λ ( i ) t +1 from a trial distribution g (Λ t +1 | λ ( i ) t , KF ( i ) t , y t +1 ) (2) : run one step Kalman filter conditioning on ( λ ( i ) t +1 , KF ( i ) t , y t +1 ) and obtain KF ( i ) t +1 . (3) : calculate the incremental weight p ( λ ( i ) t , λ ( i ) t +1 | y t +1 ) u ( i ) t +1 = p ( λ ( i ) | y t ) g ( λ t +1 | λ ( i ) t , KF ( i ) t , y t +1 ) t and the new weight w ( i ) t +1 = w ( i ) t u ( i ) t +1 . 10

t t + 1 t + 1 p (Λ t | y t ) g (Λ t +1 | y t ,y t +1 ) p (Λ t +1 | y t , y t +1 ) ց ր ց ր g t ( · | Λ t , y t +1 ) w t +1 = w t u t +1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ( λ (1) F (1) t ,w (1) ( λ (1) F (1) ( λ (1) F (1) t +1 ,w (1) t ,K t ) − → t +1 , K t +1 ) − → t +1 ,K t +1 ) . . . . . . . . . ( λ ( m ) F ( m ) , w ( m ) ( λ ( m ) F ( m ) → ( λ ( m ) F ( m ) t +1 , w ( m ) , K ) − → t +1 , K t +1 ) − t +1 , K t +1 ) t t t 11

When Λ t is a discrete r.v. on a finite set, then (0) : For each j = 1 , . . . , J , run a Kalman filter to obtain u ( i ) j = p ( y t +1 | Λ t +1 = j, KF ( i ) t ) p (Λ t +1 = j | λ ( i ) t ) (1) : Sample a λ ( i ) t +1 from the set { 1 , . . . , J } with probability pro- portional to u ( i ) j . i.e. sample a λ from p (Λ t +1 | λ t , KF t , y t +1 ) (2) : Let KF ( i ) t +1 be the one with λ ( i ) t +1 . (3) : The new weight is J w ( i ) t +1 ∝ w ( i ) t p ( y t +1 | λ ( i ) t , KF ( i ) t ) ∝ w ( i ) u ( i ) � t j j =1 12

When Λ t is a continuous r.v., a simple (but not optimum) algo- rithm is (1) : Sample a λ ( i ) t +1 from p (Λ t +1 | Λ t = λ ( i ) t ) (2) : Run one step Kalman filter conditioning on ( λ ( i ) t +1 , KF ( i ) t , y t +1 ) and obtain KF ( i ) t +1 (3) : The new weight is w ( i ) t +1 = w ( i ) t p ( y t +1 | λ ( i ) t +1 , KF ( i ) t ) 13

Tracking in clutter: 40 20 error 0 -20 -40 0 20 40 60 80 100 time 40 20 error 0 -20 -40 0 20 40 60 80 100 time 40 20 error 0 -20 -40 0 20 40 60 80 100 time 14

Example: Tracking a target with non-Gaussian innovations State Equation: � x (1) � � � � x (1) � � � 1 T T/ 2 t t − 1 = + q w t x (2) x (2) 0 1 1 t t − 1 true signal y t = x (1) + rv t t where w t ∼ t 3 and v t ∼ t 3 . T = 1 , q 2 = 400 / 3 , r 2 = 1600 / 3 . 15

MSE x MSE speed 800 1500 M20 M20 M200 M200 P50 P50 P500 600 P500 1000 xx1[, 1] xx2[, 1] 400 500 200 0 0 0 20 40 60 80 100 0 20 40 60 80 100 Index Index 16

noise var # chains Particle MKF cpu time # miss cpu time # miss 20 9.49843 72 19.4277 1 50 20.1622 20 51.6061 1 q 2 = 16.0 200 80.3340 7 181.751 1 r 2 = 1600 500 273.369 4 500.157 1 1500 1063.36 3 2184.67 1 17

3.1.3 Partial CDLM state equation: x t = g t ( x t − 1 , ε t ) observation equation: y t = h t ( x t , e t ) • Extract the linear and Gaussian components out, and use Kalman Filter (integrating those components out) • Nonlinear components are dealt with standard Monte Carlo Filters • Non-Gaussian innovations are dealt with indicators and ap- proximations • Nonlinear functions are dealt with ’conditional linearization’? 18

State: ( x t , x ∗ t ) . Observations: ( y t , y ∗ t ) x t = g t ( x t − 1 , x ∗ t − 1 , ε t ) (1) x ∗ t = H x t x ∗ t − 1 + W x t w t (2) y t = h t ( x t , e t ) (3) y ∗ t = G x t x ∗ t + V x t v t (4) • x t , y t : nonlinear nonGaussian component • x ∗ t , y ∗ t : conditional linear Gaussian component • H x t , G x t , W x t , V x t : known matrices given x t • w t ∼ N (0 , I ) and v t ∼ N (0 , I ) and independent. Given the trajectory of the NLNG components { x 1 , . . . x t } , the system (2) (4) is linear and Gaussian. 19

Example: Digital Signal Extraction in Fading Channels  x t = Hx t − 1 + w t    State Equations: α t = Gx t  s t ∼ p ( · | s t − 1 )   Observation equation: y t = α t s t + v t Or State Equations: x t = Hx t − 1 + w t s t ∼ p ( · | s t − 1 ) Observation equation: y t = Gx t s t + v t 20

Example: 2-d target with GPS and IMU sensor. State: • position p 1 t , p 2 t • speed v 1 t , v 2 t • (total) acceleration a 1 t , a 2 t • IMU facing θ t • IMU rotational speed ψ t • Two motion status: – M t = 1 : (roughly) zero acceleration (constant between ob- servations) – M t = 2 : (roughly) constant acceleration 21

δ : time gap between observations State equations: ( M t = 1 ) p it = p it − 1 + v it − 1 δ + 0 . 5 δε it i = 1 , 2 v it = v it − 1 + ε it i = 1 , 2 a it = ε it i = 1 , 2 θ t = θ t − 1 + ψ t − 1 δ + 0 . 5 δε ∗ t ψ t = ψ t − 1 + ε ∗ t P ( M t = i | M t − 1 = j ) = p ij Similar for M t = 2 • One can also impose constraints (maps) on the state equations. • variance of ε it depends on platform (walking or vehicle etc) 22

Observations: • p ∗ 1 t , p ∗ 2 t : (post-processed) GPS signal • a ∗ 1 t , a ∗ 2 t : acceleration in the direction of θ t • η t : rotational acceleration Observational equations: p ∗ it = p it + e 1 t i = 1 , 2 a ∗ 1 t = cos ( θ t ) a 1 t + sin ( θ t ) a 2 t + w 1 t a ∗ 2 t = − sin ( θ t ) a 1 t + cos ( θ t ) a 2 t + w 2 t η t = ψ t − ψ t − 1 + w 3 t Give θ t , ψ t , M t , the rest of the system is linear and Gaussian. Hence, • θ t , ψ t , M t are the NLNG stat components • η t is the NLNG observation component. 23