估計理論

估計理論是統計學和信號處理中的一個分支，主要是通過測量或經驗數據來估計機率分布參數的數值。這些參數描述了實質情況或實際對象，它們能夠回答估計函數提出的問題。

例如，估計投票人總體中，給特定候選人投票的人的比例。這個比例是一個不可觀測的參數，因為投票人總體很大；估計值建立在投票者的一個小的隨機採樣上。

又如，雷達的目的是物體（飛機、船等）的定位。這種定位是通過分析收到的回聲（回波）來實現的，定位提出的問題是「飛機在哪裡？」為了回答這個問題，必須估計飛機到雷達之間的距離。如果雷達的絕對位置是已知的，那麼飛機的絕對位置也是可以確定的。

在估計理論中，通常假定信息隱藏在包含雜訊的信號中。噪聲增加了不確定性，如果沒有不確定性，那麼也就沒有必要估計了。

## 使用估計理論的領域

有非常多的領域使用參數估計理論。這些領域包括（當然不局限於以下列出的領域）:

- 信號處理
- X射線斷層成像
- 腦電圖
- 心電圖
- 核磁共振
- 醫學超聲波掃描術
- 雷達、聲納、地震學——物件的定位
- 噪聲方差
- 參數化（例如周期圖和相關圖譜）分析
- 非參數化（例如MUSIC、Root-MUSIC和ESPRIT）譜分析
- 維納濾波
- 粒子濾波器
- 臨床試驗
- 民意調查
- 質量控制
- 通訊
- 信道參數
- DC增益（請看下邊的例子）
- 控制理論
- 卡爾曼濾波
- 隨時間改變的執行器（英文：Actuator）
- 網路入侵偵查系統

測量參數包含噪聲或者其他不確定性。通過統計機率，可以求得最優化的解，用來從數據中提取儘可能多的資訊。

## 估計過程

估計理論的全部目的都是獲取一個估計函數，最好是一個可以實現的估計函數。估計函數輸入測量數據，輸出相應參數的估計。

我們通常希望估計函數能最優，一個最優的估計意味著所有的資訊都被提取出來了；如果還有資訊沒有提取出來，那就意味著它不是最優的。

一般來說，求估計函數需要三步：

- 為了實現一個預測單個或者多個參數的所期望的估計器，首先需要確定系統的模型。這個模型需要將需要建模的過程以及不確定性和和噪聲融合到一起，這個模型將描述參數應用領域的物理場景。
- 在確定模型之後，需要確定估計器的限制條件。這些限制條件可以通過如Cramér-Rao不等式這樣的方法找到。
- 下一步，需要開發一個估計器或者應用一個已知的對於模型有效的估計器。這個估計器需要根據限制條件進行測試以確定它是否是最優估計器，如果是的話，它就是最好的估計器。
- 最後，在估計器上運行試驗或者仿真以測試性能。

當實現一個估計器之後，實際的數據有可能證明推導出估計器的模型是不正確的，這樣的話就需要重複上面的過程重新尋找估計器。不能實現的估計器需要拋棄，然後開始一個新的過程。總的來說，估計器根據實際測量的數據預測物理模型的參數。

## 基礎

對於給定模型，估計器需要若干統計 "成分"才能實現。第一，統計樣本從長度為 N 的隨機向量（Random Variable，RV）中採樣獲得，觀測值構成向量：

$$\mathbf{x} = \begin{bmatrix} x[0] \\ x[1] \\ \vdots \\ x[N-1] \end{bmatrix}.$$

第二，有 M 個參數：

$$\boldsymbol{\theta} = \begin{bmatrix} \theta_1 \\ \theta_2 \\ \vdots \\ \theta_M \end{bmatrix},$$

它們的值需要被估計。第三，用於生成連續數據的機率密度函數（Probability density function，PDF）或離散數據的機率質量函數（Probability mass function，PMF）以參數值為條件（這些機率函數潛在存在），即條件概率為：

$$p(\mathbf{x} | \boldsymbol{\theta}).\,$$

參數自身可能也存在機率分布（如貝葉斯統計），此時就需要定義貝葉斯機率：

$$\pi( \boldsymbol{\theta}).\,$$

模型形成後，目標是估計參數，估計的參數通常表示為 $\hat{\boldsymbol{\theta}}$，其中 $\hat{\cdot}$ 表示估計值。

常用的估計器包括最小均方誤差（Minimum mean squared error，MMSE）估計器，它利用了估計參數和參數實際值之間的誤差：

$$\mathbf{e} = \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}$$

作為優化的基礎。該誤差項平方的期望對MMSE估計器來說是最小的。

## 估計函數（估計子）

以下是一些相關的估計函數以及相關的主題

- 最大似然估計（Maximum likelihood estimation，簡稱MLE）
- 貝葉斯估計器（Bayes estimator）
- 矩估計（Method of moments estimators，簡稱MME）
- Cramér-Rao界
- 最小二乘法（Least squares）
- 最小均方差（Minimum mean squared error，簡稱MMSE）
- 最大後驗概率（Maximum a posteriori probability，簡稱MAP）
- 最小方差無偏估計（Minimum variance unbiased estimator，簡稱MVUE）
- 非線性系統識別（Nonlinear system identification）
- 最佳線性非偏估計（BLUE）
- 非偏估計，見偏差 (統計學)。
- 粒子濾波器（Particle filter）
- 馬爾可夫鏈蒙特卡洛（Markov chain Monte Carlo，簡稱MCMC）
- 卡爾曼濾波
- 維納濾波

## 例子：高斯白噪聲中的直流增益

考慮由$N$個獨立採樣點構成的離散信號$x[n]$，它由常數$A$和零均值、方差為$\sigma^2$的加性高斯白噪聲$w[n]$（即$\mathcal{N}(0, \sigma^2)$）構成。方差已知，未知參數為$A$。

信號的模型為：

$$x[n] = A + w[n] \quad n=0, 1, \dots, N-1$$

參數$A$的兩個可能的估計器是：

- $$\hat{A}_1 = x[0]$$
- $\hat{A}_2 = \frac{1}{N} \sum_{n=0}^{N-1} x[n]$，即採樣平均（Sample mean）

通過計算兩個估計器的期望可以發現，它們的均值均為$A$：

$$\mathrm{E}\left[\hat{A}_1\right] = \mathrm{E}\left[ x[0] \right] = A$$

和

$$\mathrm{E}\left[ \hat{A}_2 \right] = \mathrm{E}\left[ \frac{1}{N} \sum_{n=0}^{N-1} x[n] \right] = \frac{1}{N} \left[ \sum_{n=0}^{N-1} \mathrm{E}\left[ x[n] \right] \right] = \frac{1}{N} \left[ N A \right] = A$$

兩個估計器的均值沒有差異，然而它們的方差不同：

$$\mathrm{var} \left( \hat{A}_1 \right) = \mathrm{var} \left( x[0] \right) = \sigma^2$$

和

$$\mathrm{var} \left( \hat{A}_2 \right) = \mathrm{var} \left( \frac{1}{N} \sum_{n=0}^{N-1} x[n] \right) = \frac{1}{N^2} \left[ \sum_{n=0}^{N-1} \mathrm{var}(x[n]) \right] = \frac{1}{N^2} \left[ N \sigma^2 \right] = \frac{\sigma^2}{N}$$

當$N>1$時，$\mathrm{var} \left( \hat{A}_1 \right) < \mathrm{var} \left( \hat{A}_2 \right)$，所以似乎採樣平均$\hat{A}_2$是一個更好的估計器。

### 最大似然估計

使用最大似然估計繼續上面的例子，噪聲在採樣點$w[n]$上的機率密度函數（pdf）為：

$$p(w[n]) = \frac{1}{\sigma \sqrt{2 \pi}} \exp\left(- \frac{1}{2 \sigma^2} w[n]^2 \right)$$

此時$x[n]$的機率為（$x[n]$服從分布$\mathcal{N}(A, \sigma^2)$）：

$$p(x[n]; A) = \frac{1}{\sigma \sqrt{2 \pi}} \exp\left(- \frac{1}{2 \sigma^2}(x[n] - A)^2 \right)$$

由於相互獨立，$\mathbf{x}$的機率為：

$$p(\mathbf{x}; A) = \prod_{n=0}^{N-1} p(x[n]; A) = \frac{1}{\left(\sigma \sqrt{2\pi}\right)^N} \exp\left(- \frac{1}{2 \sigma^2} \sum_{n=0}^{N-1}(x[n] - A)^2 \right)$$

對上式取自然對數：

$$\ln p(\mathbf{x}; A) = -N \ln \left(\sigma \sqrt{2\pi}\right) - \frac{1}{2 \sigma^2} \sum_{n=0}^{N-1}(x[n] - A)^2$$

於是最大似然估計器為：

$$\hat{A} = \arg \max \ln p(\mathbf{x}; A)$$

計算對數-最大似然函數的一階導數：

$$\frac{\partial}{\partial A} \ln p(\mathbf{x}; A) = \frac{1}{\sigma^2} \left[ \sum_{n=0}^{N-1}(x[n] - A) \right] = \frac{1}{\sigma^2} \left[ \sum_{n=0}^{N-1}x[n] - N A \right]$$

令其為0：

$$0 = \frac{1}{\sigma^2} \left[ \sum_{n=0}^{N-1}x[n] - N A \right] = \sum_{n=0}^{N-1}x[n] - N A$$

得到最大似然估計器：

$$\hat{A} = \frac{1}{N} \sum_{n=0}^{N-1}x[n]$$

它是一個簡單的採樣平均。從這個例子中可以發現，被獨立同分布的加性高斯白噪聲污染的、由未知常數構成的$N$點信號的最大似然估計其就是採樣平均。

### Cramér-Rao下限

為了找到採樣平均估計器的Cramér-Rao下限（CRLB），需要找到Fisher information數

$$\mathcal{I}(A) = \mathrm{E} \left( <pre> \left[ \frac{\partial}{\partial\theta} \ln p(\mathbf{x}; A) \right]^2 </pre> \right) = -\mathrm{E} \left[ <pre> \frac{\partial^2}{\partial\theta^2} \ln p(\mathbf{x}; A) </pre> \right]$$

從上面得到

$$\frac{\partial}{\partial A} \ln p(\mathbf{x}; A) = \frac{1}{\sigma^2} \left[ \sum_{n=0}^{N-1}x[n] - N A \right]$$

取二階導數

$$\frac{\partial^2}{\partial A^2} \ln p(\mathbf{x}; A) = \frac{1}{\sigma^2}(- N) = \frac{-N}{\sigma^2}$$

發現負的期望值是無關緊要的（trivial），因為它現在是一個確定的常數

$$-\mathrm{E} \left[ <pre> \frac{\partial^2}{\partial A^2} \ln p(\mathbf{x}; A) </pre> \right] = \frac{N}{\sigma^2}$$

最後，將Fisher information代入

$$\mathrm{var}\left( \hat{A} \right) \geq \frac{1}{\mathcal{I}}$$

得到

$$\mathrm{var}\left( \hat{A} \right) \geq \frac{\sigma^2}{N}$$

將這個值與前面確定的採樣平均的變化比較顯示對於所有的$N$和$A$來說採樣平均都是等於Cramér-Rao下限。

採樣平均除了是最大似然估計器之外還是最小變化無偏估計器（MVUE）。

這個直流增益 + WGN的例子是Kay的統計信號處理基礎中一個例子的再現。

## 相關書籍

- Fundamentals of Statistical Signal Processing: Estimation Theory by Steven M. Kay (ISBN 0-13-345711-7)
- An Introduction to Signal Detection and Estimation by H. Vincent Poor (ISBN 0-38-794173-8)
- Detection, Estimation, and Modulation Theory, Part 1 by Harry L. Van Trees (ISBN 0-47-109517-6; website)