Kalman Smoother and EM for Stochastic Constant Growth Model

by allenlu2007

 

以下的例子為 diffusion approximation for a stochastic exponential growth model.  可以用 linear space model 

NewImage

NewImage

 B 是平均 growth rate.  A = C = 1.

NewImage

NewImage

Initial condition t=1

NewImage


再來結合 Markov 的 property (1)

NewImage

 

and observed output only depends on the corresponding states (2)

NewImage

 

We can get a very simple joint pdf of all Xi and Yi in Kalman filter (or Hidden Markov Model) as!!

NewImage

以下是有趣的部份:

NewImage

Where only Yi are given. 

記住 Xi 是 hidden variables.  因些我們可以(需要)用 EM algorithm 來 optimize

 

EM Algorithm

由 Shumway and Stoffer (1982) 所建議,基本上有四個步驟。

(0) 計算 initial parameter estimates

NewImage

(1) Kalman smoother: generate an estimate of xt by estimating E(xt | Y1, … YT) given above initial values.

(2) Define  Q function as below.  Update R, Q, B, Pi, V by maximizing the Q function.

NewImage

(3) Check if si function is converged and repeat (1)

 

以下為分解動作

Step 0: Initial Parameter Estimates

By Holmes and Fagan (2002)

NewImage

NewImage

 

Step 1: Kalman Smoother (Kalman-Rauch Recursion)

NewImage

(a)  先用 Kalman filter forward pass to estimate Gaussian pdf of p(xt | y1, .. yt).  

   Gaussian pdf 由 mean and variance 完全決定。Mean (xt = E(xt | y1, .. yt) 和 variance (Vt) 結果如下

NewImage

   Vt 在 prediction phase 時會增加 (by Q), 而在 update phase 時會減小 (by Kt Vt^(t-1))

(b) 再用 Rauch backward pass to estimate p(xt| y1, .. yT).  同樣 mean (xt) and variance (Vt) 結果如下

NewImage

   最後一項的右手第二項一般是負的,所以整體的 variance 可以再減少。 

(c) 再一個 backward pass 是計算

NewImage   

NewImage

 

Step 2: EM Algorithm 

接下來定義 E-step and M-step.  E-step 如下

NewImage

利用在開始所推導的 log P 公式,其實我們不必算出 psi function.

重點在於 taking partial derivative of psi w.r.t. each parameters, set the derivative to zero and solve the maximizing parameter:

NewImage

(好像少了一個 1/T log |R| 的導數? and what is C?)

NewImage

NewImage

Step 3: 檢查收歛性 

 

Matlab code example

clear all

close all

 

A = 1

C = 1

B = 0.2

D = 0

 

t = [0:200];

Q = 0.02;

R = 1;

 

x0 = 15;

 

sys = ss(A, B, C, D, -1);

 

u = ones(length(t), 1);     % external force B*u

v = sqrt(Q)*randn(length(t), 1); %state noise

w = sqrt(R)*randn(length(t), 1); %measurement noise

 

y = lsim(sys, u+v, t, x0); 

logdata = y + w;

 

plot(logdata)

 

[BB, QQ, RR, initx, V1, LL] = PVAKalman(logdata’)

 

BB = 0.2

RR = 0.9748

QQ = 0.0013

initx = 15.1161

 

NewImage

 

EM Algorithm 2:

考慮另一個 model

NewImage

A, B, and C 是已知,但是 Q and R 未知:

NewImage

 

 

Advertisements