漫談 HMM之三:Kalman/Particle Filtering

NO IMAGE

漫談 HMM之三:Kalman/Particle Filtering

上次我們講了 HMM 的 Forward-Backward 演算法,得到了關於 α\alpha 和 β\beta 的遞推公式。不過由於中間需要進行 marginalization,這些式子裡有麻煩的積分存在。如果是離散型的隨機變數,積分實際上是求和,一般來說就沒有什麼問題了。但是對於上一次舉的火星車移動的例子,實際上隨機變數是連續的,比如火星車的位置 zz 實際上是一個 R3\mathbb{R}^3(或者其他任何用於火星表面的座標系)下的連續隨機變數,於是我們就面臨了求積分的問題。

積分是一件很困難的事情,就算是裝備了 Wolfram Mathematica 的 PhD 也不能隨便誇海口的,因為有很多積分式本來就沒有一個式子可以把結果寫下來。當然你可以像 Γ\Gamma-函式 那樣自己定義新的符號來隱藏後面算不出來的一大坨式子,不過這對於最終需要做數值計算來說好像是躲得了和尚躲不了廟啊。所以我們就不說這些歪門邪道了,正統的解決方法主要有兩種。第一種是把概率分佈限制到高斯分佈上,由於高斯分佈像開了外掛似得,求 marginalization 和 conditioning 之後都還是高斯分佈,並且本身只需要用均值和方差兩 “個” 量即可描述,所以我們在這裡碰到的積分問題都可以解析地計算出來,這裡介紹一種代表性的演算法叫做 Kalman Filtering;另外一種方法是放棄精確計算,採用近似計算的方法,這裡也介紹一種代表性方法叫做 Particle Filtering。

於是一切要從高斯分佈說起,高斯分佈有幾個很好的性質:

一個服從多維聯合高斯分佈的隨機變數 XX 的線性投影 AXAX 仍然服從高斯分佈
一個服從多維聯合高斯分佈的隨機變數 XX 在對其中一部分維度進行 marginalization 和 conditioning 的時候,結果仍然是多維聯合高斯分佈

並且這些結果的高斯分佈的引數可以通過 XX 本身的分佈引數帶入公式直接得到,這裡就不一個一個寫出來了,感興趣的同學可以自行推導或者參考 Wikipedia。於是,我們回顧一下上次推導得到的 Forward-Backward 演算法裡關於 α\alpha 的遞推式

α(zt 1)=∫ztα(zt)p(zt 1|zt)dztp(xt 1|zt 1) \alpha(z_{t 1}) = \int_{z_t}\alpha(z_t)p(z_{t 1}|z_t)\;dz_t\;p(x_{t 1}|z_{t 1})

以及 β\beta 的遞推式:

β(zt)=∫zt 1p(xt 1|zt 1)β(zt 1)p(zt 1|zt)dzt 1 \beta(z_t) = \int_{z_{t 1}} p(x_{t 1}|z_{t 1})\beta(z_{t 1})p(z_{t 1}|z_t)\;dz_{t 1}

為了用到高斯分佈的性質來求積分,我們需要把積分式裡的那些東西都限制為高斯分佈,其中主要就是有狀態轉移 p(zt 1|zt)p(z_{t 1}|z_t) 和觀察值 p(xt|zt)p(x_t|z_t) 兩種情況。當然使用高斯分佈的原因是除了高斯分佈之外都沒法算,但是這個理由多少有點不夠響亮,所以我們得找點其他的 justification:看看用高斯分佈來進行建模到底有沒有道理。

具體來說,對於火星車的例子,把 p(xt|zt)p(x_t|z_t) 建模成一個高斯分佈實際上是很自然的,因為在給定真實值 zz 的情況下,測量誤差一般會被建模為一個均值為零,方差為 Σ\Sigma 的高斯分佈,理由當然可以從中心極限定理啊或者各種經驗啊之類的扯一大堆,總之這是一個標準做法,這樣一來實際上觀察值的條件分佈就是

p(xt|zt)=N(zt,Σ) p(x_t|z_t) = \mathcal{N}(z_t,\Sigma)

接下來是狀態轉移,如果 tt 時刻的狀態 ztz_t 已經是服從高斯分佈的了(當然從最初的其實狀態 z0z_0 是可以由我們設定的),如何使得 zt 1=f(zt)z_{t 1}=f(z_t) 也服從高斯分佈呢?根據高斯分佈隨機變數的性質,如果 ff 是一個線性函式,亦即存在矩陣 AA 使得 f(zt)=Aztf(z_t) = Az_t,那麼 zt 1z_{t 1} 仍然是一個高斯分佈。當然,考慮到之前提過的機械裝置本身的操作誤差,我們並不能完美地得到想要的結果,所以實際上 ff 是這樣子的:f(zt)=Azt ϵf(z_t)=Az_t \epsilon,這裡 ϵ\epsilon 是一個獨立的誤差隨機變數,我們再一次將 ϵ\epsilon 也建模為一個零均值 Σ′\Sigma’ 方差的高斯分佈,原因和剛才一樣:高斯分佈一向被用來建模誤差,更重要的是,如此一來,zt 1z_{t 1} 還是一個高斯分佈。並且,在這種情況下:

p(zt 1|zt)=N(Azt,Σ′) p(z_{t 1}|z_t) = \mathcal{N}(Az_t,\Sigma’)

值得注意的一點是,這裡 ff 是線性的這一點非常重要,如果 ff 是一個 general 的非線性函式,即使能保證 zt 1|ztz_{t 1}|z_t 這個條件分佈是高斯的,也沒法做到 marginal zt 1z_{t 1} 是高斯的,那樣的話 Forward-Backward 的遞推就沒有辦法走下去。

所以,總結起來,把一切都建模成高斯分佈的做法實際上在各個方面都是相當自然和有道理的,唯一的一個缺陷就是狀態轉移這裡被限制成了只能是線性函式。這也是為什麼 Kalman Filter 被稱為是線性方法。這在有些時候是一個非常大的限制,所以不得不再探索其他的方法和擴充套件,其中有一個擴充套件叫做 Extended Kalman Filter,是將非線性的狀態轉移進行線性化近似,據說這個擴充套件演算法被用在了阿波羅登月計劃裡。不過就最土的 Kalman Filter 本身已經在非常多的問題中得到了廣泛和成功的應用了。

說了半天,其實還沒有講 Kalman Filter 到底是什麼,其實 Kalman Filter 就是在我們剛才的描述下將所有的分佈全部取成高斯的情況下的連續隨機變數 HMM 的 Forward-Backward 演算法……的變種。之所以說是一種變種,是因為它並不是計算 α\alpha 和 β\beta,而是計算了 α\alpha 和一個叫做 γ\gamma 的東西。其中 α\alpha 的定義和之前一樣的,而 γ\gamma 其實也不陌生:

γ(zt)≜p(zt|xT0)=α(zt)β(zt)p(xT0) \gamma(z_t) \triangleq p(z_t|x_0^T) = \frac{\alpha(z_t)\beta(z_t)}{p(x_0^T)}

也可以推匯出一個像 β\beta 一樣的從後往前的遞推公式。於是整個演算法還是 Forward 和 Backward 兩輪,Forward 的時候是計算 α(zt)=p(zt|xt0)\alpha(z_t)=p(z_t|x_0^t),也就是根據到目前為止的觀測資料所能得到的對於狀態值 ztz_t 的分佈估計,其中每一步迭代又被分為兩步完成,第一步計算 p(zt|xt−10)p(z_t|x_0^{t-1}) 叫做 “Prediction” ,也就是直接根據狀態轉移預測從 t−1t-1 時刻到 tt 時刻之後可能所處的狀態了,接下來會計算 p(zt|xt0)p(z_t|x_0^t),也就是加入 tt 時刻的測量值 xtx_t,這一步叫做 “Update” 。整個正向的迭代合在一起叫做 “Filtering” 。而 Backward 的迭代則是根據所有時刻 0,…,T0,\ldots,T 的觀察值來對 ztz_t 時刻的狀態分佈的計算進行修正,也就是計算 p(zt|xT0)p(z_t|x_0^T),這個步驟也有個名字,叫做 “Smoothing” 。

那麼,如果是任意的連續的隨機變數,到底應該怎麼做呢?既然精確計算不行了,那麼久只有求助於近似計算。首先面臨的一個問題是如何去表示一個分佈,最直接的近似表達一個分佈的方法莫過於所謂的 empirical distribution 了:如果 x1,…,xNx_1,\ldots,x_N 是服從於 PxP_x 的 IID 樣本 1 1 突然想起之前有過一次有趣的關於 “樣本” 是什麼的討論,結論是,所謂 xx 是分佈 PxP_x 的一個樣本,這句話在數學上的意義其實就是在說 xx 是個隨機變數並且它的分佈是 PxP_x。 的話,那麼

P^x(⋅)=1N∑Ni=11xi(⋅) \hat{P}_x(\cdot) = \frac{1}{N}\sum_{i=1}^N \mathbf{1}_{x_i}(\cdot)

就是關於原來分佈的一個近似(如果是概率密度函式的話則要使用 狄拉克 δ\delta 函式 ),特別地,關於原來分佈的期望可以通過這個 empirical distribution 近似為

Ex[f(x)]≈1N∑Ni=1f(xi) \mathbb{E}_{x}[f(x)] \approx \frac{1}{N}\sum_{i=1}^N f(x_i)

根據強大數定理,當 N→∞N\rightarrow \infty 的時候右邊會以概率 1 收斂到左邊。也就是說,一個分佈可以由一堆 sample 來近似表示,不過這裡我們要用一個改良的表示方法,除了 sample 之外,每個 sample xix_i 還被賦予一個權值 wiw_i。具體來說,這個東西是來自於 Importance Sampling——本身也是一個非常重要的 sampling 方法,它可以在不知道 partition function 的情況下對一個分佈進行採用。

具體來說,如果一個分佈

p(x)=p∗(x)Zp p(x) = \frac{p^*(x)}{Z_p}

其中 ZpZ_p 是 partition function(或者叫做 normalization factor),在概率圖模型中經常都會碰到 ZpZ_p 很難算的情況,這樣取樣也會變得很困難,不過 Importance Sampling 可以解決這個問題。具體的辦法是先找另一個比較容易採用的分佈 qq(比如是高斯分佈或者平均分佈之類的),並從 qq 那裡採 NN 個樣本 x1,…,xNx_1,\ldots,x_N,然後對每個樣本定義權值

wi=p∗(xi)q(xi) w_i = \frac{p^*(x_i)}{q(x_i)}

在對期望進行近似的時候,不再使用簡單的平均而是加權平均:

1N∑Ni=1wif(xi)1N∑Ni=1wi→Ep[f(x)]a.s.,as N→∞ \frac{\frac{1}{N}\sum_{i=1}^N w_if(x_i)}{\frac{1}{N}\sum_{i=1}^N w_i} \rightarrow \mathbb{E}_p[f(x)]\;a.s., \quad \text{as } N\rightarrow\infty

這裡的收斂性同樣是根據大數定理得到的,對分子分母分別求期望就可以很容易得到想要的結果 2 2 注意這裡必須在強大數定理的那種 “以概率 1 收斂” 意義下才能分子分母分別求期望之後再相除,只用弱大數定理似乎不足以得到這樣的結論。感謝 XH 同學提供的參考。 。而整個過程有效地避開了對 ZpZ_p 的計算。當然 Importance Sampling 也並不是隨便抓一個平均分佈就可以解決世間一切問題的銀彈,對取樣分佈 qq 的選取有時候會對近似的好壞收斂速率等產生非常大的影響。

而把這裡的權重和樣本組成的對叫做 Particle,就成了 Particle Filter 演算法,其整體結構和 Kalman Filter 如出一轍,其中 Prediction 是從 p(zt|xt0)p(z_t|x_0^t) 計算 p(zt 1|xt0)p(z_{t 1}|x_0^t),換句話說,我們已經有了 p(zt|xt0)p(z_t|x_0^t) 分佈的 particle {zti,wti}Ni=1\{z^t_i,w^t_i\}_{i=1}^N,現在要通過這些 particle 來生成 p(zt 1|xt0)p(z_{t 1}|x_0^t) 的 particle。由於 HMM 結構的特殊性,我們可以令每個 zt 1iz^{t 1}_i 是由 p(zt 1|zti)p(z_{t 1}|z_i^t) 中取樣出來的一個 sample。我們假設這個步驟是容易實現的,否則可以再用一次 Importance Sampling 或者用其他的 Metropolis-Hastings 之類的近似取樣方法。而每個樣本的權重則保持不變 w~t 1i=wti\tilde{w}^{t 1}_i = w^t_i。

假設原來的 ztiz^t_i 本身是取樣自分佈 qq,而權重遞迴地等於 wti=p(zti|xt0)/q(zti)w^t_i=p(z^t_i|x_0^t)/q(z^t_i)。這裡上下標有點混亂了…… -.-bb tt 是用於表示時刻的量,而 ii 則用於對 particle 進行計數。於是我們有

E[w~t 1i]E[w~t 1if(zt 1i)]=∫w~t 1iq(zti)dzti=∫p(zti|xt0)dzti=1=E[E[w~t 1if(zt 1i)|zti]]=∫∫w~t 1if(zt 1i)p(zt 1i|zti)q(zti)dzt 1idzti=∫f(zt 1i)∫p(zti|xt0)p(zt 1i|zti)dztidzt 1i=∫f(zt 1i)p(zt 1i|xt0)dzt 1i=E[f(zt 1i)|xt0] \begin{aligned} \mathbb{E}[\tilde{w}^{t 1}_i] &= \int \tilde{w}^{t 1}_i q(z^t_i) \;d z^t_i \\ &= \int p(z_i^t|x_0^t)\;dz_i^t \\ &= 1 \\ \mathbb{E}[\tilde{w}_i^{t 1}f(z_i^{t 1})] &= \mathbb{E}\left[\mathbb{E}[\tilde{w}_i^{t 1}f(z_i^{t 1})|z^t_i]\right] \\ &= \int\int \tilde{w}_i^{t 1}f(z_i^{t 1}) p(z^{t 1}_i|z^t_i)q(z^t_i)\;dz^{t 1}_idz^t_i \\ &= \int f(z^{t 1}_i) {\color{red}{\int p(z^t_i|x_0^t)p(z^{t 1}_i|z_i^t)\;dz^t_i}}\;dz^{t 1}_i \\ &= \int f(z_i^{t 1}) {\color{red}{p(z_i^{t 1}|x_0^t)}}\;dz_i^{t 1} \\ &= \mathbb{E}[f(z^{t 1}_i)|x_0^t] \end{aligned}

也就是說 particle zt 1i,w~t 1iz_i^{t 1}, \tilde{w}^{t 1}_i 確實是在對分佈 p(zt 1|xt0)p(z_{t 1}|x_0^t) 在做近似。接下來是 Update step,這次的 particle 是保持 zt 1iz^{t 1}_i 不變,而修改權重

wt 1i=w~t 1ip(xi 1|zt 1i)=wtip(xi 1|zt 1i) w^{t 1}_i = \tilde{w}^{t 1}_i p(x_{i 1}|z^{t 1}_i) = w^t_i p(x_{i 1}|z^{t 1}_i)

和剛才一樣我們可以證明對這樣的 particle 加權平均其實是在對 p(zt 1i|xt 10)p(z^{t 1}_i|x_0^{t 1}) 在做近似,具體推導就不專門寫出來了。總而言之這樣不斷地迭代就可以一直得到關於 p(zt|xt0)p(z_t|x_0^t) 的一個近似表示,而又不受限於狀態轉移的具體形式——當然,必須也要使得對 p(zt 1|zti)p(z_{t 1}|z_i^t) 的工作可以完成才行。

結束之前再提一句的是,隨著迭代的進行有些 particle 的權重可能會變得非常非常小,所以有時候需要進行一下 resampling。而用 particle 對分佈進行近似這個思想當然不止適用於 HMM,而是可以用在 general 的 graphical model 的 inference 上,但是對於 general 的 graph,就沒有 HMM 那麼良好的結構有上面的巧妙的演算法可以直接得到迭代下一步的 particle 了,而是需要顯式地把幾個分佈相乘再 marginalize 然後取樣新的 particle,其中每一步都並不是 trivial 的,比如 “兩個由 particle 表示的分佈相乘之後究竟應該會是什麼樣的” 。還有就是 particle filter 屬於 Sequential Monte Carlo 方法的一種。

至於為什麼這些演算法都叫做什麼什麼 Filter 嘛……其實我也不知道,估計是因為是一步一步迭代的,並且每次輸入 xtx_t 會得到 p(zt|xt0)p(z_t|x_0^t) 所以看起來像一個 filter 一樣的?

轉自:http://freemind.pluskid.org/machine-learning/hmm-kalman-particle-filtering/