统计参数估计

统计参数估计

在进行统计决策时，我们假设后验概率是已知的，这往往不现实。在有些情况下，我们假设数据近似服从某种分布(最常见的如高斯分布)，估计数据的分布参数。

最大似然估计

我们一般假设所观测到的变量是最有可能出现的，因此，直观上最大似然估计就是估计出使分布函数最契合所观测到数据分布的参数。形式化的语言表示为：选择θ，使得：
$$\hat{\theta }=\mathrm{argmax}[\mathrm{p}(\mathrm{X}|\theta )]$$
上述公式对应只有一个观测变量的情况，当有多个观测变量时，最大似然估计的目标是使每个观测变量概率都最大，也就是这些变量的概率乘积最大化：
$$\hat{\theta }=\mathrm{argmax}\left[\log \prod _{\mathrm{k}=1}^{\mathrm{N}}\mathrm{p}\left({\mathrm{x}}^{(\mathrm{k}}|\theta \right)\right]=\mathrm{argmax}\left[\sum _{\mathrm{k}=1}^{\mathrm{N}}\log \mathrm{p}\left({\mathrm{x}}^{(\mathrm{k}}|\theta \right)\right]$$
通过一个取对数操作，就讲乘积变成了加法，这样求导就方便了很多。最大似然估计，其实就是一个求最优化的问题。

多元高斯分布为例

对于一个多元高斯分布(就是有多个变量的高斯分布)：
$$p(x;\mu ,\Sigma )=\frac{1}{(2\pi {)}^{\frac{D}{2}}|\Sigma {|}^{\frac{1}{2}}}\exp \left(-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )\right)$$
其中：
$$\mu =\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ {\mu }_N\end{array}\right]$$

$$\Sigma =\frac{1}{m}{X}^{T}X$$

分别为均值和协方差。对于N个观测变量，最大似然估计需要估计均值和协方差使得：
$$E(\theta )=\sum _{k=1}^N\log p\left(x_k|\theta \right)=\sum _{k=1}^N-\log (2\pi {)}^{\frac{D}{2}}-\frac{1}{2}\log |\Sigma |-\frac{1}{2}{\left(x_k-\mu \right)}^T{\Sigma }^{-1}\left(x_k-\mu \right)$$
最大，即：
$$\hat{\theta }=\arg \max E(\theta )$$
1.首先估计均值。对μ求偏导：
$$x\frac{\partial E(\theta )}{\partial \mu }=\sum _{k=1}^N{\Sigma }^{-1}\left(x_k-\mu \right)=0(注:\frac{\partial (X^{T}AX)}{\partial X}=(A+A^T)X,{\Sigma }^{-1^T}={\Sigma }^{-1})$$

$$\mu =\frac{1}{N}\sum _{k=1}^N{x}_k$$

2.接下来估计协方差，对Sigma求偏导(实际上是对Sigma的逆求偏导，这并不影响结果，但可以是推导更简便)：
$$\frac{\partial E(\theta )}{\partial {\Sigma }^{-1}}=\frac{\partial \left(\sum _{k=1}^N\frac{1}{2}\log \left|{\Sigma }^{-1}\right|-\frac{1}{2}{\left(x_k-\mu \right)}^T{\mathbf{\Sigma }}^{-1}\left(x_k-\mu \right)\right)}{\partial {\Sigma }^{-1}}$$
上式求解有些复杂，分别对式中两项拆开进行求解：

第一项：
$$对于对称矩阵A,有:\left|A^{-1}\right|={\left|A\right|}^{-1};\frac{\partial \log \left|A\right|}{\partial A}=2A^{-1}-Diag(A^{-1})$$

$$于是就有：\frac{\partial \log \left|{\Sigma }^{-1}\right|}{\partial {\Sigma }^{-1}}=2\Sigma -Diag(\Sigma )$$

第二项：
$$\begin{array}{l}对于对称矩阵A、T，有：X^{T}AX=Tr(AX{X}^T);\frac{\partial Tr(AB)}{\partial A}=B+B^T-Diag(B);\\ 从而可以得出：\frac{\partial X^{T}AX}{\partial A}=2X{X}^T-Diag(X{X}^T)\end{array}$$

$$于是就有：\frac{\partial {\left(x_k-\mu \right)}^T{\Sigma }^{-1}\left(x_k-\mu \right)}{\partial {\Sigma }^{-1}}=2\left(x_k-\mu \right){\left(x_k-\mu \right)}^T-Diag(\left(x_k-\mu \right){\left(x_k-\mu \right)}^T)$$

综合这两项有：
$$\begin{array}{l}\frac{\partial E(\theta )}{\partial {\Sigma }^{-1}}=\frac{1}{2}\sum _{k=1}^{N}2\Sigma -2\left(x_k-\mu \right){\left(x_k-\mu \right)}^T-Diag(\Sigma )+Diag(\left(x_k-\mu \right){\left(x_k-\mu \right)}^T)\\ =\frac{1}{2}\sum _{k=1}^{N}2(\Sigma -\left(x_k-\mu \right){\left(x_k-\mu \right)}^T)-Diag(\Sigma -\left(x_k-\mu \right){\left(x_k-\mu \right)}^T)=0\\ （注:直观可得Diag(A)+Diag(B)=Diag(A+B)）\end{array}$$
求解这个优化问题可得：
$$\sum _{k=1}^N\Sigma -\left(x_k-\mu \right){\left(x_k-\mu \right)}^T=0\Rightarrow \Sigma =\frac{1}{N}\sum _{k=1}^N\left(x_k-\mu \right){\left(x_k-\mu \right)}^T$$

估计的偏差和方差

这部分内容部分参考了知乎@大海啊你全是水的文章。

估计的偏差

所谓偏差就是估计值的期望和真实值的差距。公式描述为：
$$\mathrm{bias}\left({\hat{\mathbf{\theta }}}_m\right)=\mathbb{E}\left({\hat{\mathbf{\theta }}}_m\right)-\mathbf{\theta }$$
以高斯分布为例。我们知道高斯分布的参数估计分别为：
$$\mu =\bar{x},{\sigma }^2=\frac{1}{n}\sum _{i=1}^n{\left(x_i-\mu \right)}^2$$
推导过程可见@我是8位的 博客。

均值估计的偏差为：
$$\begin{array}{rl}\mathrm{bias}\left({\hat{\mu }}_m\right)& =\mathbb{E}\left({\hat{\mu }}_m\right)-\mu \\ & =\mathbb{E}\left[\frac{1}{m}\sum _{i=1}^m{x}^{(i)}\right]-\mu \\ & =\left(\frac{1}{m}\sum _{i=1}^m\mathbb{E}\left[x^{(i)}\right]\right)-\mu \\ & =\left(\frac{1}{m}\sum _{i=1}^m\mu \right)-\mu \\ & =\mu -\mu =0\end{array}$$
所以均值是无偏估计。方差估计的偏差为：
$$\mathrm{bias}\left({\hat{\sigma }}_m^2\right)=\mathbb{E}\left[{\hat{\sigma }}_m^2\right]-{\sigma }^2=\frac{m-1}{m}{\sigma }^2-{\sigma }^2=-\frac{{\sigma }^2}{m}$$
这个期望E是这样计算来的：
$$\begin{array}{rl}\mathbb{E}\left[{\hat{\sigma }}_m^2\right]& =\mathbb{E}\left[\frac{1}{m}\sum _{i=1}^m{\left(x^{(i)}-{\hat{\mu }}_m\right)}^2\right]\\ & =\mathbb{E}\left\{\frac{1}{m}\sum _{i=1}^m\left[{\left(x^{(i)}\right)}^2-2x^{(i)}{\hat{\mu }}_m+{\hat{\mu }}_m^2\right]\right\}\\ & =\mathbb{E}\left\{\frac{1}{m}\sum _{i=1}^m\left[{\left(x^{(i)}\right)}^2\right]-\frac{2{\hat{\mu }}_m}{m}\sum _{i=1}^m\left[x^{(i)}\right]+\frac{1}{m}\sum _{i=1}^m\left({\hat{\mu }}_m^2\right)\right\}\\ & =\mathbb{E}\left\{\frac{1}{m}\sum _{i=1}^m\left[{\left(x^{(i)}\right)}^2\right]-2{\hat{\mu }}_m^2+{\hat{\mu }}_m^2\right\}\\ & =\mathbb{E}\left\{\frac{1}{m}\sum _{i=1}^m\left[{\left(x^{(i)}\right)}^2\right]-{\hat{\mu }}_m^2\right\}\\ & =\frac{1}{m}\sum _{i=1}^m\mathbb{E}\left[{\left(x^{(i)}\right)}^2\right]-\mathbb{E}\left({\hat{\mu }}_m^2\right)\end{array}$$

$$根据方差定义:\mathrm{Var}(x)=\mathbb{E}\left\{{[x-\mathbb{E}(x)]}^2\right\}=\mathbb{E}\left\{x^2-2x\mathbb{E}(x)+\mathbb{E}{[(x)]}^2\right\}=\mathbb{E}(x^2)-[\mathbb{E}(x){]}^2$$

$$\mathbb{E}\left(x^2\right)=\mathrm{Var}(x)+[\mathbb{E}(x){]}^2$$

故有：
$$\begin{array}{rl}\mathbb{E}\left\{{\left[x^{(i)}\right]}^2\right\}=& \mathrm{Var}\left[x^{(i)}\right]+{\left\{\mathbb{E}\left[x^{(i)}\right]\right\}}^2={\sigma }^2+{\mu }^2\\ \mathbb{E}\left({\mu }_m^2\right)& =\mathrm{Var}\left[{\mu }_m\right]+{\left[\mathbb{E}\left({\mu }_m\right)\right]}^2=\mathrm{Var}\left\{\frac{1}{m}\sum _{i=1}^m\left[x^{(i)}\right]\right\}+{\mu }^2\\ & =\frac{1}{m^2}\sum _{i=1}^m\mathrm{Var}\left[x^{(i)}\right]+{\mu }^2\\ & =\frac{{\sigma }^2}{m}+{\mu }^2\end{array}$$
最终可得方差的期望：
$$\mathbb{E}\left[{\hat{\sigma }}_m^2\right]={\sigma }^2+{\mu }^2-\frac{{\sigma }^2}{m}-{\mu }^2=\frac{m-1}{m}{\sigma }^2$$
所以方差是有偏估计。

估计的方差

方差顾名思义即：
$$\mathrm{Var}(\hat{\theta })$$
在大多数情况下，偏差和方差就行家庭和事业不可得兼：

贝叶斯估计

最大似然估计把待估计参数θ看作是一个确定但未知的参数，而贝叶斯估计把参数θ看作为一个随机变量。我们可能对θ的分布有一个粗略的认识，即我们可能粗略知道一个先验p(θ)，通过利用观测到的样本使我们能够得出θ的后验概率密度p(θ|X)。

下面以简单的高斯分布为例。

先给出高斯分布公式：
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$
假设我们观测到了下面的一系列样本，它们来自一个方差已知的高斯分布，我们需要估计这个分布的均值μ。
$$X=\left\{x^1,x^2,\dots x^N\right\}$$
并且我们粗略地知道μ服从一个高斯分布(即先验)：
$${\mathrm{p}}_0(\theta )=\frac{1}{\sqrt{2\pi }{\sigma }_0}\exp \left(-\frac{1}{2{\sigma }_0^2}{\left(\theta -{\mu }_0\right)}^2\right)$$
根据贝叶斯公式就可以获得μ的后验概率密度p(θ|X)：
$$\begin{array}{l}p(\theta |X)=\frac{p(X|\theta )p_0(\theta )}{p(X)}\\ =\frac{p_0(\theta )}{p(X)}\prod _{k=1}^{N}p\left(x^{(k}|\theta \right)\\ =\frac{1}{\sqrt{2\pi }{\sigma }_0}\exp \left(-\frac{1}{2{\sigma }_0^2}{\left(\theta -{\mu }_0\right)}^2\right)\frac{1}{p(X)}\prod _{k=1}^N\left[\frac{1}{\sqrt{2\pi \sigma }}\exp \left(-\frac{1}{2{\sigma }^2}{\left(x^{(k}-\theta \right)}^2\right)\right]\end{array}$$
我们选择参数θ时，一般选择后验概率最大的那个。即：
$$x\hat{\theta }=\arg \max p(\theta |X)$$
通过求解最优化问题获得参数θ。接着上面的例子，对p(θ|X)取对数后求偏导(取对数是为了避免连乘求导的困难)：
$$\frac{\partial }{\partial \theta }\mathrm{logp}(\theta |X)=\frac{\partial }{\partial \mu }\left(-\frac{1}{2{\sigma }_0^2}{\left(\mu -{\mu }_0\right)}^2+\sum _{k=1}^N\frac{-1}{2{\sigma }^2}{\left(x^{(k}-\mu \right)}^2\right)=0$$

$${\mu }_{\mathrm{N}}=\underset{\mathrm{P}\mathrm{R}\mathrm{I}\mathrm{O}\mathrm{R}}{\underbrace{\frac{{\sigma }^2}{{\sigma }^2+N{\sigma }_0^2}{\mu }_0}}+\underset{\mathrm{M}\mathrm{A}\mathrm{X}\mathrm{I}\mathrm{M}\mathrm{U}\mathrm{M}\mathrm{L}\mathrm{I}\mathrm{K}\mathrm{E}\mathrm{L}\mathrm{H}\mathrm{O}\mathrm{O}\mathrm{D}}{\underbrace{\frac{\mathrm{N}{\sigma }_0^2}{{\sigma }^2+\mathrm{N}{\sigma }_0^2}\frac{1}{N}\sum _{k=1}^N{x}^{(k}}}$$

其中，前面项是先验项，后面是似然项。当观测变量比较少时，先验信息为主要贡献；当观测变量比较多时，似然信息为主要贡献。