\[
\newcommand{\bs}{\boldsymbol}
\newcommand{\bsX}{\boldsymbol{X}}
\newcommand{\bf}{\mathbf}
\newcommand{\msc}{\mathscr}
\newcommand{\mca}{\mathcal}
\newcommand{\T}{\text{T}}
\newcommand{\rme}{\mathrm{e}}
\newcommand{\rmi}{\mathrm{i}}
\newcommand{\rmj}{\mathrm{j}}
\newcommand{\rmd}{\mathrm{d}}
\newcommand{\rmm}{\mathrm{m}}
\newcommand{\rmb}{\mathrm{b}}
\newcommand{\and}{\land}
\newcommand{\or}{\lor}
\newcommand{\exist}{\exists}
\newcommand{\sube}{\subseteq}
\newcommand{\lr}[3]{\left#1 #2 \right#3}
\newcommand{\intfy}{\int_{-\infty}^{+\infty}}
\newcommand{\sumfy}[1]{\sum_{#1=-\infty}^{+\infty}}
\newcommand{\vt}{\vartheta}
\newcommand{\ve}{\varepsilon}
\newcommand{\vp}{\varphi}
\newcommand{\Var}{\text{Var}}
\newcommand{\Cov}{\text{Cov}}
\newcommand{\edef}{\xlongequal{def}}
\newcommand{\prob}{\text{P}}
\newcommand{\Exp}{\text{E}}
\newcommand{\t}[1]{\text#1}
\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\versionofnewcommand}{\text{260125}}
\]
Cramer-Rao Lower Bound (CRLB)
这是科学研究的另一种思路. 一些人利用各种优化过程由一个估计的构造不断寻找 MVUE, 一些人在寻求一个不依赖于具体的估计构造的下界. 这样就把底托住了, 或者说勾勒出我们讨论范围的轮廓, 这当然是非常 fundamental 的.
CRLB 是对估计的方差划定了一个下界, 且这个下界不依赖于具体的估计构造, 仅依赖于模型本身, 记为 \(I^{-1}(\theta)\). 也即:
\[
\forall \hat{\theta}, E(\hat{\theta})=\theta \Rightarrow Var(\hat{\theta})\geq I^{-1}(\theta)
\]
该不等式称为 Cramer-Rao 不等式, \(I^{-1}(\theta)\) 就是 Cramer-Rao 下界. 其中 \(I(\theta)\) 称为 Fisher Information.
\[
I(\theta)=E\left(\frac{\partial}{\partial\theta}\log f(X,\theta)\right)^{2}=-E\left( \frac{\partial^2}{\partial\theta^2}\log f(X,\theta) \right)
\]
注意, 这里我没写 \(x\), 而是写 \(X\), 也即把随机变量代入到分布模型里, 使得括号里那一坨依然是个随机变量, 从而取期望这件事才有意义.
来解释一下. 首先, 为了将 MSE 与 \(\frac{\partial}{\partial\theta} f(x,\theta)\) 这两者联系起来, 我们需要在一个式子里同时纳入 \((\hat{\theta}-\theta)\) 和 \(\frac{\partial}{\partial\theta} f(x,\theta)\) . 于是, 对:
\[
f(x,\theta)\geq0, \ \int_{-\infty}^{+\infty}f(x,\theta)dx=1
\]
两边求导, 有:
\[
\frac{\partial}{\partial\theta}\int_{-\infty}^{+\infty}f(x,\theta)dx=0
\]
根据工科思维, 不严谨地交换积分微分次序:
\[
\int_{-\infty}^{+\infty}\frac{\partial}{\partial\theta}f(x,\theta)dx=0 \tag{1}
\]
其次,
根据无偏估计:
\[
E(\hat{\theta})=\int_{-\infty}^{+\infty}\hat{\theta}\cdot f(x,\theta)dx=\theta
\]
两边求导, 再交换一下次序, 于是有:
\[
\int_{-\infty}^{+\infty}\hat{\theta}\cdot\frac{\partial}{\partial\theta} f(x,\theta)dx=1 \tag{2}
\]
比较 \((1)\) 式和 \((2)\) 式, 我们做这样的操作, 给 \((1)\) 式两端都乘上 \(\theta\) :
\[
\int_{-\infty}^{+\infty}\theta\cdot\frac{\partial}{\partial\theta} f(x,\theta)dx=0 \tag{3}
\]
这样就对齐了! 于是 \((2)-(3)\) 立刻得出:
\[
\int_{-\infty}^{+\infty}(\hat{\theta}-\theta)\cdot\frac{\partial}{\partial\theta} f(x,\theta)dx=1
\]
Fisher 天才的地方就在这里了, 做一个 \(\log\) 的变换, 把 \(f(x,\theta)\) 从微分里解救出来, 从而构造出求期望的形式.
\[
\int_{-\infty}^{+\infty}(\hat{\theta}-\theta)\cdot f(x,\theta)\cdot\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)dx=1
\]
利用 Cauchy-Schwarz 不等式,
\[
\begin{aligned}
1=&\int_{-\infty}^{+\infty}(\hat{\theta}-\theta)\cdot f(x,\theta)\cdot\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)dx\\
=& \int_{-\infty}^{+\infty}(\hat{\theta}-\theta)\sqrt{f(x,\theta)}\cdot\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)\sqrt{f(x,\theta)}\ dx \\
\leq& \left( \int_{-\infty}^{+\infty}(\hat{\theta}-\theta)^2f(x,\theta)\ dx\cdot\int_{-\infty}^{+\infty}\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)^2f(x,\theta)\ dx\right)^{\frac{1}{2}}
\end{aligned}
\]
平方消掉开根:
\[
\begin{aligned}
1\leq& \int_{-\infty}^{+\infty}(\hat{\theta}-\theta)^2f(x,\theta)\ dx\cdot\int_{-\infty}^{+\infty}\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)^2 f(x,\theta)\ dx\\
=& E(\hat{\theta}-\theta)^2 \cdot E\left(\frac{\partial}{\partial\theta}\log f(X,\theta)\right)^2
\end{aligned}
\]
因此有！
\[
MSE=Var(\hat{\theta})=E(\hat{\theta}-\theta)^2\geq 1/E\left(\frac{\partial}{\partial\theta}\log f(X,\theta)\right)^{2} = I^{-1}(\theta)
\]
其中,
\[
I(\theta)=E\left(\frac{\partial}{\partial\theta}\log f(X,\theta)\right)^{2}
\]
进一步地,
\[
1=\int_{-\infty}^{+\infty} f(x,\theta)dx \Rightarrow 0
=\int_{-\infty}^{+\infty} \frac{\partial}{\partial\theta}f(x,\theta)dx
=\int_{-\infty}^{+\infty} \left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)\cdot f(x,\theta)\ dx
\]
再求导, 则有:
\[
\begin{aligned}
0=&\int_{-\infty}^{+\infty}
\frac{\partial}{\partial\theta}
\left(\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)\cdot f(x,\theta)\right) dx\\
=&\int_{-\infty}^{+\infty}\left(\frac{\partial^2}{\partial\theta^2}
\log f(x,\theta)\right)\cdot f(x,\theta)\ dx
+\int_{-\infty}^{+\infty}\left(\frac{\partial}{\partial\theta}
\log f(x,\theta)\right)\left(\frac{\partial}{\partial\theta}f(x,\theta)\right)dx\\
=& E\left( \frac{\partial^2}{\partial\theta^2}\log f(X,\theta) \right)
+\int_{-\infty}^{+\infty}\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)^2 f(x,\theta)\ dx\\
=& E\left( \frac{\partial^2}{\partial\theta^2}\log f(X,\theta) \right)
+ E\left(\frac{\partial}{\partial\theta}\log f(X,\theta)\right)^2
\end{aligned}
\]
从而这两坨互为相反数, 即:
\[
I(\theta)=E\left(\frac{\partial}{\partial\theta}\log f(X,\theta)\right)^{2}=-E\left( \frac{\partial^2}{\partial\theta^2}\log f(X,\theta) \right)
\]
回顾一下刚才使用的柯西不等式, 注意等号成立的条件是两项只差一个常数倍. 即:
\[
\frac{\partial}{\partial\theta}\log f(x,\theta)=K\cdot(\hat\theta-\theta)
\]
其中, \(K\) 不依赖于 \(X\). 吔, 但是 \(K\) 依赖于 \(\theta\) 是一点问题没有的! 毕竟 \(\theta\) 也是一个常数. 并且实践中, 我们也更多地写成这样:
\[
\frac{\partial}{\partial\theta}\log f(x,\theta)=K(\theta)\cdot(\hat\theta-\theta)
\]
此时, Cramer-Rao 下界可以达到. 这就是最好的估计了.
先举个例子体验一下吧:
假设有高斯分布, \(\sigma^2\) 已知:
\[
f(x,\theta)=\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{(x-\theta)^2}{2\sigma^2})
\]
一般来讲, 我们当然希望 \(\sigma\) 越小越好, 因为这样只需要取样距离 \(\theta\) 有微小的变化, 就会使结果有很大不同, 这对于统计过程来说当然是有利的. 也即, 我们希望曲率 (curvature) 大一点. 基于这一动机, 我们来算算二阶导. 为了方便, 取对数计算一下:
\[
\begin{aligned}
\log (\sqrt{2\pi}\sigma\cdot f(x,\theta)) &=-\frac{(x-\theta)^2}{2\sigma^2}\\
\frac{\partial}{\partial\theta}\left(-\frac{(x-\theta)^2}{2\sigma^2}\right) &=\frac{1}{\sigma^2}(x-\theta)\\
\frac{\partial^2}{\partial\theta^2}\left(-\frac{(x-\theta)^2}{2\sigma^2}\right) &=-\frac{1}{\sigma^2}
\end{aligned}
\]
首先来看一阶导的结果, 容易发现, 这压根就是刚才讨论的柯西不等式取等条件的结果: \(\frac{\partial}{\partial\theta}\log f(x,\theta)=K(\theta)\cdot(\hat\theta-\theta)\). 这样来看, \(x\) 本身就已经是 \(\theta\) 最好的估计了, 这当然非常符合我们的直观.
进一步地, 很容易得出高斯分布的 Fisher Information
\[
I(\theta)=-E\left(\frac{\partial^2}{\partial\theta^2}\log f(x,\theta)\right)=\frac{1}{\sigma^2}
\]
Cramer-Rao 界自然就是 \(\sigma^2\) 了.
诶, 等等! 很巧的是, 刚才算出的 \(K(\theta)\) 在高斯分布的例子中等于 \(I(\theta)\). 这难道是恒成立的? 事实上确实如此: 对取等条件式
\[
\frac{\partial}{\partial\theta}\log f(x,\theta)=K(\theta)\cdot(\hat\theta-\theta)
\]
两边平方, 再取期望:
\[
\begin{aligned}
\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)^2 & =K^2(\theta)\cdot(\hat\theta-\theta)^2\\
E\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)^2 &=E\left(K^2(\theta)\cdot(\hat\theta-\theta)^2\right)\\
&=K^2(\theta)\cdot E(\hat\theta-\theta)^2\\
\end{aligned}
\]
而取等条件满足时, 方差就是\(I^{-1}(\theta)\), 即:
\[
I(\theta)=K^2(\theta)\cdot I^{-1}(\theta)
\]
这个巧妙, 这样 \(K(\theta)\) 压根就是 \(I(\theta)\) 了.
接着再考虑一个例子:
假设随机变量\(X_1,\cdots,X_n\sim_{i.d.d.} N(\theta,\sigma^2)\). 对于
\[
f(x_1,\cdots,x_n)=\left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n\cdot\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\theta)^2\right)
\]
我们来算算它的 Cramer-Rao 下界:
\[
\begin{aligned}
\log f(x_1,\cdots,x_n) &=-n\log(\sqrt{2\pi}\sigma)-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\theta)^2\\
\frac{\partial}{\partial\theta}\log f(x_1,\cdots,x_n) &=\frac{1}{\sigma^2}\sum_{i=1}^n(x_i-\theta)
=\frac{n}{\sigma^2}\left({\frac{\sum_{i=1}^nx_i}{n}}-\theta\right)
\end{aligned}
\]
写出来这个形式, 就足够说明样本均值就是此例的 MVUE, 而且 Fisher Information 和 Cramer-Rao 下界也出来了.
进一步延伸, 当估计并非无偏:
如果 \(E(\hat\theta)=g(\theta)\), Cramer-Rao 界会有什么变化?
\[
\begin{aligned}
g(\theta) &=\int_{-\infty}^{+\infty}{\hat\theta}(x)f(x,\theta)dx\\
\Rightarrow g\prime(\theta) &=\int_{-\infty}^{+\infty}\hat\theta(x)\frac{\partial}{\partial\theta}f(x,\theta)dx
\end{aligned}
\]
比较一下之前得到的 (1) 式和 (2) 式, 容易发现:
\[
g\prime(\theta)=\int_{-\infty}^{+\infty}(\hat\theta-\theta)\frac{\partial}{\partial\theta}f(x,\theta)dx
=\int_{-\infty}^{+\infty}(\hat\theta-\theta)f(x,\theta)\frac{\partial}{\partial\theta}\log f(x,\theta)dx
\]
从而有
\[
\begin{aligned}
g^\prime(\theta) &=E\left( (\hat\theta-\theta)\frac{\partial}{\partial\theta}\log f(x,\theta) \right)\\
\Rightarrow (g^{\prime}(\theta))^2 &\leq E(\hat\theta-\theta)^2\cdot E\left(\frac{\partial}{\partial\theta}\log f(x,\theta)\right)^2
\end{aligned}
\]
这样, 就有了新的 Cramer-Rao 下界喽.
\[
E(\hat\theta-\theta)^2 \geq (g^{\prime}(\theta))^2\cdot I^{-1}(\theta)
\]
再一个延伸, 当参数变多了以后:
首先有个问题, 当参数变多了, 我们该如何理解 MSE? 很自然地, 可以把它写成向量形式, 在多维空间里考虑. 也即 \(\boldsymbol\theta,\boldsymbol{\hat\theta}\in \mathbb{R}^m\).
一个自然的考虑是, 把期望写成 \(\sum_{i=1}^m(\hat\theta_i-\theta_i)^2\), 但这样有个重要的问题, 就是无法反映各个 \(\hat\theta_i-\theta\) 之间的关系. 因此, 我们考虑用矩阵这一手段表现其丰富的内涵.
定义 协方差矩阵:
\[
C(\boldsymbol{\hat\theta}-\boldsymbol\theta)=E(\boldsymbol{\hat\theta}-\boldsymbol\theta)(\boldsymbol{\hat\theta}-\boldsymbol\theta)^\mathrm T
\]
其中:
\[
\boldsymbol{\hat\theta}-\boldsymbol\theta=(\hat\theta_1-\theta_1,\cdots,\hat\theta_m-\theta_m)^\mathrm T
\]
有些地方会写作 \(C(\boldsymbol{\hat\theta})=E(\boldsymbol{\hat\theta}-E(\boldsymbol{\hat\theta}))(\boldsymbol{\hat\theta}-E(\boldsymbol{\hat\theta}))^\text T\), 其实都一样的.
但随之而来又有一个问题, 所谓 "界" 是要比大小的, 一个矩阵如何比大小咧? 试图回忆一下, 有个叫 "正定性" 的东西:
对于 \(A \in \mathbb{R}^{m\times m}\), 若任取非零向量 \(\boldsymbol x\in \mathbb{R}^m\), \(\boldsymbol x^\mathrm TA\boldsymbol x\geq0\), 则称矩阵 \(A\) 是正定的. 一般记作 \(A\geq0\). 那么, 通过正定性就能跟0 "比大小", 从而能让两个矩阵 "比大小" 了.
由此引入 Fisher Information Matrix:
\[
\exist\ I(\boldsymbol\theta)\in\mathbb{R}^{m\times m}, s.t. C(\boldsymbol{\hat\theta}-\boldsymbol\theta)-I^{-1}(\boldsymbol\theta)\geq0
\]
约定俗成地, 一般记作
\[
C(\boldsymbol{\hat\theta}-\boldsymbol{\theta})\geq I^{-1}(\boldsymbol\theta)
\]
自然地, 想到对 \(\log f(\boldsymbol x,\boldsymbol\theta)\in \mathbb{R}^m\) 的每一维求导, 也即求梯度:
\[
\nabla_\theta\ \log f(\boldsymbol x,\boldsymbol\theta)=(\frac{\partial}{\partial\theta_1}\log f,\cdots,\frac{\partial}{\partial\theta_m}\log f)^\mathrm T
\]
接着自然地猜到,
\[
I(\boldsymbol\theta)=E(\nabla_\boldsymbol\theta\ \log f)(\nabla_\boldsymbol\theta\ \log f)^\mathrm T=-E(\nabla_\boldsymbol\theta^2\ \log f(\boldsymbol x,\boldsymbol\theta))
\]
这里, \(\nabla_\theta^2\ \log f(\boldsymbol x,\boldsymbol \theta)\) 是个 Hessian Matrix.
\[
\nabla^2_\boldsymbol\theta\ g=\left(\frac{\partial^2}{\partial\theta_i\cdot\partial\theta_j}g\ \right)_{i,j}
\]
现在开始尝试证明!
Lemma: \(\forall \boldsymbol X\in \mathbb{R}^n\), \(E(\boldsymbol X\boldsymbol X^\mathrm T)\geq 0\). 即, 任何的相关矩阵都是正定的.
Proof: 只需证: \(\forall \boldsymbol\alpha\in \mathbb{R}^n\), \(\boldsymbol\alpha^\mathrm T E(\boldsymbol X\boldsymbol X^\mathrm T)\boldsymbol\alpha=0\).
\[
\boldsymbol\alpha^T E(\boldsymbol X\boldsymbol X^T)\boldsymbol\alpha=E(\boldsymbol\alpha^T\boldsymbol X\boldsymbol X^T\boldsymbol\alpha)=E((\boldsymbol\alpha^T\boldsymbol X)(\boldsymbol X^T\boldsymbol\alpha))=E(\boldsymbol\alpha^T\boldsymbol X)^2\geq0
\]
注意, 因为 \(\boldsymbol\alpha\) 不具有随机性, 所以能随意地放到期望里.
接下来, 对于 \(\boldsymbol\theta\in\mathbb{R}^m\), 作估计 \(\boldsymbol{\hat\theta}=(\hat\theta_1,\cdots,\hat\theta_m)^T\in\mathbb{R}^m\). 另有 \(\nabla_\boldsymbol\theta\ \log f(\boldsymbol x,\boldsymbol\theta)=(\frac{\partial}{\partial\theta_1}\log f,\cdots,\frac{\partial}{\partial\theta_m}\log f)^\mathrm T\). Rao 有个好主意, 构造:
\[
X=(\hat\theta_1-\theta_1,\cdots,\hat\theta_m-\theta_m,\frac{\partial}{\partial\theta_1}\log f,\cdots,\frac{\partial}{\partial\theta_m}\log f)^T\in \mathbb{R}^{2m}
\]
来算它的相关阵:
\[
E(XX^T)=
\begin{bmatrix}
C(\hat\theta,\theta) & H \\
H^T & I(\theta)
\end{bmatrix}
\]
计算一下 \(H\) 到底是个啥:
\[
\begin{aligned}
H_{ij}&=E\left((\hat\theta_i-\theta_i)\frac{\partial}{\partial\theta_j}\log f(x,\theta)\right)\\
&=\int_{-\infty}^{+\infty}(\hat\theta_i-\theta_i)\frac{\partial}{\partial\theta_j}\log f(x,\theta)\cdot f(x,\theta)\ dx\\
&=\int_{-\infty}^{+\infty}(\hat\theta_i-\theta_i)\frac{\partial}{\partial\theta_j}f(x,\theta)\ dx\\
&=\int_{-\infty}^{+\infty}\hat\theta_i\frac{\partial}{\partial\theta_j}f(x,\theta)\ dx-\theta_i\cdot\frac{\partial}{\partial\theta_j}\int_{-\infty}^{+\infty}f(x,\theta)\ dx\\
&=\int_{-\infty}^{+\infty}\hat\theta_i\frac{\partial}{\partial\theta_j}f(x,\theta)\ dx-\theta_i\cdot\frac{\partial}{\partial\theta_j}1\\
&=\int_{-\infty}^{+\infty}\hat\theta_i\frac{\partial}{\partial\theta_j}f(x,\theta)\ dx\\
&=\frac{\partial}{\partial\theta_j}\int_{-\infty}^{+\infty}\hat\theta_i\cdot f(x,\theta)\ dx
\end{aligned}
\]
记
\[
g_i(\theta)=\int_{-\infty}^{+\infty}\hat\theta_i\cdot f(x,\theta)\ dx
\]
则有
\[
E(\hat\theta)=g(\theta)=(g_1(\theta),\cdots,g_m(\theta))
\]
因此,
\[
H_{ij}=\frac{\partial}{\partial\theta_j}g_i(\theta)
\]
tips: 当估计无偏时, \(H\) 就是个单位矩阵了.
回到原来的那个大矩阵
\[
E(XX^T)=
\begin{bmatrix}
C(\hat\theta,\theta) & H \\
H^T & I(\theta)
\end{bmatrix}\\
\]
构造一个矩阵 ( \(A\) 待定 )
\[
\begin{bmatrix}
I & A\\
0 & I
\end{bmatrix}
\]
我们希望通过调整 \(A\) 的结构, 使得该矩阵和 \(E(XX^T)\) 相乘后, 可以让第一象限变成 0. 于是计算:
\[
H+A\cdot I(\theta)=0\Rightarrow A=-HI^{-1}(\theta)
\]
这样, 相乘后的结果就是:
\[
\begin{bmatrix}
I & A\\
0 & I
\end{bmatrix}
\begin{bmatrix}
C(\hat\theta,\theta) & H \\
H^T & I(\theta)
\end{bmatrix}
=
\begin{bmatrix}
C(\hat\theta-\theta)-HI^{-1}(\theta)H^T & 0\\
H^T & I(\theta)
\end{bmatrix}
\]
再构造一个矩阵 ( \(B\) 待定 ), 让刚得到的大矩阵被右乘后, 第三象限也变成 0.
\[
\begin{bmatrix}
I & 0\\
B & I
\end{bmatrix}
\]
\[
H^T+I(\theta)\cdot B=0 \Rightarrow B=-I^{-1}(\theta)H^T
\]
于是有:
\[
\begin{bmatrix}
I & -HI^{-1}(\theta)\\
0 & I
\end{bmatrix}
\begin{bmatrix}
C(\hat\theta,\theta) & H \\
H^T & I(\theta)
\end{bmatrix}
\begin{bmatrix}
I & 0\\
-I^{-1}(\theta)H^T & I
\end{bmatrix}=
\begin{bmatrix}
C(\hat\theta-\theta)-HI^{-1}(\theta)H^T & 0\\
0 & I(\theta)
\end{bmatrix}
\]
通过向量乘法的概念, 可以发现, 两边的矩阵互为转置, 这也就意味着对中间的矩阵进行了一次合同变换. 而这并不会对正定性有影响. 由于等号左侧中间的矩阵是正定的, 因此等号右边的分块对角矩阵也是正定的. 又: 分块对角矩阵正定 \(\Leftrightarrow\) 每个对角块正定. 因此:
\[
C(\hat\theta-\theta)-HI^{-1}(\theta)H^T \geq 0
\]
假如估计是无偏的, \(H\) 就是单位矩阵, 此时自然就有
\[
C(\hat\theta-\theta)\geq I^{-1}(\theta)
\]
该证明由 Rao 精妙地给出, 不仅给出了无偏情况的证明, 还更扩展地给出了非无偏的结果.