卡爾曼濾波器詳解——從零開始(2) Kalman Filter from Zero

來源 · 知乎

作者 · 三腳貓Frank

編輯 · 睿小妹

原文 · https://zhuanlan.zhihu.com/p/179480833

我們繼續接上一篇講Kalman filter的基礎預備知識。把基礎部分擴展成三篇。這篇算第二篇，我們主要討論經典估計理論，下一篇我們主要討論貝葉斯估計理論和Wiener filter。值得一提的是，我們應該從後者的角度更多地去理解Kalman filter——從他1960年寫的關於KF的文章中，你能明顯找到Minimum Mean Square Error Estimator(MMSE)的使用，而這正是Bayesian estimation裡面一個常用的標準。而且狀態估計和參數估計不是一樣的，狀態作為被估計量，並不是一個固定值，也是一個隨機變量。

除了擴充了很多擴展參考資料供你們參考，Kalman本人的原文[6]被加入到了reference之中，感興趣或者著急的朋友可以直接去看他文章中是怎麼推導的。（一般都讀的懂，前提是你懂概率，統計和信號處理等。）不過你要是沒什麼基礎，就跟著文章一點點了解基礎知識吧。有很多方法和思路可以引導我們最終得到經典的Kalman filter方程，但是我會把它放到後面幾篇文章中去討論。在此之前，我希望帶領大家回顧和熟悉一些信號處理和經典的估計算法的內容作為知識的儲備。主要介紹一些思路，幫助理解，他日深入的研究就會容易。

本篇目錄：

1. Supplement to linear system theory and signal processing 線性系統和信號處理理論補充

2. Classic estimation theory 經典估計理論

2.1 Unbiasedness and consistency 無偏性與相合性

2.2 method of moments 矩法

2.3 Maximum Likelihood Estimation 最大似然估計(MLE)

2.4 Mean Square Error and Minimum-Variance Unbiased Estimator 均方差、最小方差無偏估計（MVUE）

2.5 Efficiency 有效性

2.6 Best Linear Unbiased Estimator 最佳線性無偏估計(BLUE)

2.7 Ordinary Least Squares 最小二乘算法 (OLS)

2.8 Weighted least squares 加權最小二乘算法 (WLS)

2.9 Recursive least squares 迭代最小二乘法 (RLS)

3. Summary 小結

Supplement to linear system thoery and signal processing 線性系統和信號處理理論補充 

1.1 信號處理知識補充

簡單一點的線性時不變系統：

這個系統的全解是以下面的式子來計算的：

前面一部分是初值影響的狀態量，後一部分是輸入引起的狀態量。我們之前定義過LTI系統的狀態轉移矩陣(state transition matrix)為  ，所以(2)也就被重寫成了：

我們以前也講過LTI的零狀態響應在時域中可以由以下convolution表示：

所以對比(3)和(4)，我們可以定義脈衝響應函數  . 注意卷積中  的積分變量  是可以和  的  互換而保持卷積值不變的。

現在假設系統的輸入是一個random process，統一把輸入改為寬平穩(Wide-sense stationary，WSS)隨機過程  ，因為實際中我們多數信號都滿足WSS要求。我們已經知道WSS隨機過程的均值是常數，自相關函數(autocorrelation function) 只與延時  有關，而與起始時間無關，且二階矩存在(為有限值)。兩個信號  被稱為jointly WSS, 如果它們都是WSS的，並且互相關函數(cross-correlation function)是時不變的，即只與信號之間的時間間隔有關。

註：我們再強調一下這個定義與stationary process的差別。前者要求更弱，並不要求二階矩以上的矩是否存在，只需要均值為常數，自相關函數隻與時間間隔有關；後者要求過程的概率密度函數，包括一維和聯合分布函數，都必須是嚴格時不變的，而對各階矩的存在性沒有任何要求。故而是平穩隨機過程的期望，方差等等是不一定存在的（如柯西分布），但是如果存在的話，應當也是時不變的。

知道了  的統計特性和功率譜函數，下面來計算  經過系統  後的輸出  的統計特性以及功率譜函數(power spectral density, PSD)：

我們考慮LTI系統的輸出

輸出的期望  滿足：

注意  在頻域中就是頻率特性或者傳遞函數  ，利用終值定理可知其穩態時的Gain為 

輸出  和輸入  的互相關函數(cross-correlation function)  :

註：不同定義中關於  的定義有略微不同，同理也可以定義  ，但本質上我們需要指定一個固定的公式，比如  ，那麼 。這裡假設了互相關函數僅與時間間隔有關。(7)中使用了  的定義。

由於  是WSS過程，我們知道其自相關函數  只與時間間隔  有關(這是個dummy variable，只要當做是一個變量決定了整個函數值就行了)。因此上式也就是：

 的自相關函數(autocorrelation function) 根據定義為：

我們發現：

借用(8)，我們知道：

根據維納-辛欽定理(Wiener–Khinchin theorem)，  都是WWS random process，我們把(8)兩端進行Fourier transformation，得到  與  的互功率譜密度函數(cross power spectral density, CPSD)  ：

根據(11) 我們就得到了  的自功率譜密度函數(power spectral density)的兩種表達式：

註： 不就是所謂的我們在經典控制理論中的頻率特性嗎？換言之，也就是傳遞函數，或者在system identification領域叫frequency response function, FRF。(12)表明，如果我們知道  的PSD,  ，以及頻率特性或者傳遞函數  ，那麼我們就可以得到CPSD  。同時這個結論還有一個非常重要的應用，那就是傳遞函數估計。如果能通過序列估計出CPSD和PSD，我們可以通過：
來估算傳遞函數，這種FRF的估計方法稱為  estimator。有  就有 ，即採用(13)來計算FRF，這兩種估計有不同的假設，即前者假設噪音與輸入  是uncorrelated，後者假設噪音與輸出  是uncorrelated。(14)我們放著後面再解釋。互譜和自譜的估計，我們在MATLAB中就有cpsd命令幫我們計算，而專門用於傳遞函數估計的tfestimate函數，也是調用了cpsd，以及採用了上面的公式作為其中一種估計方法。這是一種估計線性系統FRF的非常簡單又實用的辦法，只需要兩個時間序列  和  就可以了。

From MATLAB documentation - tfestimate

輸出Y(t)的平均功率，滿足各態歷經假說的前提下，由二階中心矩給出：

當我們把序列看成是一個n維實數向量  時，我們有自相關矩陣(auto-correlation matrix)是任意兩個時間  之間定義的矩陣：

自相關矩陣

而自協方差矩陣(autocovariance matrix)定義為：

這兩者的差別其實就在於autocovariance matrix減掉了  的期望。如果它的均值為0，即  ，那麼自相關矩陣和自協方差矩陣是等價的。至於為什麼要寫成矩陣的形式，那當然是因為矩陣運算有其獨特的優越性。同理我們可以兩個隨機過程之間的互相關矩陣和協方差矩陣，這個很容易類推就不贅述了。更多閱讀見[1,p116-119]。

1.2 隨機線性系統、白化處理

一個連續的線性隨機系統(linear stochastic differential equation)一般可以由以下方程表示：

這是一個時變線性系統，所以各矩陣都是時變的。其中  是測量輸出，與我們以前定義的  ，名義上的輸出，還不一樣（雖然這裡還是用  矩陣符號了）。在第一行中的  就是一個噪音向量，我們考慮 zero-mean uncorrelated noise process，即均值為0的不自相關的噪音隨機過程。顯然高斯白噪聲(AGWN)是滿足條件的。  是噪聲影響動態的耦合矩陣。而  是加入在測量中的測量噪聲，我們同樣假設是zero-mean uncorreleated的，同樣也可以是AGWN。

而實際中我們經常會碰到離散的(discrete)線性隨機系統，根據上面的結論，我們可以寫出：

其中  為transition matrix，

 ， .

為什麼我們在模型中總假設噪音都是Gussian white noise呢？實際中的噪聲常常都是有色噪聲，即power spectral density並不總是為常數的隨機噪聲。上一篇中我提到了，當白噪聲通過一個線性系統時，其輸出是有色噪聲。因此，已知有色噪聲的PSD或者autocorrelation function，假定它是高斯白噪聲通過某個線性系統後產生的輸出，是否可以反向推出這個線性系統的模型呢？答案是肯定的。通過一個特定的線性系統，稱之為"shaping filter"[1,p131],使得白噪聲變成了需要的有色噪聲，我們就可以把這部分動態的模型包含到原有系統的建模中。這使得我們可以把原來系統動態中的有色噪聲，替換成白噪聲和shaping filter的模型。

我們舉一個時域連續系統的例子來說明。

當高斯白噪聲通過一個一階穩定線性系統時，其輸出為一個自相關函數為指數函數的有色噪聲。現在我們假設有色噪聲的autocorrelation為：

其對應的PSD，滿足寬平穩過程WSS的條件下，就是它的Fourier Transformation：

我們把這個有色噪聲當做  ，把高斯白噪聲當做  ，服從標準正態分布，即方差  為1，均值  為0。那麼  的自功率譜密度函數  ，其自相關函數為  。根據公式(11)：

其中  代表了convolution運算。這個結論正是公式(14)的來源。我們把(22)兩端進行Fourier Transform，得到了線性系統輸出的自譜和輸入的自譜關係：

那麼現在代入  ，我們得到了：

取穩定的頻率特性，便有了這個我們需要的shaping filter，即一個想像中假定存在的系統：

這就是把標準正態分布的高斯白噪聲轉化為autocorrelation為(20)，或者自譜為(21)的特定filter，我們可以把它寫成時域微分方程：

我們說shaping filter (26)的輸出  正是有色噪聲  ，而輸入則是高斯白噪聲  。注意這裡我們替換了符號，採用了熟悉的  作為狀態，  作為噪聲。有了(26)之後，我們便可以對有色噪聲進行替換，比如狀態空間方程：

其中  是有色噪聲，其autocorrelation滿足(20)。那麼我們 incorporate (26) into (27)，得到：

把(27)重新寫成成新的狀態空間表達式就不是什麼難事了，所謂的augmented state space。我們這種將有色噪聲輸入轉化為等效白噪聲輸入的過程也稱之為白化處理。離散系統的做法並不一樣，可以參考[1,p129, Example 4.4]，不過背後的想法是一樣的。

Classic estimation theory 經典估計理論

估計是用過去到現在的數據來計算現在的某個量的大小。估計理論主要還是屬於概率與統計學領域。我們知道統計學有兩大學派：頻率學派(frequentist statistics)和貝葉斯學派(bayesian statistics)。估計理論也可以根據這兩個學派不同的觀點分為經典估計理論和貝葉斯估計理論。

經典的估計理論我們在以前上統計學課程的時候應該已經接觸過了，主要在參數估計那一塊，包括點估計和區間估計兩大類，其中有比如最大似然估計算法，一致最小方差無偏估計等。我們抽裡面的重點進行回顧和講解，目的是為了和後面的貝葉斯估計有一個對比。Kalman filter的推導用到的估計知識，主要是後者。對比兩種估計理論的不同點，有利於加深記憶。

統計推斷(statistical inference)的任務是要通過研究樣本數據來分析總體的特徵。總體是由分布描述的，現實中很多時候容易知道分布的類型，但卻不知道其參數。那麼如何從樣本數據中得到分布的參數是關鍵問題，即參數估計(parameter estiamtion)。

值得注意的，參數估計和狀態估計問題還不能完全等同起來，這點隨著知識不斷擴展慢慢體會。

2.1 Unbiasedness and consistency 無偏性與相合性

我們下面都以  來作為未知參數。頻率學派認為，未知參數  的真值是一個確定的值(deterministic)。根據總體中的抽樣樣本  ，構造一個統計量(statistic)：

稱之為  的估計量(estimator)，或叫  的點估計(point estimator)，簡稱估計——點估計的名字是因為最後估計量的取值是一個數，與區間估計相對。統計量  是一個樣本的函數(統計量的定義)，也是一個隨機變量(因樣本具有隨機性)。

註：統計量  本身是一個函數，但很多地方符號本身也表示了其取值，類似於  這麼個意思。注意英語中的用詞，estimator表示的是估計量，即統計量，其為樣本的函數。estimate是estimator在某個樣本處的取值，是根據取得的樣本給出的估計值。

 的估計可以構造無數種，估計之間孰優孰劣，需要一些標準來評判。我們說一個估計是unbiased(無偏的)，是指這個估計的期望值與對所有可能的被估計量  的真實值一樣：

 是  可能的取值空間，稱參數空間。無偏性(unbiasedness)告訴了我們這個估計  是在真值  周圍波動的，但如果我們求平均(期望)之後，最後等於真值。除非是  的線性函數，一般  的函數  是不會繼承其無偏性的，比如修正樣本方差  是無偏的，但是開平方根的 卻是有偏的(biased)。 在其他評判標準滿足時，存在無偏估計時，應當儘量使用無偏估計，這意味著估計從原理上就不存在偏差。

有些估計，如樣本方差  本身是有偏的，但隨著樣本數增加，它可以滿足漸進無偏性(asymptotic unbiasedness)，在大樣本時是可以近似認為無偏的。一個擁有較大bias的估計認為是有缺陷的，一般是不會採用的。不是所有參數都存在無偏估計，若  存在unbiased estimator，則稱其是可估的(estimable)。

隨著樣本量  的增加，我們希望估計量不斷逼近參數真值，這種性質稱為consistency(相合性)。把  看成是隨著n變化形成的隨機變量序列，則  依概率收斂到  。判斷相合性一個常用的充分條件是：若  則  是  是consistent estimator. 相合性是統計量的一種大樣本性質。無偏和相合是兩個不同的概念。相合性被認為是很基本的要求，如果某些場合只存在有偏估計，那麼至少也要採用一個相合的估計。

2.2 method of moments 矩法

常用和常見的一種估計方法是矩法(method of moments)，其思想使用樣本矩(sample moment)去替換總體矩(population moment)，然後通過樣本矩反推總體分布的參數值。其背後的思想基于格裡汶科定理——經驗分布函數隨著樣本增加接近總體分布。

根據這個方法，即便有些時候分布形式是未知的，我們依舊可以估計一些參數：

比如用樣本均值  和樣本方差  估計總體均值  (一階原點矩  )與總體方差  (二階中心矩  ) 。

在分布和概率函數形式已知時，我們可以寫出該隨機變量  的  階矩表達式(以原點矩為例)：  ，它們應當是概率函數中未知參數  的函數。

用樣本原點矩  代替  就可以列出  個方程來求解  。

此法簡單實用，不過不一定得到無偏估計量，但一般滿足相合性。

2.3 Maximum Likelihood Estimation(MLE) 最大似然估計

MLE(最大似然估計，或極大似然估計)是非常傳統的頻率學派估計方法，當然學過了Bayesian estimation之後，我們又把它當作了一種maximum a posteriori (MAP) estimation的一種特例。MLE背後的思想非常簡單易懂的：認為已經發生的事情是最有可能發生的。

相信各位已經看厭了投硬幣的例子，那我編個讓各位感同身受的例子：你通常喜歡在打遊戲時反鎖家門。上個月由於天天把家門反鎖打遊戲，被你爸發現罵的狗血淋頭。這個月開始你洗心革面，決定不打遊戲了，天天在家裡看知乎學習。結果今天不小心你又把家裡門給反鎖了，傍晚爸爸回家看到了門被反鎖，問他是相信你當天在打遊戲還是沒打遊戲？

這個問題的題設有很多隱藏的信息需要解讀。「當天家裡門反鎖狀態」是一個隨機變量  ，滿足伯努利分布  ， 是家裡門反鎖的概率。根據你打遊戲喜歡鎖家門的習慣，我們認為  的取值受到「當天你是否在打遊戲」這個參數  的影響，  參數空間只有兩個取值「打」或者「沒打」。

很明顯，當天打沒打遊戲成為了當天家門反鎖狀態  分布列的一個至關重要的參數 。假設真實的概率情況：你當天打遊戲99%概率門是反鎖的，1%門不是反鎖的。不打遊戲則1%概率是反鎖的，99%不反鎖。那麼真實的分布列只有兩種可能：

 當天打遊戲 ——  ="當天反鎖」  ； ="當天沒反鎖」 

 當天沒打遊戲 ——  ="當天反鎖」  ； ="當天沒反鎖」 

老爸當天晚上來到門前，這個參數  「打」還是「沒打」已經是確定的事情了，只不過在他看來是未知的（其實你認真在看知乎）。今天門是否反鎖的概率為  。他從上個月整整30天的觀察中猜測了  與  之間的關係：打遊戲，家裡門一定反鎖。這意味著他認為  "打"，  "當天反鎖" ，據此此時可知，打遊戲時門反鎖概率  。而如果  「不打」， 他想了想，以前也沒有遇到過其他門反鎖的原因，於是認為  "當天反鎖"，此時  。他伸手去開門後，發現今天門被反鎖了！

綜合以上，這兩種情況，哪種可能性大一目了然！他微微一笑，晚餐的氛圍當然不會很融洽——這小子又TM在打遊戲。結果巧了，超小概率事情今兒發生了——咱們真的在認真看知乎（那你鎖個毛線門啊，喂！）

爸爸的推斷用到了最大似然估計的思想。我們假定一個隨機變量  的分布函數包含  個未知參數需要估計，它們組成了向量  ，它們都是確定值。所以這  取值必然也直接影響到  的  個樣本  的聯合概率密度函數，即  。將這一個joint probability density function看作是  的函數，我們把定義為likelihood function(似然函數)  ：

最大似然估計就是取estimator滿足似然函數的取值最大化： 

老爸運用了MLE給了判斷，不過運氣不好小概率事情發生了！我們再次分一下整個過程：當天家門反鎖情況  的真實分布是受到  這個參數的影響。當天你沒有打遊戲，所以  的真值其實是「no game today」。爸爸根據伯努利分布知道今天門被反鎖的概率為  ，馬上檢查了一下門，被反鎖了！

據此寫出了樣本的聯合概率分布函數（這裡只有當天的一個樣本可檢查，所以實際上只是單變量概率分布函數），亦為此時的似然函數  :

我們尋找一個  使得(32)最大化，則它就是我們給出的MLE。這裡有一個很重要的問題，  與  之間的關係是什麼？往往在丟硬幣的例子中，我們把這個  當做了  ，這是我們對確定分布形式得到的一個必然結果，因此我們只要選擇大的  就可以令(32)最大化。這裡「打遊戲」和「當天家門是反鎖的」的關係不是數值關係。

對於既定的問題，我們必須要明確給出或者假定分布函數和參數的關係，才能讓我們選擇參數來最大化似然函數。爸爸在對  關係的猜測是基於過去30天經驗給出的—— 所以他得到的結論是：  =「打」會讓  ，  ="不打"，則反鎖概率  。顯然「今天兒子打遊戲」這個判斷能夠使得「今天門上鎖」這件事的概率達到最大。所以他自信推開門，想要開始教訓你。

很有意思的是，真實的情況，你打遊戲時有99%的可能會鎖門（剩下的1%不知道你在幹什麼鎖門，難道是在看...？），所以真實的關係應該是，這樣得出的概率 =「打」  ， ="不打"，則反鎖概率。無論是那種情況，只要爸爸相信打遊戲時門上鎖的概率更大，那麼他一發現當天門的反鎖了，就會得出一樣的結論：你今天打了遊戲。

從條件概率的角度看，MLE試圖找到最大化條件概率密度  的條件  。

這個例子與常見的丟硬幣用MLE估計正面概率有微小區別。因為成功概率  和參數  的關係是可能改變的，取決於爸爸如何看待，而不像投硬幣時  是確定的。不管是什麼關係，在做出MLE前，它們之間的關係應該是明確的，這樣才有優化似然函數的方向，就像爸爸相信打遊戲那天門反鎖的概率遠比不打時的大。

這其實意味著，這個問題中「當天家門反鎖」這個隨機變量的分布形式不僅僅受到  取值的影響，還受到了  關係的影響。這與一般MLE中我們假定分布形式確定，而參數未知還是有區別的。一旦  關係確定，我們就和投硬幣問題沒什麼差別了。分布形式的確定是至關重要的，這會影響整個估計結果。  關係也算是根據經驗得到的信息，是先驗信息的一種。只有這個確定了，我們才能最終總體  分布的具體形式。

貝葉斯估計中，認為參數  本身也是個隨機變量，即: 那天我打沒打遊戲都不是個確定的事情，應當根據經驗給出一個  ，這是參數的先驗分布函數。老爸可以根據過去經驗猜測，預先猜測了今天我大概率在打遊戲。我們後面會再提到。

最後總結一下，MLE邏輯是這樣的：

有一個或多個固定的參數  影響隨機變量 樣本的聯合概率分布函數(因而影響了聯合概率密度函數)。也就是說，  的在參數空間中的不同取值會影響某一個樣本 出現的概率。

既然如此，我們寫出樣本的聯合概率分布函數或者PDF，那它就是樣本和參數  共同函數。

對於我們抽樣得到的某個樣本，我們認為既然它發生了，它聯合概率就應該是最大的。把這個樣本代入到聯合概率分布函數或者PDF後，它就只是參數  的函數了，定義它為似然函數  。

既然聯合概率要最大，找到  使得  最大，即得MLE。

對數似然函數  用的比較多，因為它本身也是單調函數，且解析解求取比較方便。實際中似然函數多是用數值求極值的多，會用到牛頓法。根據似然函數對參數的一階導數為0，MLE的解並不一定存在(解跑出參數空間之外)，也可能存在多解的問題。

理論上MLE並不總是能給出無偏估計，不過仍可能是漸進無偏的。再加上MLE實際上需要你對統計模型有一個很好的建立，實際中的無偏也沒法實現。MLE通常被認為是consistent的，雖然因為模型誤差實際也不會得到，而且某些概率模型還真會讓MLE不是consistent的[c]。但是consistency依舊是能保證就要保證的，即便是理論上的。關於MLE我們就說到這裡，我們在後面的貝葉斯估計中可能還會call back。

關於樣本的聯合密度概率函數，樣本抽樣不一定要滿足i.i.d，即獨立同分布。當然很多時候滿足的話，可以讓我們直接將所有樣本的概率密度函數直接相乘，是個不錯的情況。我看很多人寫文章都假定了i.i.d，尤其是講機器學習的。故而提示一下，其實這不是一般的情況。

2.4 Mean Square Error and Minimum Variance Unbiased Estimator 均方差、最小方差無偏估計

估計量的unbiasedness和consistency並不意味著滿足了我們對估計的全部要求。我們看來自[b]中的下圖：

無偏 vs 有偏 來自[b]

上面左圖中的點可以代表一個參數，來自不同樣本的無偏估計量分布。我們可以看到，取不同樣本時，無偏估計的波動可以非常大，雖然我們可以大概看出其在紅點(真值)附近比較均勻的波動的，即期望仍然是真值。右圖是可以代表有偏估計的估計量分布。我們明顯能看到這些估計量都十分靠近彼此，意味著估計量隨著樣本採樣的變化波動很小，雖然我們一眼就看出來了其與真值之間總是存在某個方向的偏差。

我想有學過機器學習基礎課的朋友看到這裡的圖，應該一下就想起來了，這很像是討論所謂的 Bias–variance tradeoff 問題，可能會用到accuracy來描述bias大小，用precision來描述variance大小。沒錯，統計學和機器學習本來就有很多內容是一樣的，這裡就是個很典型的例子。

評價點估計，除了使用無偏性之外，我們還引入了均方誤差(Mean square error, MSE)來評價參數估計的優劣。我們發現有偏還是無偏，考慮的是估計誤差  的一階矩(期望)，而MSE考慮的就是  的二階原點矩：

注意！此處我們假定了  的真值是固定值，所以這種將MSE分為Variance和Bias平方的思想是frequentist的思想。這裡是一個很容易和Bayesian estimation搞混的地方——它們認為  本身也是個隨機變量，於是我們就不能令  成為一個const了，於是也不會去討論將MSE分成Var與Bias的平方兩部分。事實上，  都是隨機變量時，MSE的計算就要考慮它們之間的聯合概率密度函數。

MSE的意義是明確的：估計參數  與真值  之間的距離的平方。當認為  是固定值時，它是由方差與偏差的平方兩項組成的。同理，我們可以定義均方根誤差，直接開平方就可以了。當估計量本身是無偏的時候，我們就有 ，即估計量的均方差在估計無偏時就等於估計量的方差。

MSE是一個非常一般的估計評判標準，我們自然希望MSE能夠越小越好。比較有偏的估計，方差和偏差的大小都應該考慮進去。這也啟示我們：無偏估計雖然只有方差一項，但仍然可能使得MSE大於有偏估計。[see a page 324 ,example 6.4.1, b example 2.13] 此時有偏估計將比無偏估計在MSE評判下更優，或稱之為更有效。

看到MSE，很多在學或者學過機器學習基礎知識的同學可能會好奇，這裡的MSE定義和我最後看到的形式有點不一樣啊。比如在regression問題中檢驗data fitting好壞時常定義的MSE是：

其中  是訓練集或者測試集中  的target，  是根據predictor預測的值。通過計算預測值與target的MSE，我們評判predictor的優劣。所以實際上那裡的MSE是求已知樣本值與predictor預測值之間誤差的二階原點矩，或者說誤差的平方的樣本均值。理論分析一下，已知訓練集或者測試集中的target值  是真值  外加uncorrelated噪聲  ，那麼理論上MSE應該寫成：

對比(54)與(52)，只是因為我們使用的不是真值，而是觀測的值，所以產生了一項  。(53) 採用誤差平方的樣本均值來近似計算MSE，是一種矩估計，是相合的。所以(53)這樣的MSE定義是針對predictor的。給定訓練集中的target的真值還是確定的，因此仍然有(52)那樣分解得到Var和Bias的平方，這裡只不過多考慮了噪聲。

講到這裡你會發現這些都是統計學裡的參數估計內容，而我們最終關心的還是偏狀態估計。有三個問題我們要注意：

 雖然是未知的，但其取值是有範圍的，我們說  的所有取值構成了參數空間  。我們在統計學中定義估計量的特性時，通常會針對參數空間中任意的  取值。比如無偏性，我們希望無論  是如何取值的，一旦通過某個分布——其由任意參數  確定——產生的樣本去估計  ，產生的  都是無偏的。比如   時產生一組樣本  ，我們對此樣本數據，使用  計算的估計量是無偏的。而  確定的一組樣本  產生的估計量  也應該是無偏的。

既然有了MSE這個標準，我們自然想知道，從眾多可能的estimator中是否能找到一個這樣的  ？對 都，都有  ，其中  是任意一個  的估計。換句話說 的對所有  可能取值的估計值，總有比其他任何estimator的MSE，即它在參數空間  中所有取值上都最小化了MSE。在頻率學派的觀點下，當  是確定值時，對所有可能形式的estimator，這樣能夠全局最小化MSE的estimator通常是不存在的！

這是一個很容易與稍後Bayesian MMSE混淆的地方，在那裡我們明明定義了一個Minimum MSE的estimator！在經典估計理論中，因為我們並不可能知道真實的參數值  ，MSE中的Bias平方項的值將取決於真實的參數值。假設存在一個能夠全局最小化MSE的estimator  ，現在我們取一個特殊的estimator  使得其恆等於參數空間中的某一個參數  ，此時對  ，  。由於  是全局最小化MSE的estimator，此時必然有  。仔細一想，你就發現這是不可能的。這種最小均方差estimator在參數確定的情況下，只存在於規定形式的一類估計器中。

當我們只考慮unbiased estimator時，MSE中的Bias項就消失了，MSE就退化為估計量的方差。很多時候我們確實可以找到一個所謂Uniformly Minimum Variance Unbiased Estimator(UMVUE, or MVUE, 一致最小方差無偏估計，或最小方差無偏估計)——它相比於所有其他的unbiased estimators，對所有  的取值都有最小化的  。UMVUE並不總是存在的，有時候一個參數的unbiased estimator壓根不存在；或者即便是存在的，也有可能沒有一個estimator能夠對所有參數實現最小Variance，如下圖：

EECS 564: Estimation, Filtering, and Detection. Umich

Cramer-Rao Lower Bound(CRLB,C-R下界) 指明了任何一個unbiased estimator的方差的下界，可以幫助我們有時確定UMVUE。為了得到CRLB，首先要介紹Fisher Information（費希爾信息量）。在滿足一些Regularity conditions時(這裡先不展開)，Fisher Information 被定義為：

其中  是log-likelihood function(對數似然函數)。此處我們沒有進行任何樣本獨立同分布的假設。為了對(55)產生一點直觀的理解，讓我們考慮一個十分簡單的參數估計問題：

 ，  是一個未知參數待估計，  是測量噪聲，其滿足正態分布  。如果我們只抽樣了一個樣本  ，它應該滿足正態分布  。那麼其log-likelihood function是：

顯然對(56)求導，我們得到：

令(57)為0，我們得到了  的MLE為  。如果我們繼續求導，就可以得到：

根據(58)知道最大似然函數應該是一個開口朝下的曲線，要尋找最大值。考察函數的曲率，由於此時一階導數(57)非常接近於0（事實上求期望可知E(x)=E(A)，即  ），那麼二階導數的絕對值就可以基本等同於曲率。

我們發現在這個例子中，log-likelihood的最大值附近的曲率等於  。這意味著如果方差很大，即測量時噪聲功率很大，那麼log-likelihood曲線在最大值附近的曲率就小，導致整個函數看起來「又矮又平」；反之，噪聲功率小，樣本方差小，曲線在最大值附近看起來會非常陡峭。

一個平緩的似然函數並有利於我們縮小參數真值的範圍（自變量是參數  ），這也是合乎常理的，因為噪聲功率大了，自然就很難分辨真值了。所以我們看到  的二階偏導的越大，表示目前這些樣本能夠更好地反映真值，即已經觀察到的樣本包含了更多的  的信息。可以證明最大似然估計MLE在特定的regularity condition下滿足漸進正態性(asymptotic normality), 即隨著樣本量增加分布函數向正態分布  漸進，此時還滿足相合性以及漸進無偏性。此時Fisher Information的倒數構成了它的漸進方差。

讓我們回到log-likelihood function  ，對待估計參數  的一階偏導被稱為score function，可以記作  。可以證明滿足概率函數滿足某些運算條件時， 。score的方差被定義為Fisher Information：

可以看到Fisher Information也可以表示為似然函數二階偏導期望的相反數，可以視作似然函數curvature(曲率)。

注意最後計算  時這裡的似然函數中的  取的觀察到的樣本。

Cramer-Rao Lower Bound指出當  滿足時，任何  的無偏估計  的方差必須滿足：

並且以上不等式的等號成立，若且唯若  . 滿足等號的unbiased estimator就是UMVUE，最小方差，即CRLB為  。對於i.i.d的樣本，可以把總體密度函數寫成各樣本密度函數乘積，則可以定義單個observation的Fisher Information  , CRLB此時為  。

注意Fisher Information有些書上定義時採用population密度，有些假設了i.i.d採用了單個observation的密度。在茆詩松的書[a]上p329，他明明定義的population的FI，但是證明定理6.4.3時又採用了observation的FI，應該是個筆誤。

我們可以看到有一些無偏估計確實達到了CRLB，比如用樣本均值  估計正態分布  的  。根據(57)有

  .

所以

  

而我們知道sample mean的方差

 

達到了CRLB，因此此時  是  UMVUE. 可以看出了我們採用樣本均值的合理性。

此外，還有vector form的CRLB定理，適用於多於一個參數的情況，這裡就不介紹了。既然已經知道了概念和用途，我們需要時自己去查閱資料就好了。

2.5 Efficiency 有效性

定義CRLB之後，我們可以更好解釋efficiency這個概念。一個unbiased estimator的efficiency被定義為：

所以實際上就是理論最小方差比上實際方差，CRLB和efficiency給了一個判斷unbiased estimator之間優劣的標準。當達到了CRLB，我們就稱這個estimator是efficent的。更多關於efficiency參考wiki裡的一些內容做一個簡單了解[e]. 這些涉及細節定義上的東西，儘量不去多展開了，有需要自己做擴展閱讀會比較好。

2.6 Best Linear Unbiased Estimator 最佳線性無偏估計(BLUE)

CRLB可能可以驗證一些UMVUE(當它們達到CRLB時)，但是更多的時候這個下界是達不到的，一個UMVUE並不容易找到。如果我們能夠把unbiased estimator的形式限制在樣本的線性函數，那麼我們便可以在linear unbiased estimator中找到UMVUE了。但是這種Variance的比較僅僅針對linear unbiased estimators之間而言，並不意味著它們達到了CRLB。所以BLUE的結果其實說是best，其實算是suboptimal的。

Fig from lecture notes [i]: Relations among Unbiased Estimators

回顧MLE和一些矩估計中我們通常需要對概率密度  或者分布形式  有一個假設或者了解，概率模型的確定對結果會有很大影響。在上面計算Fisher information時也用到了這些信息。而對於BLUE，我們只需要知道  的first moment(一階矩)和second moment(二階矩)信息，連具體是什麼分布都不需要知道，規避了完全掌握密度函數形式的要求，非常實用。

前面一直是從統計學的角度看待估計問題。我們嘗試從信號處理的角度去求解BLUE。假設我們觀測到了的一串受到噪音汙染的信號序列，記為向量  。真實的無汙染信號序列為，記為向量 。 的值受到  個參數組成的向量  的影響。  受到了噪聲序列  的影響，其均值為  ，協方差矩陣  。  的維度就是觀測的數量  ，待估計參數的數量記為 (  往往小於 )。對於第  個觀測值，它們之間一般有如下量測關係：

其中  是給定的一般的非線性函數。我們的任務是通過  構造  的估計量。

為了得到  的BLUE，我們列出所需滿足的條件限制：

針對第一條：我們得到  的形式，其中  矩陣被定為  ，  是每個參數  的線性項係數組成的列向量。故對每一個參數都有：

針對第二條：對  求期望

為了滿足(64)，我們必須構造一個矩陣滿足  。這樣我們就有

這意味著  與  之間必須是仿射關係：

我們可以看到，x的隨機性在這裡完全是由噪聲  引起的。而  的值由於參數真值是固定的不存在隨機性。我們看到噪聲的協方差矩陣  就是測量信號的協方差矩陣

 

針對第三條：我們計算  的協方差矩陣：

協方差矩陣的對角線元素之和就是我們需要的方差之和。

則我們要最小化  同時使得  。

令  ，則  對  求偏導得到：

所以BLUE已經找到：

 的協方差矩陣此時變成了估計誤差的協方差矩陣：

注意  是協方差矩陣，故而是對稱矩陣。可逆與轉置是可以交換的。

整個過程求解BLUE的過程，我們用到了幾個信息：

(69)的解是在我們認為噪聲均值為0，協方差矩陣  已知時求出來的。如果我們假定噪聲是符合正態分布的呢？假定此時我們知道噪聲是高斯白噪聲，其聯合分布符合  。由於白噪聲互不相關，於是滿足i.i.d，其互協方差均為0。於是噪聲  的協方差矩陣只有對角線上擁有非零元素，它們都等於噪聲的方差  ,即：

將(71) 代入(70)，則得到此時的BLUE擁有協方差矩陣：

我們知道CRLB告訴了我們UMVUE的最小方差，而在多變量時我們有Fisher Information Matrix(FIM)來告訴我們一個估計的協方差的下界：

此時的。因而當噪聲是高斯白噪聲時，此時的BLUE(線性無偏最佳估計)：

並且它達到了CRLB，同時也是UMVUE！即沒有其他任何無偏估計能比它產生更小的方差 Gaussian-Markov定理也表明了此時(74)的成立，但定理不保證UMVUE。所以UMVUE只是正態分布下的特殊結果。

細心的你可能已經發現了，只要噪聲滿足方差相同，不相關，均值為0，都會產生(74)的結果。沒錯！這就是Gaussian-Markov定理的內容。但是如果不滿足正態分布的話，(74)就只是BLUE，並達不到CRLB。而不是UMVUE了！你可以仔細學習一下G-M定理的內容。

有色噪聲的協方差矩陣還有對角線之外的元素。那麼問題又來了，有色噪聲不是白噪聲怎麼辦？請複習前文：白化。對於高斯有色噪聲，如果我們知道其協方差矩陣  ，那麼(69)不僅僅是  的BLUE，同時也是UMVUE！

2.7 Ordinary Least Squares 最小二乘算法

Ordinary least squares(OLS, or least squares, 最小二乘法)是最為經典的最優估計算法。實際上，經過前面一連串的分析，我們得到了一個有意思的結論：

Gauss–Markov theorem（高斯-馬爾科夫定理）：在線性回歸模型中，如果誤差滿足零均值、同方差且互不相關，則回歸係數的最佳線性無偏估計(BLUE, Best Linear unbiased estimator)就是普通最小二乘法估計。

什麼是線性回歸模型呢？關鍵是因變量和參數為線性關係。

根據定理，我們可以立刻寫出OLS的表達式：

正是(74)。但是與之前任何估計方法不同的是，它沒有用到概率密度函數的任何信息，因此適用的範圍就大大擴大了。我們常見的線性回歸問題中，OLS是最基礎且常用的。

注意到(74)是BLUE的條件，那就是：

我們在上面已經說過了，協方差矩陣此時就只有對角線上的數值，即各分量的方差。各方差都是相同的。定理中甚至都不同滿足i.i.d條件。現在假設有高斯分布的噪聲，我們知道這種情況下各條件均是滿足的，於是就有了最小二乘法此時是BLUE的結論（雖然在高斯分布下還是UMVUE）。換句話說，如果此時噪聲各分量方差不同，且還互相相關（比如有色噪聲），那麼OLS必然不是BLUE。

那BLUE是哪個呢？你想想？答案前面翻一翻。下面也提到了。

好了，下面就舉一個OLS參數估計的例子。

OLS算法是不需要概率模型的，也不必一定要滿足Gauss–Markov theorem的要求成為BLUE。它需要一個線性測量模型和足夠且合適的樣本觀察。典型的應用是可以為一個方程數比未知數多的線性模型，稱為overdetermined system，估計出近似最優解——這組解可以使得近似解的殘差(residual)的平方和是最小的。舉一個例子來說明。考慮一個參數估計問題：

這是一個典型的單自由度振動運動方程，  是質量(mass)， 是阻尼係數(damping coefficient)，  是剛度係數(stiffness)，  為作用在物體上的外力。假設我們手上有加速度、速度和位置的數據，即  ，以及外力的信息  ，試求未知常係數  。

從解方程的角度來看，這個問題非常容易，因為有3個unknowns，只要3組不同的數據建立3個方程就可以解出來了：

當然手上這些收集的數據其實都是有噪音的，所以任意選3組數據解出來的  值可都不一樣。現在手上的數據嚴重冗餘，比方說有  吧，可以任意選擇其中3組來構成方程組(76)，每一種選擇都有對應的一種估計  ，那麼問題來了： 哪一種選擇是比較好呢？有沒有辦法把所有的數據都儘量用起來，不造成信息浪費呢？有的，比如把所有可能的解都求出來，求個平均。是不是很麻煩？Least Squares就很好地解決了這個問題。考慮把所有的數據代入(75)得到矩陣方程：

記(31)為  . 我們定義一個參數估計，也就是(77)的近似解  ，並定義其代入方程後與對應測量值  的殘差，或者說誤差：

其中  是  的第  行。求令所有  的平方和最小時對應的解，即令：

其中  是  的第  行，第  列的元素。(79)對  每一個元素  求偏導，並且使偏導為0，得到：

因為  , 那麼就有：

 

等式的左邊是  第  列的轉置與  的乘積，右邊是  第  列的轉置與  的乘積。當  時，我們就有[1,p8]：

方程(82)稱之為normal equation. 式(82)的解，當  是nonsingular時，為：

這時  是least squares意義下的最優參數估計。可見我們只要把加速度信號，速度信號和位置信號排列成矩陣  ，外力測量值排列成向量  ，我們就可以得到一個唯一的解，當然前提是  是可逆的。  的可逆性取決於  是否是列滿秩的，實際上  也就是當  列滿秩時的Moore-Penrose pseudoinverse(廣義逆矩陣，或稱之為偽逆矩陣)。實際遇到的矩陣  基本都是列滿秩的(要是剛剛不列滿秩，運氣得有多好)

註：  是可逆的，那麼就有  ，我們就得到了  的左逆為：  。這種超定方程  本來可以通過求逆  的方法來直接得到。對於非方的  ,我們可以通過求廣義逆來得到  ，也就是(83)。所以在Matlab中，你可以用 y=inv(H'*H)*H^T*f 也可以直接用 y=pinv(H)*f, 結果是一樣的。Matlab是通過SVD來求解的，這裡的結論採用了滿秩分解。關於線代不深究了，用到再說。

最小二乘法作為一種線性估計方法，其適用於估計值是觀測值的線性函數。我們把之前的觀測關係放到這裡就能看到它們是一樣的：

在這個參數估計問題中，在模型中加不加上噪聲並不重要，因為我也沒有噪聲的任何信息可以利用。我們給出的根據誤差平方和最小化的估計結果是：

注意到這裡的優化標準和BLUE並不一樣，誤差平方和並不是MSE，因為我們沒有求取期望。其實  中各誤差分量的權重是一樣，下面我們細講。

2.8 Weighted least squares 加權最小二乘算法

如果我們觀察OLS的指標：

(怎麼定義  的符號並不影響，這裡與我上面的例子中正負不相同)在最小化  的過程中，每一項誤差的平方  在最小化過程中是有相同地位的，即它們的權重是一樣的。這樣的設定暗示著我們把每一個觀測值的(決定估計量時的)重要程度都是一樣的。然而不是所有時候觀測值都是靠譜的；有多傳感器融合時，不同觀測值的統計特性有很大差別，比如方差大小不同。這一塊我們在下一篇中寫誤差分析總結時會再次提到。總之，給不同的誤差賦予不同的權重可以把這種差異主動考慮進估計結果之中，即產生了所謂的Weighted least squares(WLS, 加權最小二乘算法)。考慮WLS的指標為：

其中  是加權矩陣，  是  對角線上的元素。  矩陣的形式應該很容易想像吧？就是只有對角線元素的方陣。我們直接給出WLS的估計公式：

推導過程和LS是一樣的，沒有什麼難度。

其實(88)我們已經很熟悉了，沒錯！當  時，它和BLUE的結果是一樣的！這意味著權重矩陣並不是隨意的，當權重矩陣和噪音的協方差矩陣互逆時，能使其成為線性無偏最優估計BLUE。

2.9 Recursive least squares 迭代最小二乘算法

傳統的LS和WLS說起來算是一種「離線」估計，或者說批處理——每次都會把之前所有的觀測值一起使用並直接給出估計值。計算機執行這種算法時，需要不斷開闢存儲空間來滿足不斷增加的觀測量，以及矩陣計算會越加費力。在線估計的想法是，根據當前的觀測值計算出現有的估計量，再根據新的觀測值去更新原有的觀測值，而不再需要重新使用歷史中的觀測量。這樣大大降低了計算機的計算負擔和存儲負擔，有利於算法的應用。LS和WLS都可以根據這種思想改為recursive least squares，即RLS算法。

考慮使用  個觀測值進行WLS估計，我們有觀測方程  ：

其中  是矩陣  的第  行，  是行向量。考慮一個  維權重方陣  ，根據(42)有WLS估計量:

其中  。考慮在此基礎上，新添加一個觀測方程，使用  個觀測值進行WLS估計，同樣有：

其中  。因為  , 

 ,  ，我們代入(45)中得到：

又因為其中:

即：

所以我們其實就是把兩個公式進行了某些項的拆分並建立前後的聯繫。把(44)和(47）代入(46)，我們有：

最終有Recurive WLS 迭代公式為：

其中  為迭代的增益向量，會隨著迭代的進行而發生改變最終收斂。相信做過優化的朋友應該對(50)的形式不會感到陌生，非常經典的形式。其中  正是採用前一次估計量  時產生的residual(殘差)。

我們考慮一個十分簡單的問題。如果有  個離散數據，我們嘗試用一個常數  來估計這堆數據，也就是說我們假設這對數據是圍繞一個常值  去波動的。那麼每個量測方程就滿足

很自然的，我們就知道這裡的  。LS算法給出的估計結果就是：

就是所有觀測量的平均數！如果我們考慮用RLS來計算的話，

 

根據(50）有最終結果：

可見這裡的增益  會隨著迭代進行不斷衰減，最後收斂的結果和批處理的LS結果是一致的。

2.10 About other method 其他方法

還有不少算法其實在機器學習中經常出現，我可以對目標函數加以限時，得到LASSO或者ridge regression等等都是非常實用的參數估計算法。感興趣的同學可以多搜集搜集，並且學著編編程序來實現他們。這裡和KF關係有點遠了，就不繼續往外擴展了。

本篇開始回顧了線性系統和信號處理的一些知識，重點是自功率譜函數和互功率譜函數，它們和FRF之間的關係，以及有色噪聲的白化。

第二部分我們回顧了統計學中經典估計理論和眾多實用的結果。我們從無偏和相合的基本概念出發，重點介紹了MLE，MVUE，BLUE以及最小二乘法和它的迭代算法。這一部分主要針對的還是參數估計，記住其重要的假設就是參數是確定的。這也是所說的為什麼沒有通常沒有minimum MSE的estimator存在。在後續的貝葉斯估計中我們看到，參數也可以作為隨機變量，我們可以給定或者猜想先驗概率  。貝葉斯估計又將帶領我們進入另外一個全新的世界，以及我們要去那裡迎接Wiener filter——Kalman當初的文章也是為了解決Wiener Problem。

了解來龍去脈，方可遊刃有餘。

Main Reference

[1] Grewal, M. S., & Andrews, A. P. (2015).Kalman filtering: Theory and Practice with MATLAB Fourth Edition. John Wiley & Sons.

[2] Anderson, B. D., & Moore, J. B. (1979).Optimal filtering. Prentice-Hall, Englewood Cliffs, New Jersey.

[3] Gelb, A. (1974).Applied optimal estimation. MIT press.

[4] Musoff, H., & Zarchan, P. (2009).Fundamentals of Kalman filtering: a practical approach. American Institute of Aeronautics and Astronautics.

[5] 王可東(2019).Kalman濾波基礎及MATLAB仿真，北京航空航天大學出版社.

[6] Kalman, R. E. (March 1, 1960). "A New Approach to Linear Filtering and Prediction Problems." 

Extended Reference for this article

[a] 茆詩松等，概率論與數理統計教程 第二版, 高等教育出版社

[b] 何志堅，數理統計講義

[c]  .

如何領取參考資料

1) 關注睿慕課公眾號

2) 回復「卡爾曼」獲取下載連結