.. _ch_estimate: ######################################################## 参数估计 ######################################################## 在上一章节我们介绍了概率论的一些基础知识，我们知道了什么是随机变量，什么是概率分布。 在有些场景下，我们知道一个随机变量 :math:`X` 服从什么概率分布 :math:`p(X;\theta)` ，但是这个概率分布存在一个未知的参数 :math:`\theta` 。 比如高斯分布中的均值 :math:`\mu` 或者方差 :math:`\sigma^2` 可能是未知的，亦或者两者都未知。 但是如果我们能得到这个随机变量的一些观测数据，我们就可以从这些数据中"学习(learning)"或者说"估计(estimate)" 出模型的未知参数，这些观测数据一般称为"训练数据(training data)"。 本章我们讨论如何利用这个概率分布的观测样本估计出这个未知参数。 .. hint:: 什么是观测样本？ 就是从某个概率分布 :math:`p(X;\theta)` 中取得的样本值。 .. _ch_liklihood: 极大似然估计 ################################################### 假设有一个随机变量 :math:`X` ，其概率分布是 :math:`p(X;\theta)` ，其中 :math:`\theta` 是这个概率分布的参数，其值是未知的。 我们有一些这个随机变量的观测值，这些观测值集合用符号 :math:`\mathcal{D}=\{x^{(1)},x^{(2)},\ldots,x^{N}\}` 表示。这些观测值都是从同一个概率分布 :math:`p(X;\theta)` 中得到的样本值， 并且通常我们认为这些样本是独立获取的，即每条样本值不依赖其它样本值， 所以这些样本满足独立同分布(i.i.d)的特性。 我们知道其中任意一条样本 :math:`x^{(i)}` 的发生概率是 :math:`p(x^{(i)};\theta)` ， 那么所有样本发生的联合概率是 :math:`p(\mathcal{D};\theta)=p(x^{(1)},\ldots,x^{(N)};\theta)` ， 又由于所有样本是满足i.i.d的，根据联合概率分布的分解法则有： .. math:: p(\mathcal{D};\theta) = p(x^{(1)},\ldots,x^{(N)};\theta) =\prod_{i=1}^{N} p(x^{(i)};\theta) 现在我们要思考如何得到 :math:`\theta` 的值，我们假设 :math:`\theta` 的可能取值空间为 :math:`\Theta` 。 :math:`\theta` 的取值有很多种可能，要想确定什么值合适就需要有评价的标准， 通过评价标准对比各个不同取值的"优劣"，然后选出其中最优的那个值。 这个评价的标准就是所有观测样本的联合概率 :math:`p(\mathcal{D};\theta)` 。 通常我们认为最有可能发生的是概率最大的那个事件，那么既然这些观测样本已经发生了(我们已经得到)， 我们就认为使得这些样本发生概率最大的那个 :math:`\theta` 值是最优的。 通常我们把观测样本发生的联合概率称为似然(likelihood)，一般用符号 :math:`L(\theta;\mathcal{D})` 表示， 有时也称为似然函数(likelihood function)。 .. math:: L(\theta;\mathcal{D}) = p(\mathcal{D};\theta) = \prod_{i=1}^{N} p(x^{(i)};\theta) 我们说最优的 :math:`\theta` 值是令观测样本发生概率最大的值， 也即是使得似然函数取得最大值时 :math:`\theta` 的取值是最优解。 因此，参数 :math:`\theta` 的估计值 :math:`\hat{\theta}` 为： .. math:: \hat{\theta}_{ML} = \mathop{\arg \max}_{\theta} L(\theta;\mathcal{D}) = \mathop{\arg \max}_{\theta} \prod_{i=1}^{N} p(x^{(i)};\theta) 仔细观察发现似然函数是每条样本概率 :math:`p(x^{(i)};\theta)` 的连乘， 而概率值都是在[0,1]之间的，一系列小于1的数字连乘会趋近于0。 而计算机在处理浮点数时存在精度问题，太小的值是无法表示的。 所以一般我们会为似然函数加上一个对数操作来解决计算机的精度问题， 我们把加了对数的似然函数称为 *对数似然函数(log-likelihood function)* ， 一般用符号 :math:`\ell` 表示。 .. math:: \ell(\theta;\mathcal{D}) = \log L(\theta;\mathcal{D}) 通过极大化对数似然函数 :math:`\ell(\theta;\mathcal{D})` 得到 :math:`\hat{\theta}` 和极大化似然函数 :math:`L(\theta;\mathcal{D})` 是等价的，这里不再证明，有兴趣的读者可以参考其他资料。 .. math:: \hat{\theta}_{ML} &= \mathop{\arg \max}_{\theta} \ell(\theta;\mathcal{D}) &= \mathop{\arg \max}_{\theta} \log \prod_{i=1}^{N} p(x^{(i)};\theta) &= \mathop{\arg \max}_{\theta} \sum_{i=1}^N \log p(x^{(i)};\theta) 那么如何进行极大化求解呢？通常有如下三种方法： 1. 解析法(Analytic)，又叫直接求解法。我们知道一个函数在取得极值时其一阶导数是为0的， 那么通过令对数似然函数的一阶导数为0来解得 :math:`\hat{\theta}_{ML}` 。 通过令对数似然函数的一阶导数为0，然后解等式方程的方法得到的 :math:`\hat{\theta}_{ML}` 称为解析解。 .. math:: \frac{\partial \ell}{\partial \theta} = 0 函数的一阶导数为0的点称为“驻点”(stationary point)，可能为（局部）极大或者极小值点，也可能为鞍点(saddle point)， 可以通过极值点的二阶导数判断是极大值点还是极小值点。 并不是所有情况都能得到解析解的，很多时候是无法直接求得的，在后面的章节中我们会详细讨论。 2. 网格搜索法(Grid Search)。如果我们知道 :math:`\hat{\theta}` 的值在实数域 :math:`R` 的一个子空间中，可以对这个子空间进行搜索 来得到是的似然函数最大的参数值。换句话说，就是尝试这个空间中的每个值，找到令似然函数取得最大的参数值。网格搜索方法是一种很好的方法，它表明可以通过重复逼近和迭代来找到似然函数的最大值。 但是，它在大多数情况下不切实际，并且当参数数量变多时变得更加困难。 3. 数值法(Numerical)。这是现在最常用的算法。本质上就是先为 :math:`\theta` 赋予一个初始值， 然后利用爬山法找到最大值。梯度下降(上升)法(Gradient descent)，牛顿法(Newton-Raphson)，BHHH，DFP等等都属于这类。 本章我们暂时只讨论极大似然估计的直接求解法，在后续的章节中再讨论数值法。 通过极大化似然函数估计参数的方法称为极大似然估计(maximum likelihood estimation,MLE)， 亦可以叫做最大似然估计。下面看几个具体的例子来直观的理解极大似然估计。 二值离散变量 ============================ 我们假设一个离散随机变量 :math:`X` 是伯努利分布(Bernoulli distribution)， 即只有两种可能的取值 :math:`X \in \{0,1\}`， 我们设其取值为1的概率为 :math:`p(X=1)=\theta` ， 其概率分布可以写成: .. math:: p(X;\theta) = \theta^x (1-\theta)^{(1-x)},x \in \{0,1\} 其中 :math:`\theta` 就是我们需要估计出的参数。 假设现在有变量 :math:`X` 的N次独立观测样本，用符号 :math:`\mathcal{D}=\{x^{(1)},\ldots,x^{(N)}\}` 表示， 样本集的规模为 :math:`|\mathcal{D}|=N` 。 我们利用极大似然估计法估计出参数 :math:`\theta` 。 首先我们写出观测样本的对数似然函数。 .. math:: \ell(\theta;\mathcal{D}) &= \sum_{i=1}^N \log p(x^{(i)};\theta) & = \sum_{i=1}^N \log [ \theta^{x^{(i)}} (1-\theta)^{(1-x^{(i)})} ] & = \sum_{i=1}^N \log \theta^{x^{(i)}} + \sum_{i=1}^N \log (1-\theta)^{(1-x^{(i)})} 我们定义几个统计值，在样本集 :math:`\mathcal{D}` 中， :math:`n_0` 表示随机变量 :math:`X=0` 在观测样本中的次数， :math:`\hat{p}_{\mathcal{D}}(0)=\frac{n_0}{N}` 表示随机变量 :math:`X=0` 在观测中出现的相对频次(经验分数)； :math:`n_1` 表示 :math:`X=1` 在样本中的次数， :math:`\hat{p}_{\mathcal{D}}(1)=\frac{n_1}{N}` 表示 :math:`X=1` 在样本中出现的相对频次(经验分数)。 把 :math:`n_0,n_1` 代入到对数似然函数中可得： .. math:: \ell(\theta;\mathcal{D}) &= \sum_{i=1}^N \log \theta^{x^{(i)}} + \sum_{i=1}^N \log (1-\theta)^{(1-x^{(i)})} &= \sum_{i=1}^N x^{(i)} \log \theta + \sum_{i=1}^N (1-x^{(i)}) \log (1-\theta) &= n_1 \log \theta + n_0 \log (1-\theta) 我们知道当对数似然函数的导数为0时，函数取得极大值（注意极大值不是最大值，极大值点也可能不存在，也可能存在多个）。 那么我们只需要对对数似然函数求导得到 :math:`\theta` 的估计值。 .. math:: \begin{aligned} 0 &=\frac{\partial \ell}{\partial \theta} \\ &= \frac{n_1}{\theta}- \frac{n_0}{(1-\theta)} \end{aligned} 又由于 :math:`n_0 = N - n_1` ，带入上式可得： .. math:: \frac{n_1}{N-n_1}=\frac{\theta}{1-\theta} 最后化简可得 :math:`\theta` 的极大似然估计值： .. math:: \hat{\theta}_{M L}=\frac{n_1}{N} = \hat{p}_{\mathcal{D}}(1) 发现没有，最大似然的估计值就等于经验估计值(在样本中出现的相对频次) :math:`\hat{p}_{\mathcal{D}}(1)` 。 这个值称为伯努利分布的充分统计量(sufficient statistics)。 通过强大的大数定律，经验分布最终会收敛到正确的概率，因此我们看到在这种情况下最大似然与之是一致的。 .. hint:: 充分统计量(sufficient statistic) 对于一个概率分布(或者一个分布族)都有一个固定的概率密度(质量)函数，而这个函数中又存在着一些可变参数， 只有当这些参数的值确定时才能确定一个具体的分布。对于一个参数未知的分布， 很多时候我们可以通过这个分布的一些IID样本去估计（比如最大似然估计）出参数的值。 当我们通过分布的样本去估计分布的参数时，自然是需要样本中包含的一些"信息(information)"，利用这些信息去估计参数值。 统计量(statistic)是 **样本的一个函数**，其代表着从样本中提取的一些"信息"，比如样本的均值(mean)，样本的总和(sum)等等。 很多时候这些信息可以用于确定这个分布的未知参数， 如果仅需要一个统计量就能确定这个分布的未知参数，而不再需要其它的额外"信息"，那么这个统计量就称为这个分布(或者分布族) 的 **充分统计量(sufficient statistic)** 。 在估计分布参数时，只需要保留这个充分统计量，而不再需要样本本身。 一般离散变量 ====================================== 当随机变量 :math:`X` 的取值不是只有 :math:`\{0,1\}` 两个取值，而是更一般的离散值空间时 :math:`\mathcal{X}=\{x_1,\dots,x_M\}` ， 我们就用类别分布(Categorical distribution)来表示变量 :math:`X` 的概率分布，也可以叫做单次实验的多项式分布(Multinomial Distribution)， 多项式分布是伯努利分布的扩展。 假设进行N次多项式分布的实验，其中类别m出现的次数记为 :math:`n_m` ，则有 :math:`\sum_{m=1}^M n_m=N` ，如果实验顺序相关，那么多项式分布表示 :math:`n_m` 的联合概率分布。 .. math:: p(n_1,n_2,\ldots,n_M|\theta_1,\theta_2,\ldots,\theta_M) = \frac{\Gamma(N+1)}{\prod_{m=1}^{M} \Gamma (n_m +1 )} \prod_{m=1}^M \theta_m^{n_m} 多项式分布表示的是，在N次实验中，各个类别出现次数为 :math:`(n_1,n_2,\ldots,n_M)` 的概率分布，其考虑了样本的出现顺序， 其中包含Gamma函数的子项是实验顺序相关的组合数， 如果只有单次实验就不需要这部分。 在这里假设我们的观测样本独立的从类别变量 :math:`X` 采样得到的，每一个样本都看成是独立的单次实验。 单次实验的多项式分布的概率质量函数可以表示为： .. math:: :label: 2_100_10 p(X;\theta) = \prod_{m=1}^{M} \theta_m^{\delta (x,x_m)} 其中 :math:`\theta=[\theta_1,\dots,\theta_M]` 为分布函数的参数。 为表示方便，我们定义一个指示函数 :math:`\delta(a,b)` ，当 :math:`a=b` 时，指示函数输出为1，否则输出为0。 .. math:: \delta (x,x_m) = \left \{ \begin{aligned} 1 &, \ x = x_{m} ;\\ 0 &, \ x \ne x_{m} ; \end{aligned} \right . :eq:`2_100_10` 的含义是，变量 :math:`X` 取值为类别 :math:`x_m` 的概率是 :math:`\theta_m` ， 即 :math:`p(x_m)=\theta_m` 。其中参数 :math:`\theta` 需要满足约束： .. math:: \sum_{m=1}^M \theta_m = 1 \theta_m \in [0,1] 观测样本 :math:`\mathcal{D}` 的似然函数为： .. math:: :label: 20_mutil_likelihood L(\theta;\mathcal{D}) &= \prod_{i=1}^N \prod_{m=1}^{M} \theta_m^{\delta (x,x_m)} &=\prod_{m=1}^{M} \theta_m^{n_m} 其中 :math:`n_m` 表示 :math:`x_m` 在样本中出现的次数，那么有 :math:`\sum_{m=1}^M n_m=N` 。 我们为似然函数加上对数操作，以便把连乘符号转换成加法。 .. math:: \ell( \theta;\mathcal{D} ) & \triangleq \log L(\theta;\mathcal{D} ) &= \sum_{m=1}^M {n_m} \log \theta_m 为了找到 :math:`\theta_m` 的最大似然解，我们需要关于 :math:`\theta_m` 最大化对数似然函数 :math:`\ell(\theta;\mathcal{D} )` ， 并且限制 :math:`\theta_m` 的和必须等于1。这样 **带有约束的优化问题需要用到拉格朗日乘数** :math:`\lambda` 实现，即需要最大化如下方程。 .. math:: \sum_{m=1}^M {n_m} \log \theta_m + \lambda ( \sum_{m=1}^M \theta_m -1) 对上述公式求偏导，并令偏导等于0。 .. math:: 0 &=\frac{\partial }{\partial \theta_m} &= \frac{n_m}{\theta_m} + \lambda \theta_m &= - \frac{n_m}{\lambda} 把 :math:`\theta_m = - \frac{n_m}{\lambda}` 代入到约束条件 :math:`\sum_{m=1}^M \theta_m=1` 中， 解得 :math:`\lambda = -N` ，最终可求得 :math:`\theta_m` 的值： .. math:: \theta_m = \frac{n_m}{N} 我们令 :math:`\hat{p}_{\mathcal{D}}(x_m) = \frac{n_m}{N}` ， 类别变量的参数似然估计值和伯努利的例子具有相同的形式。 .. math:: (\hat{\theta}_m)_{ML} = \hat{p}_{\mathcal{D}}(x_m) =\frac{n_m}{N}, 1\le m \le M .. _ch_2_Gaussian_ML: 高斯分布 ============================ 现在我们假设有一个高斯变量 :math:`X \sim N(\mu,\sigma^2)` ，其参数有两个，均值 :math:`\mu` 和方差 :math:`\sigma^2` 。这个变量的一个观测样本集为 :math:`\mathcal{D}=\{x^{(1)},x^{(2)},\ldots,x^{N}\}` ，我们利用极大似然估计出参数 :math:`\mu,\sigma^2` ，我们先写出高斯分布的表达式。 .. math:: p(x;\mu,\sigma^2) = \frac{1}{ (2\pi \sigma^2)^{\frac{1}{2}} } \exp \left \{ - \frac{1}{2\sigma^2} (x-\mu)^2 \right \} 这个概率分布表达式表示一条样本的发生概率，所有样本都发生的概率，也就是似然函数为： .. math:: L(\mu,\sigma^2;\mathcal{D}) = \prod_{i=1}^N p(x^{(i)};\mu,\sigma^2) 其对数似然函数，我们选择以自然对数e为底数的对数ln ，则其对数似然函数为： .. math:: \ell (\mu,\sigma^2;\mathcal{D};) &= ln \prod_{i=1}^N p(x^{(i)};\mu,\sigma^2) &= \sum_{i=1}^N \ln p(x^{(i)};\mu,\sigma^2) &= \sum_{i=1}^N \left [ -\frac{1}{2} \ln 2 \pi \sigma^2 - \frac{1}{2\sigma^2} (x^{(i)}-\mu)^2 \right ] &= \sum_{i=1}^N \left [ -\frac{1}{2} \ln 2 \pi -\frac{1}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} (x^{(i)}-\mu)^2 \right ] &= -\frac{N}{2} \ln 2 \pi -\frac{N}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} \sum_{i=1}^N (x^{(i)}-\mu)^2 然后对参数求偏导数，并令偏导数为0。 .. math:: \begin{cases} \frac{\partial \ell}{\partial \mu} = \frac{1}{\sigma^2} \sum_{i=1}^N (x^{(i)}-\mu) = 0 \frac{\partial \ell}{\partial \sigma^2}= -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4} \sum_{i=1}^N (x^{(i)}-\mu)^2=0 \end{cases} 由第一个等式可以解的： .. math:: \mu = \frac{1}{N} \sum_{i=1}^N x^{(i)} = \bar{x} 其中 :math:`\bar{x}` 表示样本的平均值，然后代入第二个等式解得： .. math:: \sigma^2= \frac{1}{N} \sum_{i=1}^N (x_i-\bar{x}) 我们发现最终似然估计的参数的值只依赖两个量： .. math:: \sum_{i=1}^N x^{(i)},\sum_{i=1}^N (x_i-\bar{x}) 这两个量被称为高斯分布的充分统计量(sufficient statistics)。 总结 ===================================== 通过上面的两个例子，我们来总结下用极大似然方法估计参数的一般模式。 令 :math:`X` 表示随机变量， 令 :math:`p(X;\theta)` 为随机变量 :math:`X` 的概率分布的参数化表示， 其中 :math:`\theta` 是概率分布函数中的未知参数向量 （注意 :math:`\theta` 表示是多个参数的集合，不一定只有一个未知参数）。 我们把 :math:`\theta` 看做是一个 **非随机变量** ，是一系列数值变量，但是其具体取值却是未知的。 - :math:`p(\cdot;\theta)` 是一个参数为 :math:`\theta` 的概率分布的pdf(概率质量函数)或者pmf(概率密度函数) - :math:`L(\cdot ; \mathcal{D}) \triangleq p(\mathcal{D} ; \cdot)` 是随机变量 :math:`X` 的观测数据集 :math:`\mathcal{D}` 的似然函数(likelihood function) 。 我们给出 :math:`N` 个观测样本(训练数据集)，:math:`\mathcal{D}=\{x^{(1)}, \ldots, x^{(N)}\}` ，上标代表样本编号。 目标是学习出参数 :math:`\theta` 。 一个在给定观测时学习未知参数的方法是使用最大似然估计(maximum likelihood estimation,MLE)，也叫极大似然估计。 .. glossary:: 最大似然估计(maximum likelihood estimation,MLE) 最直接的理解是：当使得给定观测数据的似然函数取得最大值时，此时似然函数中未知参数的取值是最优的。 似然函数的概率解释是：这些观测值(观测值是可观测到的随机变量的取值)"同时(不是时间上的同时，是联合发生)"发生的概率。 最大似然就是：既然这些样本事件已经发生了，那么我就假设他们的发生的概率是最大的，就认为使得这些样本具有最大发生概率时参数的值是最优取值。 注意，最大似然估计的解不一定存在，也不一定唯一，这取决于似然函数是否有极值点，以及有几个极值点。 我们假设观察样本都是独立同分布(i.i.d.)于 :math:`p(X;\theta)` ， 则这些观测样本 :math:`\mathcal{D}=\{x^{(1)}, \ldots, x^{(N)}\}` 的对数似然函数是： .. math:: \begin{aligned} \ell (\theta;\mathcal{D}) &= \log L(\theta;\mathcal{D}) &= \log p(x^{(1)}, \ldots, x^{(N)} ;\theta ) &= \log \prod_{i=1}^{N} p(x^{(i)} ; \theta) &= \sum_{i=1}^{N} \log p(x^{(i)} ; \theta) \end{aligned} 然后我们通过极大化对数似然函数的方式求解参数的值。当然有些时候我们可以令偏导数为0直接求得解析解， 然而有时候却无法得到解析解，这时就需要用梯度迭代的方法求解。 .. math:: \hat{\theta}_{ML} = \mathop{\arg \max}_{\theta} \ell (\theta;\mathcal{D}) 极大似然估计有一个很大的缺陷，例如投掷一个普通的硬币3次，每次都是正面朝上，这时极大似然估计正面朝上的概率时结论会是1， 表示所有未来的投掷结果都是正面向上。这明显是有问题的，当数据集较少时非常容易出现错误的结果， 贝叶斯估计可以一定程度上避免这类问题。 .. _ch_Bayesian_estimation: 贝叶斯估计 ################################################### 我们已经讨论了参数的极大似然估计，以及了解到极大似然估计在某些情况下非常容易产生过拟合现象。 现在我们讨论另一种常用的参数估计方法-贝叶斯估计。 在似然估计中，参数 :math:`\theta` 被看做一个数值对象，只有 :math:`X` 是随机变量， 随机变量 :math:`X` 的概率分布是一个以 :math:`\theta` 为参数的概率分布 :math:`p(X;\theta)` 。 变量 :math:`X` 的样本产生是由参数化的边缘概率分布 :math:`p(X;\theta)` 决定。 我们是利用观测样本得到参数的一个估计值(点估计)，然后将参数估计值代入到变量 :math:`X` 的概率函数中 :math:`p(X;\hat{\theta})` ，得到了变量 :math:`X` 的概率分布函数，之后直接用其进行新样本的预测。 .. math:: \mathcal{D} \xrightarrow{\text{MLE}} \hat{\theta}_{\text{数值}} \xrightarrow{\text{代入}} p(X=x_{new};\hat{\theta}_{\text{数值}}) **联合概率** 然而在贝叶斯学派中，参数 :math:`\theta` 也被看做是随机变量，一条样本的产生是由随机变量 :math:`\theta` 和 :math:`X` 共同决定，假设随机变量 :math:`\theta` 的概率分布是 :math:`p(\theta)` ， 随机变量 :math:`X` 是依赖于随机变量 :math:`\theta` 的，所以变量 :math:`X` 的概率分布是以变量 :math:`\theta` 为条件的条件概率分布 :math:`p(X|\theta)` （这里分号变成了竖线，竖线表示条件概率）。 多个随机变量一起组成联合概率分布 :math:`p(X,\theta)` ，这时样本应该是由联合概率分布 :math:`p(X,\theta)` 生成， 根据链式法则，联合概率可以分解成条件概率的乘积。 .. math:: :label: eq_2_10 p(X,\theta) = p(\theta) p(X|\theta) 其中条件概率分布 :math:`p(X|\theta)` 是变量 :math:`X` 的概率分布，是已知的。 边缘概率分布 :math:`p(\theta)` 是变量 :math:`\theta` 的概率分布， :math:`\theta` 是我们的未知参数，显然它的概率分布是未知的。 我们的目标是找到 :math:`X` 的概率分布，然后可以用其预测新样本。此时有两种方法： 1. 找到 :math:`\theta` 的一个估计值 :math:`\hat{\theta}` ，然后就得到条件概率分布 :math:`p(X|\hat{\theta})` ，把条件概率分布 :math:`p(X|\hat{\theta})` 作为 :math:`X` 的概率分布进行后续的预测与分析，似然估计就属于这种。 2. 利用条件概率 :math:`p(X,\theta)` 进行边缘化，得到 :math:`X` 的边缘概率分布 :math:`p(X)` ，把边缘概率分布 :math:`p(X)` 作为变量 :math:`X` 的概率分布。 变量 :math:`X` 的边缘概率分布为： .. math:: :label: eq_2_12 p(X) &= \int p(X,\theta) d\theta &= \int p(\theta) p(X|\theta) d \theta 第二种方法的一个关键问题就是我们需要知道 :math:`\theta` 的边缘概率分布 :math:`p(\theta)` ，不知道 :math:`p(\theta)` 就不能得到 :math:`X` 的边缘概率 :math:`P(X)` 。贝叶斯估计属于第二种方法，其利用贝叶斯定理， 找到 :math:`\theta` 的概率分布 :math:`p(\theta)` 。 **贝叶斯定理** 我们知道变量 :math:`\theta` 和变量 :math:`X` 组成联合概率 :math:`p(\theta,X)` ，并且二者不是相互的独立的， 变量 :math:`\theta` 影响着变量 :math:`X` ，二者存在"因果关系"， :math:`\theta` 是"因"，:math:`X` 是"果"。 联合概率可以通过链式法则分解成一系列条件概率的乘积形式， 链式法则并没有限定变量的顺序(只受到变量间独立性影响)， 所以联合概率 :math:`p(X,\theta)` 有两种分解方式： .. math:: :label: eq_2_13 p(X,\theta) = p(\theta) p(X|\theta) = p(X) p(\theta|X) 通过移项可以得到： .. math:: :label: eq_2_14 p(\theta|X) = \frac{p(\theta) p(X|\theta)}{p(X)} :eq:`eq_2_14` 就是贝叶斯定理，贝叶斯定理的核心就是如下的转换: .. math:: :label: eq_2_15 p(\text{因}|\text{果}) = \frac{p(\text{因})p(\text{果}|\text{因}) }{p(\text{果})} 很多场景下，我们可以看到"果"，也就是我们有变量 :math:`X` 的观测值，但我们不知道导致这个"果"的"因"是什么，也就是不知道变量 :math:`\theta` 是什么。这时我们就可以利用贝叶斯定理推断出"因"，而这就是通常所说的贝叶斯推断(Bayesian inference)， 很多资料中会把"结果"(观测值)称之为证据(evidence)，把"果"变量称为证据变量。 **贝叶斯推断** 变量 :math:`\theta` 是"因"变量，变量 :math:`X` 是"果"变量，而其观测值 :math:`\mathcal{D}` 就是看到的"结果"， 我们把变量 :math:`X` 的观测样本 :math:`\mathcal{D}` 和变量 :math:`\theta` 写成贝叶斯定理的形式： .. math:: :label: eq_2_16 p(\theta|\mathcal{D}) =\frac{p(\mathcal{D}|\theta) p'(\theta)}{p(\mathcal{D})} - :math:`p(\theta|\mathcal{D})` 表示基于"结果" :math:`\mathcal{D}` 推断出的"因"变量 :math:`\theta` 的概率分布， 通常被称为 **后验概率分布** ，这里"后验"就表示有了经验之后，这里的经验就是指"结果"，也就是观测值。 - :math:`p'(\theta)` 表示的是在没有任何证据(观测样本集)时，经验上对 :math:`\theta` 的认知， 称为先验概率分布(prior probability distribution)。 先验一般是根据具体的应用场景凭借经验为之设定一个常见的概率分布， 如果你对 :math:`p(\theta)` 一无所知那可以设定为均匀分布。 注意这里的 :math:`p'(\theta)` 和 :eq:`eq_2_10` 中的 :math:`p(\theta)` 虽然都表示参数变量的边缘概率， 但它们是在贝叶斯估计中不同阶段的表示，所以这里我们加了一个上标"'"进行区分， 后面我们会说明。 - :math:`p(\mathcal{D}|\theta)` 就是在有 :math:`\theta` 的条件下生成观测样本的的概率， 我们知道观测样本集是符合独立同分布(i.i.d)的，所以展开后具有如下形式： .. math:: p(\mathcal{D}|\theta) &= p(\{x^{(1)},\ldots,x^{(N)}\}|\theta) &= \prod_{i=1}^{N} p(x_i|\theta) 我们发现这其实就是样本的似然，所以 :math:`p(\mathcal{D}|\theta)` 就是样本的似然值。 - :math:`p(\mathcal{D})` 是"果"的观测，直观的讲就是观测样本集的概率分布，通常被称为证据(evidence)。 :math:`p(\mathcal{D})` 作为分母，本质上就是归一化因子，是分子所有可能取值的求和，保证输出的 :math:`[0,1]` 区间内合法概率值，可以通过对分子积分(求和)得到。 .. math:: p(\mathcal{D}) = \int p(\mathcal{D}|\theta) p'(\theta) d\theta :math:`p(\mathcal{D})` 作为归一化因子，通过对分子中参数变量积分得到，消除了参数的影响，其不再受到参数的 影响。换句话说，只要样本集 :math:`\mathcal{D}` 确定了，那么 :math:`p(\mathcal{D})` 的值就确定了，不再变化，在确定了样本集后，其是一个固定值。 综上，贝叶斯推断可以表述成如下方式， 其中符号 :math:`\propto` 表示正比关系。 .. math:: :label: eq_2_12 \text{后验概率} &= \frac{\text{似然(likelihood)} \times \text{先验(prior)}}{\text{证据(evidence)}} & \propto \text{似然} \times \text{先验} 我们可以用贝叶斯推断找到参数变量 :math:`\theta` 的后验概率分布 :math:`p(\theta|\mathcal{D})` ，然后把 :math:`p(\theta|\mathcal{D})` 作为参数的"真实"概率分布， 然后代入到 :eq:`eq_2_12` 中， 这样我们就确定了变量 :math:`X` 和变量 :math:`\theta` 的联合概率分布，并且依此得到 :math:`X` 的边缘概率分布： .. math:: p(X) &= \int p(X,\theta) d\theta &= \int p(\theta|\mathcal{D}) p(X|\theta) d \theta 但是要利用 :eq:`eq_2_16` 推断出 :math:`p(\theta|\mathcal{D})` 还存在两个难点： 1. 先验分布 :math:`p'(\theta)` 如何确定。 2. 分母 :math:`p(\mathcal{D})` 需要计算积分，并且是对 :math:`p'(\theta)` 进行积分， :math:`p'(\theta)` 的形式会影响积分的难度。 理论上参数的先验分布应该根据我们其认知信息确定， 但实际上多数情况下我们对参数是一无所知的，没有任何信息， 这时，我们就需要一种被称为无信息先验（noninformative prior）的先验分布。 这种先验分布的目的是尽可能的对后验分布产生小的影响（Jeffries, 1946; Box and Tao, 1973; Bernardo and Smith, 1994）。 这有时也被称为“让数据自己说话”。 除无信息先验外，另外一种确定先验分布的方法为共轭先验(conjugate prior)， 共轭先验是一种使用非常广泛的确定先验分布的方法， 本节我们只讨论共轭先验法。 **共轭先验** .. _ch_conjugate_prior: .. topic:: 共轭先验(conjugate prior) 在贝叶斯推断中，如果后验分布与先验分布属于同一种概率分布，则此先验分布称为共轭先验。 注意，由于后验分布是由先验与似然相乘得到的，所以共轭指的是先验与似然共轭， 共轭先验与似然相乘后，不改变分布的函数形式，所以得到后验与先验具有相同的形式。 共轭先验使得后验分布和先验分布拥有相同的形式， 很多时候可以直接给出后验的结果， 而不必计算分母 :math:`p(\mathcal{D})` ，这极大的降低了后验分布的计算复杂度。 高斯分布的似然函数的共轭分布仍然是高斯分布，伯努利分布的似然函数的共轭先验是beta分布， 类别分布的似然函数的共轭分布是狄利克雷分布， 稍后我们会举例说明。 共轭先验也是有缺点的，**其一是只有指数族分布才存在共轭先验，在下一章我们会详细讨论指数族。** **其二是，选取共轭先验更多是为了计算简单，而不是为了更精确的估计参数。** 选取了合适的参数先验分布后，就可以利用贝叶斯推断 :eq:`eq_2_16` 得到参数的后验概率分布 :math:`p(\theta|\mathcal{D})` ，后验概率分布就是我们在观测样本集的条件下对参数变量 :math:`\theta` 概率分布的估计。 然后就可以用后验概率分布 :math:`p(\theta|\mathcal{D})` 替代 :eq:`eq_2_10` 中参数变量的边缘概率分布。 .. math:: p(\theta) \triangleq p(\theta|\mathcal{D}) :eq:`eq_2_10` 表示的联合概率也就变成： .. math:: p(X,\theta) \triangleq p(\theta|\mathcal{D}) p(X|\theta) 此时变量 :math:`X` 的边缘概率分布为： .. math:: :label: eq_2_12 p(X) &= \int p(X,\theta) d\theta &= \int p(\theta|\mathcal{D}) p(X|\theta) d \theta &= p(X|\mathcal{D}) 有了 :math:`X` 的边缘概率分为，就可以预测新样本的概率： .. math:: p(X=x_{new}) &= p(X=x_{new}|\mathcal{D}) &= \int p(\theta|\mathcal{D}) p(X=x_{new}|\theta) d \theta 伯努利变量 =============================== 我们假设随机变量 :math:`X` 服从伯努利分布(Bernoulli distribution)， 其参数化的概率质量函数为： .. math:: p(X;\theta) = \theta^x (1-\theta)^{(1-x)},x \in \{0,1\} 其中 :math:`\theta` 是需要估计的未知参数， 现在我们认为 :math:`\theta` 也是一个随机变量，并且其概率分布为 :math:`p(\theta)` ， 变量 :math:`X` 的概率质量函数就变成了条件概率分布 :math:`p(X|\theta)` 。 .. math:: p(X|\theta) = \theta^x (1-\theta)^{(1-x)},x \in \{0,1\} **先验分布** 变量 :math:`X` 是伯努利分布，而伯努利分布的似然函数的共轭先验是 *Beta分布* ， 一般可记为 :math:`\theta \sim Beta(\theta|a,b)` 。 .. math:: p(\theta;a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{a-1}(1-\theta)^{b-1} 其中 :math:`\Gamma(\cdot)` 称为Gamma函数，并有如下性质： .. math:: :label: eq_estimate_2 \Gamma(x+1)=x\Gamma(x) \Gamma(n)=(n-1)! \, \text{n是整数} Beta分布的期望和方差为： .. math:: E[\theta] = \frac{a}{a+b} var[\theta] = \frac{ab}{(a+b)^2(a+b+1)} Beta分布是一个连续值的概率分布( :math:`\theta \in [0,1]` 是连续值)，对于一个连续值的概率分布满足积分为1。 .. math:: \int_{0}^1 \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{a-1}(1-\theta)^{b-1} = 1 这个先验分布中有两个参数a,b，一般情况我们会根据经验直接给定这两个参数的值，也就是其值已知的。 那么其中的Gama函数部分 :math:`\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}` 是一个常量。 积分符号内部的常量可以提到积分外面去，我们把这个积分等式做个变形，稍后会用到。 .. math:: :label: eq_estimate_3 \int_{0}^1 \theta^{a-1}(1-\theta)^{b-1} = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} **后验分布** 有了先验分布后，我们把这个先验分布代入到后验分布中，由于a,b的值是确定的，所以先验分布中的Gamma函数部分 :math:`\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}` 是一个常数，与参数 :math:`\theta` 无关，也与观测样本无关。 此外，我们用 :math:`n_1` 表示观测样本中1的次数，用 :math:`n_0` 表示观测样本中0出现的次数， 则有 :math:`n_1+n_0=N` 。 .. math:: :label: eq_estimate_4 p(\theta|\mathcal{D}) &= \frac{L(\theta;\mathcal{D}) p(\theta)}{p(D)} &= \frac{ \left [ \prod \theta^{x^{(i)}} (1-\theta)^{(1-{x^{(i)}})} \right ] \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{a-1}(1-\theta)^{b-1}}{p(D)} &= \frac{ \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \left [ \theta^{n_1} (1-\theta)^{n_0} \right ] \theta^{a-1}(1-\theta)^{b-1}}{p(D)} &= \frac{\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1}}{p(D)} 现在我们来看分母 :math:`p(D)` ，我们知道分母其实是分子的归一化，由于 :math:`\theta` 是连续值，所以分母其实就是分子的积分。 另外借助 :eq:`eq_estimate_2` 和 :eq:`eq_estimate_3` 可以进行简化。 .. math:: p(D) &= \int_0^1 \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} &= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \int_0^1 \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} &= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \frac{\Gamma(a+n_1)\Gamma(b+n_0)}{\Gamma(a+b+n_1+n_0)} 我们把分母代入到后验分布 :eq:`eq_estimate_4` ，可得： .. math:: p(\theta|\mathcal{D}) &= \frac{\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1}}{p(D)} &= \frac{\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1}} {\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \frac{\Gamma(a+n_1)\Gamma(b+n_0)}{\Gamma(a+b+n_1+n_0)}} &= \frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} 发现没有，后验概率分布 :math:`p(\theta|\mathcal{D})` 也是一个Beta分布。 后验与先验具有相同的概率分布形式，这反映出先验关于似然函数共轭的性质。 :math:`\theta` 的先验分布是 :math:`Beta(\theta|a,b)` ， 后验分布是 :math:`Beta(\theta|n_1+a,n_0+b)` ， 而且后验分布就是在先验分布的基础上加上的观测样本的一些统计值。 .. math:: a_{\text{后验}} &= a_{\text{先验}} + n_1 b_{\text{后验}} &= b_{\text{先验}} + n_0 也就是说对于某些概率分布，如果我们选取共轭先验作为参数的先验分布，那么只需要对观测数据集进行一些统计， 就能直接给出参数的后验概率分布。 .. important:: 有了参数 :math:`\theta` 的后验概率分布 :math:`p(\theta|\mathcal{D})` ， 就相当于得到了参数 :math:`\theta` 的"估计值"，和最大似然估计不同的是， 最大似然估计得到的是点估计(参数一个数值估计)。 而贝叶斯估计是把参数看做一个随机变量，得到的是参数的后验概率分布，类似于区间估计。 **预测新样本** 有了参数的估计（后验概率分布）后，就相当于确定了变量 :math:`X` 和 :math:`\theta` 的联合概率分布 :math:`p(X,\theta)=p(\theta|\mathcal{D}) p(X|\theta)` ， 通过对联合概率的边缘化得到变量 :math:`X` 的边缘概率分布： .. math:: p(X|\mathcal{D}) = \int_0^1 p(\theta|\mathcal{D}) p(X|\theta) d\theta 利用 :math:`p(X)` 我们可以预测新样本的概率: .. math:: p(x_{new}|\mathcal{D}) &= \int_0^1 p(x|\theta) p(\theta|\mathcal{D}) d \theta &= \int_0^1 \theta^x (1-\theta)^{(1-x)} \frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} &=\frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} \int_0^1 \theta^x (1-\theta)^{(1-x)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} &=\frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} \int_0^1 \theta^{n_1+a-1+x} (1-\theta)^{n_0+b-1+1-x} &= \frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} \frac{\Gamma(a+n_1+x)\Gamma(b+n_0+1-x)}{\Gamma(a+b+n_1+n_0+1)} &= \frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} \frac{\Gamma(a+n_1+x)\Gamma(b+n_0+1-x)}{ (a+b+n_1+n_0) \Gamma(a+b+n_1+n_0)} &= \frac{\Gamma(a+n_1+x)\Gamma(b+n_0+1-x)}{ (a+b+n_1+n_0){\Gamma(a+n_1)\Gamma(b+n_0)} } :math:`x_{new}=1` 的概率为： .. math:: p(x_{new}=1|\mathcal{D}) &=\frac{\Gamma(a+n_1+1)\Gamma(b+n_0+1-1)}{ (a+b+n_1+n_0){\Gamma(a+n_1)\Gamma(b+n_0)} } &=\frac{\Gamma(a+n_1+1)}{ (a+b+n_1+n_0){\Gamma(a+n_1)} } &=\frac{(a+n_1)\Gamma(a+n_1)}{ (a+b+n_1+n_0){\Gamma(a+n_1)} } &=\frac{a+n_1}{ a+b+n_1+n_0 } &= \frac{a+n_1}{ a+b+N } 我们发现用积分法去计算新样本的概率太复杂，实际上这个过程是可以简化的。 我们知道对概率分布求积分相当于其期望，所以我们可以求出参数的后验概率分布的期望值， 然后把期望值作为参数的一个点估计值，代入到变量 :math:`X` 的条件概率中，通常称为 **后验期望法(mean of the posterior)** 。 .. math:: \hat{\theta}_{Bayes} = \int \theta p(\theta|\mathcal{D}) d \theta = \mathbb{E}_{p(\theta|\mathcal{D})} [\theta] Beta分布的期望值我们可以直接写出： .. math:: \hat{\theta}_{Bayes} = \mathbb{E}_{p(\theta|\mathcal{D})} [\theta] = \frac{a+n_1}{a+b+n_1+n_0} 把这个估计值直接带入到 :math:`X` 的条件概率 :math:`p(X|\theta)` 分布中，同样可以预测下一个样本的值。 .. math:: p(x_{new}=1|\hat{\theta}_{Bayes} ) = \hat{\theta}_{Bayes} (1-\hat{\theta}_{Bayes})^{(1-1)} = \hat{\theta}_{Bayes} = \frac{a+n_1}{a+b+n_1+n_0} = \frac{a+n_1}{a+b+N} 我们发现这和上面通过积分法进行预测是等价，实际上 :math:`\int_0^1 p(x|\theta) p(\theta|\mathcal{D}) d \theta` 就相当于在求期望。 最后，回顾一下伯努利分布的极大似然估计的结果 :math:`\hat{\theta}_{ML}=\frac{n_1}{N}` ， 和贝叶斯估计的结果对比一下，发现贝叶斯估计的结果就是在极大似然估计的基础上加入了先验知识， *加入先验类似于似然函数加上惩罚项，可以起到防止过拟合的效果。* 类别变量 ====================================== 上一节我们以二值离散变量(伯努利变量)为例讨论了贝叶斯参数估计的方法，本节我们扩展到更一般的离散变量--类别变量。 假设随机变量 :math:`X` 一个多值离散变量，可以称之服从类别分布，或者叫多项式分布， 其参数向量用符号 :math:`\theta` 表示，这里 :math:`\theta` 不再是一个标量，而是一个向量。 在最大似然参数估计中，我们将参数 :math:`\theta` 建模为 *非随机变量(nonradom)* 。 相反的，在贝叶斯参数估计(Bayesian parameter estimation)中， 我们把未知参数 :math:`\theta` 看成是 *随机变量* ，并且我们认为边缘分布是 :math:`p(\theta)` ， 类别变量 :math:`X` 的概率分布为条件概率分布 :math:`p(X|\theta)` 。 注意，我们用符号 :math:`p(X|\theta)` 取代了符号 :math:`p(X;\theta)` ， 因为 :math:`\theta` 是一个随机变量。 随机变量 :math:`X` 和 :math:`\theta` 组成联合概率分布 :math:`p(X,\theta)` ，这时边缘概率分布 :math:`p(X)` 需要通过边际化的方法得到， 观测变量 :math:`X` 的边缘概率分布通过下式给出： .. math:: :label: fm_20_0010 p(X) = \int p(X,\theta) d \theta = \int p(\theta) p(X|\theta) d \theta 在极大似然估计中，我们尝试寻找一个最优的 :math:`\theta` 值， 而在贝叶斯估计中，我们并不是求解出参数变量 :math:`\theta` 具体的取值。 而是利用训练数据计算出参数变量 :math:`\theta` 的后验概率分布(posterior distribution)。 我们继续用符号 :math:`\mathcal{D}=\{x^{(1)},x^{(2)},\dots,x^{(N)}\}` 表示观测数据集， 符号 :math:`x^{(1)}` 的上标数字代表观测数据集(训练数据集)的样本编号， 数据集的大小是 :math:`|\mathcal{D}|=N` ， 参数变量 :math:`\theta` 的后验概率分布可以表示为： .. math:: p(\theta|\mathcal{D})=p(\theta|x^{(1)},x^{(2)},\dots,x^{(N)}) 后验概率分布可以根据贝叶斯定理求得： .. math:: :label: 20_bayes p(\theta|\mathcal{D}) =\frac{p'(\theta) p(\mathcal{D}|\theta)}{p(\mathcal{D})} 我们先看分子部分，其中 :math:`p(\mathcal{D}|\theta)` 就是观测数据集的似然函数， :math:`p'(\theta)` 是参数变量 :math:`\theta` 的先验概率分布。 先验分布的选择，会影响着计算后验分布的复杂程度，一个不合适的先验分布，甚至导致后验分布无法进行积分计算。 选取共轭先验可以令后验分布拥有与先验相同的形式，这使得计算得以简化。 上文我们知道，类别分布的似然函数具有多项式分布的形式，而多项分布的共轭先验是狄利克雷(Dirichlet)分布。 分母部分其实就是对分子的归一化，在 :math:`\mathcal{D}` 已经确定的情况下是个常数。 .. hint:: 后验概率分布正比于似然函数和先验的乘积，如果我们选择一个拥有合适函数形式的先验分布， 使得两者乘积的结果(后验概率分布)拥有和先验分布相同的函数形式，这样的先验分布就叫做共轭先验，这个性质称为共轭性(conjugacy)。 比如对于具有多项式分布形式的似然函数是 :math:`\theta_m` 的幂指函数 :eq:`20_mutil_likelihood` ， 如果选择一个同样是关于 :math:`\theta_m` 的幂指函数的先验分布，则二者乘积后的后验概率分布也就具有同样函数形式。 这里我们假设随机变量 :math:`X` 是一个以 :math:`\theta` 为参数变量的类别分布， 其概率质量函数(pmf)如下式。注意，这里 :math:`\theta` 不是一个数值参数， 而是一个随机变量，所以这个公式是个条件概率分布。 .. math:: p(x|\theta) = \prod_{m=1}^{M} \theta_m^{\delta (x,x_m) } 其中 :math:`\delta (x,x_m)` 是一个指示函数，当 :math:`x=x_m` 时， :math:`\delta (x,x_m)=1` ；反之， :math:`\delta (x,x_m)=0` 。 我们令 :math:`n_m` 表示 :math:`x_m` 的样本在全部数据集中出现的次数， 则似然函数(未加对数log)为： .. math:: :label: 20_liklihood L(\mathcal{D};\theta) =p(\mathcal{D}|\theta) &= \prod_{i=1}^N \prod_{m=1}^M \theta_m^{\delta (x,x_m) } &=\prod_{m=1}^M \theta_m^{n_m} **先验分布** 我们看到，类别分布的似然函数是一个多项式分布的形式，多项式分布的共轭先验是狄利克雷分布， 所以这里我们选取狄利克雷分布作为参数变量 :math:`\theta` 的先验分布。 注意先验的参数 :math:`\alpha` 的值我们根据经验直接给出即可，这里看做是已知量。 狄利克雷分布的概率函数为： .. math:: :label: 20_prior p'(\theta;\alpha) &= \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{\alpha_m -1 } & \propto \prod_{m=1}^M \theta_m^{\alpha_m - 1} .. hint:: 狄利克雷分布是一个 *多元连续变量* 的分布，一个概率分布同时输出多个子变量 :math:`\theta_m(1\le m \le M)` 的概率值， 并满足约束 :math:`\sum_m \theta_m = 1` 。 如果不是很理解狄利克雷分布，请看考 :numref:`补充狄利克雷分布的章节引用` 狄利克雷分布是连续值分布，所以满足积分为1的约束。 .. math:: \int \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{\alpha_m -1 } d \theta= \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \int \prod_{m=1}^M \theta_m^{\alpha_m -1 } d \theta = 1 我们把这个积分式稍微变换一下，稍后会使用到。 .. math:: :label: 20_integrated_change \int \prod_{m=1}^M \theta_m^{\alpha_m -1 } d \theta = \frac{\prod_m \Gamma(\alpha_m)}{\Gamma(\sum_m \alpha_m)} **后验分布** 我们把先验分布 :eq:`20_prior` 、似然函数 :eq:`20_liklihood` 代入到求后验分布的贝叶斯公式 :eq:`20_bayes` 中： .. math:: :label: 20_posterior_distribution p(\theta|\mathcal{D}) &= \frac{p'(\theta;\alpha)p(\mathcal{D}|\theta)}{p(\mathcal{D})} &= \frac{ p'(\theta;\alpha) L(\mathcal{D};\theta)}{p(\mathcal{D})} &= \frac{ \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{\alpha_m -1 } \prod_{m=1}^M \theta_m^{n_m} } {p(\mathcal{D})} &= \frac{ \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } }{p(\mathcal{D})} 我们再看分母部分 :math:`p(\mathcal{D})` ，分母是对分子的归一化， 由于这里 :math:`\theta` 是连续值变量，所以分母是对分子的积分。 也可以理解成是对联合概率分布 :math:`p(\mathcal{D},\theta;\alpha)` 进行边际化求得边缘概率 :math:`p(\mathcal{D})` 。 .. math:: p(\mathcal{D}) &= \int p(\mathcal{D},\theta;\alpha) d \theta &= \int p(\theta;\alpha)p(\mathcal{D}|\theta) d \theta &= \int \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{\alpha_m -1 } \prod_{m=1}^M \theta_m^{n_m } d \theta &= \int \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } d \theta &= \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \int \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } d \theta 参考一下积分变换 :eq:`20_integrated_change` ，其中的积分部分可以改写一下得到 :math:`p(\mathcal{D})` 。 .. math:: p(\mathcal{D}) = \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \frac{\prod_m \Gamma(n_m+ \alpha_m)}{\Gamma(\sum_m n_m + \alpha_m)} 我们把这个代入回后验概率分布 :eq:`20_posterior_distribution` 的分母部分。 .. math:: p(\theta|\mathcal{D}) &= \frac{ \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } } {p(\mathcal{D})} &= \frac{ \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } } { \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \frac{\prod_m \Gamma(n_m+ \alpha_m)}{\Gamma(\sum_m n_m + \alpha_m)} } &= \frac{\Gamma(\sum_m n_m + \alpha_m)}{\prod_m \Gamma(n_m+ \alpha_m)} \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } 我们看到后验概率分布仍然是一个狄利克雷分布(PS:是不是很神奇！)， 多项式分布的参数进行贝叶斯估计时，参数的共轭先验是狄利克雷分布，得到的参数后验概率分布也是狄利克雷分布 .. math:: \theta|\mathcal{D} \sim Dirichlet(\alpha_1+n_1,\dots, \alpha_M+n_m) **预测新样本** 有了参数的后验概率分布后，我们就可以预测多项式随机变量 :math:`X` 新的样本概率， 在推导预测概率的过程中需要利用几个性质。 - 上文的积分变换 :eq:`20_integrated_change`， - Gamma函数的一个性质 :math:`\Gamma(x+1)=x\Gamma(x)` 。 - :math:`1= \sum_m \mathbb I \{x_m=1\}` ，因为在类别变量 :math:`\mathbf{x}` 的一个采样中 :math:`[x_1,x_2,\ldots,x_M]` 只有一个元素为1，其它都是0 。 - :math:`N= \sum_m n_m ` ，:math:`n_m` 表示类别 :math:`x_m` 在观测样本集 :math:`\mathcal{D}` 中出现的次，所有类别出现次数求和自然是总样本数N。 - 还有一个需要好好理解的，注意指示函数 :math:`\mathbb I \{x_m=1\}` 表示在类别变量 :math:`\mathbf{x}` 的一个采样中 :math:`[x_1,x_2,\ldots,x_M]` :math:`x_m` 为1，其它都是0 。 .. math:: \prod_m \Gamma(n_m+ \alpha_m+\mathbb I \{x_m=1\} ) &= \Gamma(n_m+ \alpha_m + 1 ) \prod_{\mathbb{m} \in {M \text{去掉} m}} \Gamma(n_\mathbb{m}+ \alpha_\mathbb{m}) &= (n_m+ \alpha_m) \Gamma(n_m+ \alpha_m) \prod_{\mathbb{m} \in {M \text{去掉} m}} \Gamma(n_\mathbb{m}+ \alpha_\mathbb{m}) &=(n_m+ \alpha_m) \prod_m \Gamma(n_m+ \alpha_m) 在观测集 :math:`\mathcal{D}` 的条件预测一个新样本 :math:`x_{new}` 的预测分布为： .. math:: p(x_{new} = x_m|\mathcal{D} ) &= \int \underbrace{p(\theta|\mathcal{D})}_{\text{后验概率分布}} \underbrace{p(x^{'} | \pmb{\mathrm{\theta}})}_{\text{类别分布}} d \theta & = \frac{\Gamma(\sum_m n_m + \alpha_m)}{\prod_m \Gamma(n_m+ \alpha_m)} \int \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } \prod_{m=1}^{M} \theta_m^{\mathbb I \{x_m=1\} } d \theta &=\frac{\Gamma(\sum_m n_m + \alpha_m)}{\prod_m \Gamma(n_m+ \alpha_m)} \int \prod_{m=1}^M \theta_m^{n_m + \alpha_m - 1 } \prod_{m=1}^{M} \theta_m^{\mathbb I \{x_m=1\} } d \theta &=\frac{\Gamma(\sum_m n_m + \alpha_m)}{\prod_m \Gamma(n_m+ \alpha_m)} \int \prod_{m=1}^M \theta_m^{n_m + \alpha_m + \mathbb I \{x_m=1\} - 1 } d \theta &= \frac{\Gamma(\sum_m n_m + \alpha_m)}{\prod_m \Gamma(n_m+ \alpha_m)} \frac{\prod_m \Gamma(n_m+ \alpha_m+\mathbb I \{x_m=1\} )}{\Gamma(\sum_m (n_m + \alpha_m+\mathbb I \{x_m=1\}))} &= \frac{\Gamma(\sum_m n_m + \alpha_m)}{\prod_m \Gamma(n_m+ \alpha_m)} \frac{(n_m+ \alpha_m) \prod_m \Gamma(n_m+ \alpha_m) }{\Gamma( \sum_m (n_m + \alpha_m) + \sum_m \mathbb I \{x_m=1\} ))} &= \frac{(n_m+\alpha_m) \Gamma(\sum_m n_m + \alpha_m)}{\Gamma(\sum_m (n_m + \alpha_m)+ 1 ))} &= \frac{(n_m+\alpha_m) \Gamma(\sum_m n_m + \alpha_m) }{\sum_m (n_m + \alpha_m) \Gamma(\sum_m n_m + \alpha_m)} &= \frac{(n_m+\alpha_m) }{\sum_m (n_m + \alpha_m) } &= \frac{\alpha_m +n_m } {N + \sum_{m=1}^M \alpha_m } 同理，我们可以先计算出后验分布 :math:`p(\theta|\mathcal{D})` 的期望值 :math:`\hat{\theta}_{E}`， 然后将期望值带入到条件概率 :math:`p(X|\theta)` 。 .. math:: \hat{\theta}_{m} = \mathbb{E}[p(\theta_{m}|\mathcal{D}) ] .. math:: p(x_{new} = x_m|\mathcal{D} ) = p(x_{new} = x_m |\theta=\hat{\theta}_{m}) = \hat{\theta}_{m} 贝叶斯估计计算后验概率分布的过程是困难的，需要在整个参数空间求和或者求积分，这在通常情况下是非常困难的(采用共轭先验会简化)， 然后在做预测或者模型比较时又要再次积分(求期望需要积分)。 此外，当数据集规模较小时，后验分布的结果接近先验的结果，当数据集足够大时，贝叶斯估计的结果就会逐渐偏离先验，等价于极大似然估计的结果。 当数据集规模趋近于无穷时，贝叶斯估计的结果和极大似然的结果是一致的。 **在实际应用中，贝叶斯估计先验的选择通常是为了计算方便(共轭先验)而不是为了反映出任何真实的先验知识，** **然而当先验选择不好的时候，贝叶斯方法有很大可能得到错误的结果。** 最大后验估计 ##################################### 贝爷估计有个很大的难点就是计算 :math:`p(\mathcal{D})` ，计算 :math:`p(\mathcal{D})` 需要对参数空间进行积分， 而积分操作的成本很多时候是非常高昂的，甚至都无法计算。 如果我们仅仅是为了预测 :math:`X` 的新样本，而不需要对参数变量本身进行过多的探索， 那么我们不需要得到完整的后验分布，而是只得到参数的一个点估计即可，类似于似然估计。 贝叶斯估计是通过计算参数后验概率分布的期望值作为参数的估计值， 而计算期望值就需要完整的计算出概率分布。 然而除了期望值，我们还可以用后验概率分布最大概率的参数值作为参数的估计值， 称之为最大后验估计(Maximum a posteriori estimation,MAP)。 .. math:: \hat{\theta}_{MAP} = \mathop{\arg \max}_{\theta} p(\theta|\mathcal{D}) 回顾一下公式 :eq:`eq_2_12` ，后验概率的分母 :math:`p(\mathcal{D})` 是一个定值， 后验概率是正比于先验乘以似然的。 .. math:: \text{后验概率} &= \frac{\text{似然} \times \text{先验}}{evidence} & \propto \text{似然} \times \text{先验} 在进行最大后验估计时，其实并不需要先计算出后验概率分布 :math:`p(\theta|\mathcal{D})` 的具体形式， 因为后验概率是正比于分子的 :math:`p(\theta|\mathcal{D}) \propto L(\theta;\mathcal{D}) p(\theta)` ，所以极大化求解后验概率分布时，只需要极大化后验分布的分子即可。 .. math:: \hat{\theta}_{MAP} &= \mathop{\arg \max}_{\theta} p(\theta|\mathcal{D}) &\triangleq \mathop{\arg \max}_{\theta} \text{似然} \times \text{先验} &= \mathop{\arg \max}_{\theta} L(\theta;\mathcal{D}) p(\theta) &\triangleq \mathop{\arg \max}_{\theta} \log L(\theta;\mathcal{D}) p(\theta) &= \mathop{\arg \max}_{\theta} \underbrace{\log L(\theta;\mathcal{D})}_{\text{对数似然}} + \underbrace{\log p(\theta)}_{\text{对数先验}} 发现没有， **最大后验估计就是在极大似然估计的基础上多了一个参数的先验！！！** 所以最大后验估计很多方面是和极大似然估计类似的，但由于多了先验和极大似然估计又有些不同。 .. hint:: 其实最大后验估计加的先验和损失函数加正则项是等价的。对参数引入拉普拉斯先验等价于L1正则化，高斯先验相当于L2正则。 PS:如果你不知道什么是损失函数、正则项，没关系可以暂时无视这句话，以后就会懂的。 伯努利变量 ======================================== 我们截取 :eq:`eq_estimate_4` 的分子部分 .. math:: p(\theta|\mathcal{D}) \propto \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} 最大后验估计的结果为： .. math:: \hat{\theta}_{MAP} &= \mathop{\arg \max}_{\theta} p(\theta|\mathcal{D}) &\triangleq \mathop{\arg \max}_{\theta} \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} 我们通过目标函数加对数，并且令导数为0的方法进行求解。 .. math:: \frac{ \partial}{\partial \theta} \log J(\theta) &= \frac{ \partial}{\partial \theta} \log \frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} \theta^{n_1+a-1} (1-\theta)^{n_0+b-1} &= \frac{ \partial}{\partial \theta} \log \frac{\Gamma(a+b+n_1+n_0)}{\Gamma(a+n_1)\Gamma(b+n_0)} + \frac{ \partial}{\partial \theta} \log \theta^{n_1+a-1} + \frac{ \partial}{\partial \theta} \log (1-\theta)^{n_0+b-1} &= \frac{ \partial}{\partial \theta} (n_1+a-1) \log \theta + \frac{ \partial}{\partial \theta} (n_0+b-1) \log (1-\theta) &= \frac{(n_1+a-1)}{\theta} - \frac{ (n_0+b-1)}{1-\theta} &= 0 解得： .. math:: \hat{\theta}_{MAP} = \frac{n_1+a-1}{n_1+n_0+a+b-2} = \frac{n_1+a-1}{N+a+b-2} 类别变量 ========================================= 同理，我们截取 :eq:`20_posterior_distribution` 的分子部分作为极大化的目标函数。 .. math:: \hat{\theta}_{MAP} = \mathop{\arg \max}_{\theta} \frac{\Gamma(\sum_m \alpha_m)}{\prod_m \Gamma(\alpha_m)} \prod_{m=1}^M \theta_m^{n_m+\alpha_m -1 } 同理加对数，求偏导，然后令偏导数为0，可解得： .. math:: \hat{\theta}_{m} = \frac{\alpha_m +n_m - 1}{\sum_m \alpha_m +N - M} 与最大似然估计的结果进行对比， MAP估计结果考虑了训练样本的数量。特别是，当N很小时，MAP值接近先验的结果;当N很大时，MAP值接近经验分布。 从这个意义上讲，MAP估计可以看作是通过惩罚更复杂的模型来控制过度拟合，即那些离先验更远的模型。 最大似然估计与贝叶斯估计的对比 ##################################### 在概率论一直存在着两者学派，一个是频率学派，一个是贝叶斯学派。 这里我们不讨论这两个学派的本质差别，只关注它们在参数空间上的差别。 通常我们用概率分布(probability distribution)去描述一个随机变量， 我们会说一个随机变量会服从于什么概率分布，比如一个随机变量 :math:`X` 服从于伯努利分布。而一个概率分布都包含一个或多个参数，只有当参数的值确定时才能唯一确定一个分布。 当一个概率分布的参数值未知时，我们需要找到它来确定这个概率分布，然后利用这个概率分布去做一些有价值的事情。 频率学派和贝叶斯学派在参数的认知上存在着很大差异。 频率学派认为概率分布中的参数值就仅仅是一个数值，所以用参数化的方法定义概率分布 :math:`p(X;\theta)` ，并且他们认为参数空间中只有一个值是最优的(或者说是真理)，需要做的就是想办法找到它。 因此在这个基础上提出了最大似然估计法，目标是找出那个最优的参数值。 当然要想估计出参数 :math:`\theta` 的值，我们需要有随机变量 :math:`X` 的一些观测样本， 我们通过这些样本去估计这个概率分布的未知参数。这些样本都是同一个概率分布 :math:`p(X;\theta)` 的样本，所以它们是同分布的， 而且样本与样本之间通常没有什么关系，所以观测样本集都是满足IID(独立同步分布)的。 我们用符号 :math:`\mathcal{D}=\{x_1,x_2,\cdots,x_N \}` 表示这个样本集， 其中每一条样本的发生概率是 :math:`p(x_i;\theta)` ，那么所有样本都发生的概率是一个联合概率： .. math:: p(\mathcal{D};\theta) = p(x_1,x_2,\cdots,x_N ;\theta) = \prod_{i=1}^N p(x_i;\theta) :math:`p(\mathcal{D};\theta)` 通常被称为 *似然函数(likelihood function)* ，习惯上我们用符号 :math:`L(\theta;\mathcal{D})` 表示似然函数。 最大似然估计的思想就是：使得这个样本集发生的联合概率(似然函数)最大的那个参数值是最优的。所以最大似然的估计值为： .. math:: \hat{\theta}_{ML} = \mathop{\arg \max}_{\theta} p(\mathcal{D};\theta) = \mathop{\arg \max}_{\theta} L(\theta;\mathcal{D}) 有了参数 :math:`\theta` 的估计值，就确定了变量 :math:`X` 的概率分布 :math:`p(X;\hat{\theta}_{ML})` ，然后就可以预测新的样本。 .. math:: p(X=x_{new}) = p(X=x_{new};\hat{\theta}_{ML}) 然而贝叶斯学派的观点却恰恰相反，他们认为未知参数既然是未知，那么这个参数取值为参数空间中任意一个值都是有可能的， 所以参数本身也是一个随机变量，也需要用一个概率分布去描述(贝叶斯派的核心一切未知的变量都是随机变量)， 因此他们把带参数的概率分布定义成一个 **条件概率** :math:`p(X|\theta)` （注意这里和频率派有了本质差别）。 同时，他们利用贝叶斯定理把随机变量 :math:`X` 和参数变量 :math:`\theta` 两者之间的关系变成"可逆"的。 .. math:: p(\theta|X) = \frac{p(X|\theta)p(\theta)}{p(X)} 通过贝叶斯定理我们把变量 :math:`X` 和参数变量 :math:`\theta` 的关系定义出来了， 公式中的 :math:`p(\theta)` 表示参数变量 :math:`\theta` 的边缘概率分布， 是在随机变量 :math:`X` 发生之前 :math:`\theta` 的概率分布， 所以我们称之为 :math:`\theta` 的先验分布(prior distribution)。 但实际上我们并不知道参数 :math:`\theta` 的真实概率分布是什么， 所以通常我们会为其假设一个概率分布。 我们假设 :math:`\theta` 的先验概率分布为某一个已知的分布， 然后在这个先验分布 :math:`p(\theta)` 以及条件概率分布 :math:`p(X|\theta)` 情况下， 观测到了变量 :math:`X` 的一些样本 :math:`\mathcal{D}=\{x_1,x_2,\cdots,x_N \}` ， 这个样本集中的所有样本都是从联合概率 :math:`p(X,\theta)=p(X|\theta)p(\theta)` 中采样得到的， 现在我们希望能从这个样本集中反推出 :math:`\theta` 的真实概率分布。 也就是在观测样本集的条件下 :math:`\theta` 的概率分布，这些样本都是随机变量 :math:`X` 的采样， 可以把每个样本点都看成随机变量 :math:`X` 的一个副本， 所以有 :math:`p(\theta|X)\Rightarrow p(\theta|x_1,x_2,\cdots,x_N)= p(\theta|\mathcal{D})` 。 .. math:: :label: eq_2_20 p(\theta|\mathcal{D}) &= \frac{p(\mathcal{D}|\theta)p(\theta)}{p(\mathcal{D})} &=\frac{p(\theta) \prod_{i=1}^N p(x_i|\theta)}{p(x_1,x_2,\cdots,x_N)} 条件概率分布 :math:`p(\theta|\mathcal{D})` 称为参数 :math:`\theta` 的 *后验概率分布(posterior distribution)* ， 因为是在观测样本的条件下 :math:`\theta` 的概率分布，所以称为后验。 后验概率分布是我们在样本集的基础上对参数 :math:`\theta` 的估计分布， 我们把后验概率分布 :math:`p(\theta|\mathcal{D})` 作为参数 :math:`\theta` 的估计。 有了参数的估计分布后，我们就可以用来预测新的样本。在贝叶斯理论的前提下，随机变量 :math:`X` 的样本是由 联合概率 :math:`p(X,\theta)=p(\theta)p(X|\theta)` 产生的，其中 :math:`\theta` 的概率分布我们用估计的后验概率分布替换，所以新的样本的预测分布为： .. math:: p(X=x_{new}) = \int p(\theta|D)p(X=x_{new}|\theta) d \theta 这个方式其实等价于把 :math:`\theta` 的期望值 :math:`\hat{\theta}_{E}` 作为估计值，然后把估计值代入条件概率 :math:`p(X|\theta)` 进行预测。 .. math:: \hat{\theta}_{E} = \mathbb{E}_{p(\theta|\mathcal{D})}[\theta] =\int \theta p(\theta|\mathcal{D}) d \theta p(X=x_{new}) = p(X=x_{new}|\theta=\hat{\theta}_{E}) 后验概率的期望值通常被称为参数的贝叶斯估计(Bayes estimate)： .. math:: \hat{\theta}_{Bayes} = \mathbb{E}_{p(\theta|\mathcal{D})}[\theta] 然而，并不是所有情况下都能求出后验概率分布的期望值的， 要想求得后验概率分布的期望值，就需要求出后验概率分布 :math:`p(\theta|\mathcal{D})` 的具体形式， 后验概率分布 :eq:`eq_2_20` 中的分母是对分子的积分，很多时候这个积分的计算复杂度是很高的，以至于无法计算出来。 .. math:: p(\mathcal{D}) = \int p(\mathcal{D}|\theta)p(\theta) d\theta 因此有时候我们是无法得到后验概率分布的期望的。而且就算我们得到了后验概率分布的具体形式， 要计算后验概率分布的期望有需要对后验概率分布进行积分，这在很多时候也是无法达成的。 所以贝叶斯推断中还有另外一种参数估计方法，*最大后验估计(maximum a posterior)* ： .. math:: \hat{\theta}_{MAP} = \mathop{\arg \max}_{\theta} p(\theta|\mathcal{D}) 最大后验估计的思想是令后验概率中概率最大的那个值作为参数的估计值，而不是期望值。 我们发现后验概率 :eq:`eq_2_20` 是正比于分子部分的。 .. math:: :label: eq_2_21 p(\theta|\mathcal{D}) \propto p(\mathcal{D}|\theta)p(\theta) 我们只需要通过极大化分子就能得到 :math:`\theta` 的最大后验估计值 :math:`\hat{\theta}_{MAP}` ， 所以我们不需要计算积分。 .. math:: \hat{\theta}_{MAP} &= \mathop{\arg \max}_{\theta} p(\theta|\mathcal{D}) &\triangleq \mathop{\arg \max}_{\theta} p(\mathcal{D}|\theta)p(\theta) 我们用 :math:`\hat{\theta}_{MAP}` 作为参数的一个具体估计值，然后用于预测新的样本。 .. math:: p(X=x_{new}) = p(X=x_{new}|\theta=\hat{\theta}_{MAP}) 此外，我们发现其中的 :math:`p(\mathcal{D}|\theta)` 和似然函数 :math:`L(\theta;\mathcal{D})=p(\mathcal{D}|\theta)` 是等价的， :eq:`eq_2_21` 可以表示成： .. math:: \text{后验概率(posterior)} \propto \text{似然(likelihood)} \times \text{先验(prior)} 最大后验估计相当于一个带惩罚(约束)的最大似然估计。 .. math:: \hat{\theta}_{MAP} &= \mathop{\arg \max}_{\theta} p(\mathcal{D}|\theta)p(\theta) &= \mathop{\arg \max}_{\theta} \{ \log p(\mathcal{D}|\theta) + \log p(\theta) \} 统计量和充分统计量 ######################################################## 在本节中，我们讨论充分性的重要概念。 假设我们有一个随机样本集 :math:`\mathcal{D} = \{x_1,x_2,\dots,x_N\}` ，样本集中的样本都是从同一个概率分布 :math:`p(x|\theta)` 采样得到，其中 :math:`\theta` 是这个分布的未知参数，参数空间为 :math:`\Theta` 。 我们可以从这个样本集中估计出未知参数 :math:`\theta` 的值，常见的参数估计方法有：矩估计(method of moments)， 最大似然估计(maximum likelihood)，贝叶斯估计(Bayes estimation)。 我们先回顾一下极大似然估计，极大似然估计是通过极大化似然函数求得参数的估计值，似然函数其实就是样本的联合概率。 比如我们有一个概率分布 :math:`p(x|\theta)` ，其中 :math:`\theta` 是这个分布的未知参数， 这个分布的i.i.d观测样本为 :math:`\mathcal{D}=\{x_1,x_2,\dots,x_N\}` ， 所有样本都发生的联合概率为 :math:`p(x_1,x_2,\dots,x_N|\theta)=\prod_{i=1}^N p(x_i|\theta)` ， 这个联合概率就是似然函数，通过极大化求得参数的估计值 :math:`\hat{\theta}` 。 从参数估计章节的例子中我们可以看到，最后参数的估计量可以表示为随机样本的一个函数，这样的函数称为统计量(statistic)。 .. glossary:: 统计量 正式的，任意观测样本的实值函数 :math:`T=f(\mathcal{D})` 都称为一个 **统计量(statistic)** 。 一个统计量就是一个关于样本集的函数（允许是向量形式的函数)，**在这个函数中不能有任何未知参数** 。 比如，样本的均值 :math:`\bar{x}=\frac{1}{N}\sum_i^N x_i` ，最大值 :math:`max(\mathcal{D})` ， 中位数 :math:`median(\mathcal{D})` 以及 :math:`f(\mathcal{D})=4` 都是统计量。 但是 :math:`x_1+\mu` （ :math:`\mu` 是未知参数）就不是统计量。 在进行参数估计时，我们能利用的只有观测样本集，因此观测样本是我们进行参数估计的唯一信息源。 也就是说，我们能利用的有关参数的所有可用信息都包含在观察样本中。 **因此，我们获得的参数估计量始终是观测值的函数，即参数估计量是统计量。** 从某种意义上讲，该过程可以被认为是“压缩”原始观察数据：最初我们有N个数字， 但是经过这个“压缩”之后，我们只有1个数字 。 这种“压缩”总是使我们失去有关该参数的信息，决不能使我们获得更多的信息。 最好的情况是，该“压缩”结果包含的信息量与N个观测值中包含的信息量相同， 也就是该“压缩”结果包含的信息量已经是关于参数的信息的全部。 .. glossary:: 充分统计量 假设有一个统计量 :math:`T(\mathcal{D})` ，并且 :math:`t` 是 :math:`T` 的一个特定值， 如果在给定 :math:`T=t` 的条件下，我们就能计算出样本的联合概率 :math:`p(x_1,x_2,\dots,x_N|T=t)` ， 而不再依赖参数 :math:`\theta` ，这个统计量就是 **充分统计量(sufficient statistic)** 。 换种说法，在给定充分统计量 :math:`T=t` 条件下，就能确定参数 :math:`\theta` 的值，而不再需要额外的信息， 我们可以设想只保留 :math:`T` 并丢弃所有 :math:`x_i`，而不会丢失任何信息！ 从上面的直观分析中，我们可以看到充分统计量“吸收”了样本中包含的有关 :math:`\theta` 的所有可用信息。 这个概念是R.A. Fisher在1922年提出的。 充分性的概念是为了回答以下问题而提出的： 是否存在一个统计量，即函数 :math:`T(x_1,\dots,x_n)` ，其中包含样本中有关 :math:`\theta` 的所有信息？ 如果这样，则可以将原始数据减少或压缩到该统计信息而不会丢失信息。 例如，考虑一系列成功概率未知的独立伯努利试验。 我们可能有一种直觉的感觉，成功次数包含样本中有关 :math:`\theta` 的所有信息， 而成功发生的顺序没有提供有关 :math:`\theta` 的任何其他信息。 对于高斯分布，(样本)期望和(样本)协方差矩阵就是它的充分统计量，因为如果这两个参数已知， 就可以唯一确定一个高斯分布，而对于高斯分布的其他统计量，例如振幅、高阶矩等在这种时候都是多余的。 .. topic:: 示例： 令 :math:`x_1,x_2,\dots,x_N` 是N次独立伯努利实验的结果，其中 :math:`p(x_i=1)=\theta` 。 我们将验证 :math:`T=\sum_{i=1}^N x_i` 是 :math:`\theta` 的一个充分统计量。 .. topic:: 证明： 由于 :math:`x_i` 只能取值为0或者1，所以 :math:`T=t` 可以看作是在N条样本中 :math:`x_i=1` 的次数。 根据贝叶斯定理有： .. math:: p(x_1,x_2,\dots,x_N|T=t) &= \frac{p(x_1,\dots,x_N)}{p(T=t)} &= \frac{\prod_i \theta^{x_i} (1-\theta)^{(1-x_i)} }{p(T=t)} &= \frac{ \theta^t (1-\theta)^{N-t} }{p(T=t)} 现在看分母部分，T的含义是在N次实验中1的数量，很明显这是二项式分布，有N次试验，单词成功(为1)的概率为 :math:`\theta` ， 一共成功t次(1的数量为t)的概率分布为 :math:`T=\binom{N}{t}\theta^t (1-\theta)^{N-t}` ，其中 :math:`\binom{N}{t}` 是组合数，从N个结果中任意选出t个的方法数。把分母代入上式： .. math:: p(x_1,x_2,\dots,x_N|T=t) = \frac{ \theta^t (1-\theta)^{N-t} }{\binom{N}{t}\theta^t (1-\theta)^{N-t}} = \frac{1}{\binom{N}{t}} 最终发现，样本在给定 :math:`T=t` 的条件下的联合概率与参数 :math:`\theta` 无关，也就是说在确定了 :math:`T` 之后， 就可以直接得到样本的联合概率，而不再依赖参数 :math:`\theta` 。 在很多问题中，参数的最大似然估计量就是一个充分统计量，比如，伯努利实验的参数估计量就是一个充分统计量 :math:`\hat{\theta}_{ML}=\frac{1}{N}\sum_{i=1}^N x_i=\bar{x}` 。 同样，贝叶斯参数估计量也是一个充分统计量。最大似然估计量和贝叶斯估计量都是充分统计量的一个函数， 它们"吸收"了观测样本中关于参数的所有有用信息。 .. note:: 根据统计量的定义：样本的一个函数可以称为统计量，样本的求和 :math:`\sum_{n=1}^N x_n` ， 样本的均值 :math:`\frac{1}{N}\sum_{n=1}^N x_n` 都可以称为统计量。 所以，似然估计量 :math:`\hat{\theta}_{ML}=\frac{1}{N}\sum_{i=1}^N x_i` 可以整体看做一个充分统计量(均值统计)， 也可以看做是充分统计量(求和) :math:`\sum_{n=1}^N x_n` 的一个函数。 .. _ch_2_Fisher_Information: Fisher Information ############################################################################## 在参数估计问题中，我们从目标概率分布的观测样本中获取有关参数的信息。 这里有一个很自然的问题是：数据样本可以提供多少关于未知参数信息？ 本节我们介绍这种信息量的度量方法。 我们还可以看到， 该信息量度可用于查找估计量方差的界限， 并可用于近似估计从大样本中获得的估计量的抽样分布， 并且如果样本较大，则进一步用于获得近似置信区间。 假设有一个随机变量 :math:`X` ，其概率分布函数为 :math:`p(x;\theta)` ， :math:`\theta` 是模型的未知参数，并且 :math:`\theta \in \Theta` ， :math:`\Theta` 是参数空间。 对于一个随机变量 :math:`X \sim p(x;\theta)` ，当 :math:`\theta` 为真实值时，其似然函数应该取得一个极大值，或者说，此时似然函数的导数为0， 这是最大似然估计的基本原理。 我们定义 :math:`\ell(\theta;\mathcal{D})=log p(\mathcal{D};\theta)` 为对数似然函数(log-likelihood function)， :math:`\mathcal{D}` 可以看做是一个包含N的随机值的随机变量。 :math:`p(\mathcal{D};\theta)` 是所有样本的联合概率分布。 在很多资料中，符号 :math:`X` 既表示独立的随即变量，也可以是观测样本集变量， :math:`p(x;\theta)` 既是单个样本的概率分布，又是多个观测样本的联合概率（似然）， 容易造成混淆，不利于理解。 为了便于大家理解，我们分两种情况讨论： 1. 当仅有一个观测样本时，样本变量的概率分布函数为 :math:`p(x;\theta)` ， 对数似然函数为 :math:`\ell(\theta;x)=\log p(x;\theta)` 。 2. 当有 :math:`N` 个观测样本时，样本集变量的概率分布函数为 :math:`p(\mathcal{D};\theta)` ， 对数似然函数为 :math:`\ell(\theta;\mathcal{D})=\log p(\mathcal{D};\theta)= \sum_{i=1}^N \log p(x_i;\theta)` 。 我们先讨论只有一个观测样本的情形，此时对数似然函数的一阶导数为： .. math:: \ell(\theta;x) = \frac{\partial }{\partial \theta} \log p(x;\theta) =\frac{p'(x;\theta)}{p(x;\theta)} .. note:: 这里利用了对数函数的求导公式: .. math:: \nabla \log f(x) = \frac{1}{f(x)} \nabla f(x) 其中 :math:`p'(x;\theta)` 表示函数 :math:`p(x;\theta)` 关于 :math:`\theta` 的一阶导数，同理， :math:`p''(x;\theta)` 表示二阶导数。 通常模型的 :math:`\theta` 是一个向量，这时对数似然函数的一阶导数是一个向量， 二阶导数是一个方阵。 .. topic:: Score function 对数似然函数关于参数的一阶导数称为得分函数(Score function)，通常用符号 :math:`S` 表示。 .. math:: S(\theta)= \frac{\partial \ell(\theta)}{\partial \theta} 其中 :math:`\theta` 是模型的参数，当模型存在多个参数时，:math:`\theta` 是参数向量， :math:`S(\theta)` 也是一个向量。 :math:`S(\theta)` 是一个关于 **观测样本和参数** 的函数。 通常如果似然函数是凹(concave)的，我们可以通过令 :math:`S(\theta)=0` 求得参数的解析解。 :math:`S(\theta)` 是对数似然函数的一阶导数，一阶导数描述的是函数在这一点的切线的斜率， 导数越大切线斜率越大，所以 :math:`S(\theta)` 表示的是对数似然函数在某个 :math:`\theta` 值时模型的敏感度(sensitive)。 我们知道观测样本也是一个随机变量，:math:`S(\theta)` 是观测样本的一个函数， 所以 :math:`S(\theta)` 也可以看做是一个随机变量， 因此我们可以研究 :math:`S(\theta)` 的期望与方差。 我们先来看一下 :math:`S(\theta)` 的期望， :math:`S(\theta)` 的期望的计算需要利用一些数学技巧。 一个函数的积分和求导是可以互换的，并且概率分布函数的积分一定是等于1的， 所以有如下等式成立。 .. math:: :label: eq_2_33 \int f'(x;\theta) dx= \frac{\partial}{\partial \theta} \int f(x;\theta) dx =\frac{\partial}{\partial \theta} 1 =0 类似地有： .. math:: :label: eq_2_34 \int f''(x;\theta) dx= \frac{\partial^2}{\partial \theta^2} \int f(x;\theta) dx =\frac{\partial}{\partial \theta} 1 =0 :math:`S(\theta)` 关于样本变量的期望一定是等于0的， 结合 :eq:`eq_2_33` 可以推导出 :math:`S(\theta)` 的期望为： .. math:: \mathop{\mathbb{E}}_{p(x ; \theta)} \left[ s(\theta) \right] &= \mathop{\mathbb{E}}_{p(x ; \theta)} \left[ \nabla \ell(\theta;x) \right] &= \int [\nabla \ell(\theta;x)] \, p(x ; \theta) \, \text{d}x &= \int [\nabla \log p(x ; \theta)] \, p(x ; \theta) \, \text{d}x &= \int \frac{\nabla p(x ; \theta)}{p(x ; \theta)} p(x ; \theta) \, \text{d}x &= \int \nabla p(x ; \theta) \, \text{d}x &= \nabla \int p(x; \theta) \, \text{d}x &= \nabla 1 &= 0 :math:`S(\theta)` 的二阶矩(second moment)，也就是其方差(Variance)，被称为 Fisher Information， 中文常翻译成费歇尔信息， 通常用符号 :math:`I(\theta)` 表示， :math:`I(\theta)` 是一个方阵，通常称为信息矩阵(information matrix)。 .. math:: :label: eq_2_35 I(\theta) &=\mathop{Var(S(\theta))}_{p(x ; \theta)} &= \mathop{\mathbb{E}}_{p(x ; \theta)} [(S(\theta)- \mathop{\mathbb{E}}_{p(x ; \theta)} [S(\theta)] )^2] &= \mathop{\mathbb{E}}_{p(x ; \theta)} [S(\theta)^2] &= \mathop{\mathbb{E}}_{p(x ; \theta)} [S(\theta)S(\theta)^T] 实际上， :math:`I(\theta)` 和对数似然函数的二阶导数的期望值是有关系的， 我们先来看下对数似然函数的二阶导数， 二阶导数可以在一阶导数的基础上再次求导得到。 .. math:: \ell''(\theta;x) &= \frac{\partial}{\partial \theta} \ell'(\theta;x) &= \frac{\partial}{\partial \theta} \left [ \frac{p'(x;\theta)}{p(x;\theta)} \right ] &= \frac{p''(x;\theta)p(x;\theta)-[p'(x;\theta)]^2}{[p(x;\theta)]^2} &= \frac{p''(x;\theta)p(x;\theta)}{[p(x;\theta)]^2} - \left[ \frac{p'(x;\theta)}{p(x;\theta)} \right]^2 &= \frac{p''(x;\theta)}{p(x;\theta)} - [\ell'(\theta;x)]^2 然后我们看下对数似然函数二阶导数的期望值： .. math:: \mathop{\mathbb{E}}_{p(x ; \theta)} \left[ \ell''(\theta;x) \right] &= \int \left [ \frac{p''(x;\theta)}{p(x;\theta)} - [\ell'(\theta;x)]^2 \right ] p(x;\theta) dx &= \int p''(x;\theta) dx -\int [\ell'(\theta;x)]^2 p(x;\theta) dx &= 0 - \int [S(\theta)]^2 p(x;\theta) dx &= - \mathop{\mathbb{E}}_{p(x ; \theta)} [ [S(\theta)]^2 ] &= - I(\theta) 因此，Fisher Information 就等于对数似然函数二阶导数的期望的负数。 .. math:: I(\theta) = - \mathop{\mathbb{E}}_{p(x ; \theta)} \left[ \ell''(\theta;x) \right] 一个标量值函数(scalar-valued function)的二阶偏导数矩阵(方阵)称为海森矩阵(Hessian matrix)， 通常用符号 :math:`H` 表示， 因此 :math:`I(\theta)` 经常也被表示成海森矩阵的期望的负数。 .. math:: :label: eq_34_41 I(\theta) = - \mathop{\mathbb{E}}_{p(x ; \theta)}[H(\theta)] 我们看到，无论是通过score function 的方差计算，还是通过Hessian矩阵计算， :math:`I(\theta)` 都是一个期望值，所以经常被称为期望化信息(expected information)。 现在我们来看下观测样本为 :math:`N` 个的情形， 我们知道多个独立同分布观测样本的对数似然值就是每个独立样本的求和： .. math:: \ell(\theta;\mathcal{D})=\sum_{i=1}^{N} \ell(\theta;x_i) 我们用符号 :math:`I_{X}(\theta)` 表示一个观测值的信息矩阵，通常称为 unit Fisher information。 用符号 :math:`I_{\mathcal{D}}(\theta)` 表示多个观测值的信息矩阵，则有： .. math:: I_{\mathcal{D}}(\theta) = N I_{X}(\theta) :math:`I_{\mathcal{D}}(\theta)` 是正比于 :math:`N` 的，也就是说样本越多，我们的到关于参数的信息量就越大。 估计量的评价 ############################################ 一个服从某个概率分布的随机变量 :math:`X` ，假设其分布函数为 :math:`f(x;\theta)` ，其中 :math:`\theta` 是未知参数，我们可以从变量 :math:`X` 的观测样本集中估计出参数的值。估计的算法有很多种，我们已经介绍最大似然估计、贝叶斯估计、最大后验估计三种参数估计算法。 那么一个参数估计值 :math:`\hat{\theta}` 到底准不准呢？本节我们讨论如何评价一个估计量的好坏。 估计量的方差与偏差 ========================================================================== 对于一个长度为N的观测样本集 :math:`\mathcal{D}=\{x_1,x_2,\dots,x_N\}` ， 其中每条样本 :math:`x_n` 有 :math:`|\mathcal{X}|` 种可能取值，样本集 :math:`\mathcal{D}` 一共有 :math:`|\mathcal{X}|^N` 种可能取值。 某个特定样本集的概率为其中所有样本的联合概率 :math:`p(\mathcal{D})=p(x_1,x_2,\dots,x_N)` ， 显然，我们可以把 :math:`\mathcal{D}` 也可看做一个随机变量。 **我们知道参数估计量是通过观测样本得到的，所以参数估计量一定是观测样本集的一个统计量(函数)，** 长度固定为N样本集的不同采样将得到不同的参数估计值 :math:`\hat{\theta}` ，所以参数估计量也是一个随机变量。 我们用符号 :math:`\hat{\theta}(\mathcal{D})` 表示基于样本 :math:`\mathcal{D}` 的参数估计量。 随着样本数量的增加，一个参数估计量应该越来越接近参数的真实值，并且最终能收敛到参数的真实值， 我们把最终能收敛到真实值的估计量称为一致估计。 .. topic:: 一致性(Consistency) 当样本数量趋近于无穷大时，估计量 :math:`\hat{\theta}` 依据某种 **概率** 收敛于参数的真实值 :math:`\theta_{\text{true}}` ，那么这个估计量 :math:`\hat{\theta}` 就是一致性估计量。 .. math:: p(\hat{\theta}=\theta_{\text{true}})=1, n \rightarrow \infty 为什么是 *依概率* 收敛，而不是 *确定性* 收敛？ 因为参数估计量本身是一个随机变量，服从某种概率分布，只能是以某种概率得到某个确定性的值， 所以这里是依概率收敛到真实值。 一致性是对参数估计的基本要求，一个参数估计要是不满足一致性基本无用。 一致性估计量是 *依概率* 收敛到真实值的，并不是一定收敛到真实值， 所示我们实际上得到的参数估计量和真实值之间还是会存在一定误差的。 我们需要对这个误差进行量化评估，以便能评估一个估计量的好坏。 最直接的方法是计算估计量 :math:`\hat{\theta}` 和参数真实值 :math:`\theta_{\text{true}}` 之间的均方误差(mean-squared error,MSE)， 由于其中估计量是一个随机变量，所以我们计算MSE的期望。 .. math:: :label: eq_2_39 MSE &= \mathbb{E}_{\mathcal{D}} [(\hat{\theta}(\mathcal{D})- \theta_{\text{true}} )^2] &= \mathbb{E}_{\mathcal{D}}[ \hat{\theta}(\mathcal{D})^2-2\hat{\theta}(\mathcal{D})\theta_{\text{true}} + \theta_{\text{true}}^2 ] &= \left [ \mathbb{E}_{\mathcal{D}}[\hat{\theta}(\mathcal{D})^2]- \mathbb{E}_{\mathcal{D}}[\hat{\theta}(\mathcal{D})]^2 \right ] +\left [ \mathbb{E}_{\mathcal{D}}[\hat{\theta}(\mathcal{D})]^2 -2\hat{\theta}(\mathcal{D})\theta_{\text{true}} + \theta_{\text{true}}^2 \right ] &= \underbrace{Var_{\mathcal{D}} (\hat{\theta}(\mathcal{D}))}_{\text{方差}} + \underbrace{\left ( \mathbb{E}_{\mathcal{D}}[\hat{\theta}(\mathcal{D})]- \theta_{\text{true}}\right)^2}_{\text{偏差}} .. attention:: 这里的方差是参数估计量的方差，不是观测变量 :math:`X` 的方差。 根据上述MSE的等式，一个估计量和真实值之间的MSE是由方差和偏差决定的， 一个好的估计量应该是方差和偏差都尽可能的小。 一个估计量的偏差被定义成估计量的期望和参数真实值误差的平方： .. math:: b(\hat{\theta}) = (\mathbb{E}_{\mathcal{D}}[\hat{\theta}(\mathcal{D}) ] - \theta_{\text{true}})^2 .. glossary:: 无偏估计 当一个估计量满足 :math:`b(\hat{\theta})=0` 时，也就是满足 **估计量的期望值等于参数的真实值** ，就称这个估计量为无偏估计。 无偏估计就是偏差最小的估计量(偏差为0)，而估计量的方差也是存在下界的，这可以通过一个定理给出： .. topic:: Cramer-Rao Lower Bound (CRLB) 定理 Cram´er–Rao Lower Bound (CRLB)定理描述了一个确定性参数(deterministic parameter) :math:`\theta` 的估计量的方差的下界 .. math:: :label: eq_2_50 Var(\hat{\theta}) \ge \frac{(\frac{\partial}{\partial \theta} \mathbb{E}[\hat{\theta}] )^2}{I(\theta)} 其中分子是估计量期望对参数真实值的一阶导的平方， 如果一个估计量是无偏估计，那么有 :math:`\mathbb{E}[\hat{\theta}]=\theta_{\text{true}}` ，这时分子就等于1。 .. math:: (\frac{\partial}{\partial \theta} \mathbb{E}[\hat{\theta}] )^2 = ( \frac{\partial}{\partial \theta} \theta)^2 = 1 因此对于无偏估计量，:eq:`eq_2_50` 可以简化为： .. math:: Var(\hat{\theta}) \ge \frac{1}{I(\theta)} :math:`I(\theta)` 是费歇尔信息(Fisher-Information)矩阵。 根据CRLB定理，可以看出一个估计量的方差是存在下界的， 并且对于无偏估计量，估计量的方差的最小值是费歇尔信息的倒数。 显然当一个估计量的方差为下界时，这个估计量是最稳定的， 这时我们称之为有效估计。 .. topic:: 有效估计量(Efficient Estimator) 任意一个估计量，如果其方差为CRLB的下限，那么这个估计量是有效估计量。 最好的估计量应该是偏差和方差都尽可能的小，偏差最小为无偏估计， 所以我们定义出最小方差无偏估计。 .. topic:: 最小方差无偏估计 当参数 :math:`\theta` 存在多个无偏估计时，其中方差最小的估计量就称为最小方差无偏估计 (Minimum Variance Unbiased Estimator,MVUE)。 显然，MVUE是使得MSE最小的估计量。 然而，最小方差无偏估计量并不总是存在的， 即使存在，我们也可能找不到，没有任何一种方法会始终产生MVUE。 查找MVUE的一种有用方法是为参数找到充分统计量。 大数定律和中心极限定理 ===================================== 要理解后续的内容，依赖两个很重要的定理，大数定律(Law of Large Numbers,LLN)和中心极限定理(Central Limit Theorem,CLT)。 所以这里我们首先回顾一下这两个定理的内容。 假设有一个随机变量 :math:`X` ，假设其期望为 :math:`\mu` ，方差为 :math:`\sigma^2` 。 我们对这个随机变量进行N次独立采样，得到长度为N的观测样本集 :math:`\mathcal{D}=\{x_1,x_2,\dots,x_N\}` 。 样本集中的每个样本点都是从同一个随机变量(同样的概率分布)中采样得到，并且是独立进行的， 所以这个样本集是独立同分布的(i.i.d)的。 **样本集** :math:`\mathcal{D}` **可以看成是一个随机变量的N次采样，也可以把其中每个样本点看做是一个独立的随机变量，** **即** :math:`\mathcal{D}` **是N的同分布的随机变量集合。** **此时，我们表示成** :math:`\mathcal{D}=\{X_1,X_2,\dots,X_N\}` 。大数定律和中心极限定理都是描述的 :math:`\mathcal{D}` 的性质。 大数定律 ------------------------------------------------------- 我们定义样本集 :math:`\mathcal{D}` 的平均值为： .. math:: Z_N = \frac{X_1+X_2+\cdots+X_N}{N} 由于 :math:`x_1,x_2,\dots,x_N` 都是随机值，所以它们计算得到的 :math:`Z_N` 也是一个随机值， 所以 :math:`Z_N` 也是一个随机变量。 既然 :math:`Z_N` 也是一个随机变量， 那么我们就研究一起它的期望和方差。 :math:`Z_N` 的期望值可以通过下式计算得到： .. math:: \mathbb{E}[Z_N] &= \mathbb{E} \left [ \frac{X_1+X_2+\cdots+X_N}{N} \right ] &= \frac{\mathbb{E}[X_1+X_2+\cdots+X_N ] }{N} &= \frac{\mathbb{E}[X_1] + \mathbb{E}[X_2]+\cdots+\mathbb{E}[X_N] }{N} 由于 :math:`x_1,x_2,\dots,x_N` 是独立同分布的，假设随机变量 :math:`X_i` 的期望值为 :math:`\mu` ，于是可以得到 :math:`Z_N` 的期望值为： .. math:: \mathbb{E}[Z_N] = \frac{N \mu}{N} = \mu 我们发现 :math:`Z_N` 的期望值与每个独立变量 :math:`X_i` 的期望值相同。 现在我们看下 :math:`Z_N` 的方差。 .. math:: Var[Z_N] &= Var\left[ \frac{X_1+X_2+\cdots+X_N}{N} \right ] &= \frac{Var[X_1+X_2+\cdots+X_N ] }{N^2} 此时由于 :math:`x_1,x_2,\dots,x_N` 是独立的，所以满足：和的方差等于方差的和，即 .. math:: Var[X_1+X_2+\cdots+X_N ] =Var[X_1] + Var[X_2] + \cdots + Var[X_N] 因此 :math:`Z_N` 的方差为： .. math:: :label: eq_2_40 Var[Z_N] = \frac{Var[X_1] + Var[X_2] + \cdots + Var[X_N]}{N^2} = \frac{N \sigma^2}{N^2} = \frac{ \sigma^2}{N} 至此我们发现，对于独立同分布的样本，样本均值变量的期望值就等于独立变量的期望值，样本均值的方差为独立变量的方差的 :math:`1/N` 。 :math:`N` 越大，样本均值变量的方差 :math:`Var[Z_N]` 就越小，方差 :math:`Var[Z_N]` 越小， 样本均值变量 :math:`Z_N` 就越稳定，也就是越来越收敛于其自身的期望值 :math:`\mathbb{E}[Z_N]` 。 极限情况下，当 :math:`N` 无穷大时，方差趋近于0 :math:`Var[Z_N] \rightarrow 0` ， :math:`Z_N` 就稳定在(收敛于)期望值 :math:`Z_N \rightarrow \mathbb{E}[Z_N] =\mu` 。 现在我们给出大数定律的正式定义： .. topic:: 大数定律(Law of Large Numbers,LLN) 一个独立同分布的样本 :math:`X_1,X_2,\dots,X_N` ，如果 :math:`X` 的期望值是有限值，即 :math:`|\mathbb{E}[X]| < \infty` ，那么样本的均值 :math:`Z_N` *依概率(in probability)* 收敛于变量 :math:`X` 期望值。 .. math:: Z_N = \frac{X_1+X_2+\cdots+X_N}{N} \rightarrow^d \mathbb{E}[X] 为什么是 *依概率(in probability)* ？ 时刻牢记 :math:`Z_N` （样本均值）是一个随机变量，既然是随机变量就一定服从于某种概率分布。理论上，只有当N无穷大时， :math:`Z_N` 才收敛于 :math:`\mathbb{E}[X]` ；反之，当N有限时，仅满足 :math:`\mathbb{E}[Z_N] = \mathbb{E}[X]` 。 同时 :math:`Z_N` （样本均值）也是样本的一个统计量(statistic)。 中心极限定理 ------------------------------------------------------- 大数定律描述的是样本的均值 *"依概率"* 渐近收敛于分布的期望(均值)，并没有说明是依据什么概率分布， 而中心极限定理对此做出了解释。 回顾 :eq:`eq_2_40` ，我们已经知道样本均值变量的方差为： .. math:: Var[Z_N] = \frac{ \sigma^2}{N} 这个方差不是一个固定值，而是随着 :math:`N` 的增加趋近于0的，没有讨论的意义。 我们重新定义一个新的变量 :math:`W_N` ： .. math:: W_N = \sqrt{N}(Z_N-\mathbb{E}[X]) 依据大数定律，显然随着 :math:`N` 的增加 :math:`W_N` 是趋近于0的。 并且 :math:`W_N` 的方差为： .. math:: Var[W_N] = (\sqrt{N})^2 Var[Z_N] = \sigma^2 .. topic:: 中心极限定理(Central Limit Theorem,CLT) 一个独立同分布的样本 :math:`X_1,X_2,\dots,X_N` ，如果 :math:`X` 的期望值与方差都是有限值，即 :math:`|\mathbb{E}[X]| < \infty,\sigma^2=Var[X]<\infty` ，那么 :math:`\sqrt{N}(Z_N-\mathbb{E}[X])` 渐近服从于均值为0，方差为 :math:`\sigma^2` 的正态分布。 .. math:: \sqrt{N}(Z_N-\mathbb{E}[X]) \rightarrow \mathcal{N}(0,\sigma^2) .. _ch_2_MLE_estimator: 最大似然估计的特性 ===================================== 现在我们讨论下极大似然估计量的特性， 这里我们直接给出结论，省略证明过程。 极大似然估计量通常是满足中心极限定理以及大数定律的， 通常称为渐近正态性和渐近一致性。 .. topic:: 渐近正态性(Asymptotic normality) 我们说一个估计量是渐近正态性的，如果满足： .. math:: \sqrt{n}(\hat{\theta}) \rightarrow^d \mathcal{N}(\theta_{\text{true}},\frac{1}{I(\theta)}) 或者 .. math:: \sqrt{n}(\hat{\theta}-\theta_{\text{true}}) \rightarrow^d \mathcal{N}(0,\frac{1}{I(\theta)}) **渐近正态性对应着中心极限定理，极大似然估计是满足渐近正态性的。** 似然估计量不仅是渐近服从正态分布，而且是以参数真实值为均值的正态分布， 这表明似然估计量依概率(正态分布)收敛于参数的真实值，这符合一致性的定义。 显然极大似然估计是一致性估计。 由于似然估计是一致性估计，似然估计量是渐近收敛于参数真实值的，也就是估计量的偏差渐近为0， 因此可以得出似然估计量是 **渐近无偏估计** 。 我们知道估计量的MSE是由偏差和方差组成的（ :eq:`eq_2_39` ），无偏性是对估计量的偏差的评价， 而估计量的方差影响着估计值的稳定性，方差越小估计量就越稳定， MLE估计量的方差就等于费歇尔信息的倒数， 如果 :math:`\theta` 是参数向量，估计量的方差为协方差矩阵， 此时费歇尔信息 :math:`I(\theta)` 为信息矩阵(Information matrix)。 .. math:: :label: eq_2_50 \text{Cov}(\hat{\theta}_{ML}) = [I(\theta)]^{-1} 这里我们省略证明过程，有兴趣的读者可以参考其他资料。 显然似然估计不仅仅是渐近无偏估计，而且估计量的方差就等于CLRB定理的下界， 因此似然估计量是不仅仅是有效估计量，而是其最小方差无偏估计(Minimum Variance Unbiased Estimator,MVUE)， **并且我们可以通过** :math:`I(\theta)` **量化衡量MLE估计量的方差。** 当我们用MLE估计出一个参数的估计值后，我们期望能量化评估出这个参数估计值的好坏， 大家常用的方法是计算观测值的误差： .. math:: \text{MSE} = \frac{1}{N} \sum_i^N (y_i-\hat{y}_i)^2 这种方法衡量的是整个模型的预测效果，并不能衡量出参数的估计值和参数的最优值之间的误差， 我们已经知道MLE是无偏估计，那么最终MLE估计量的误差就可以用如下公式衡量： .. math:: \text{Standard Errors} = \sqrt{Var(\hat{\theta}_{ML})} = \sqrt{ \text{diag} ([I(\theta)]^{-1})} :math:`[I(\theta)]^{-1}` 是协方差矩阵，其对角线元素是每个参数方差，开根号后得到每个参数的标准差。 最后我们总结下MLE拥有的特性： - MLE估计量符合中心极限定理，满足渐近正态性(Asymptotic normality)。 - MLE估计量符合大数定律，是一致性(Consistency)估计。 - 由于满足一致性，渐近收敛于参数真实值，渐近偏差为0，因此MLE估计量满足渐近无偏性(Invariance)。 - MLE估计量的方差是信息矩阵的逆，符合CRLB定理的下限，所以是有效估计(Efficient Estimator)，并且是最小方差无偏估计。 **MLE的这些特性都是建立在样本数量无穷大的条件下，当样本数量很小时，是没有这些特性的。** **此外，贝叶斯估计由于增加了先验信息，不再是无偏估计，而是有偏估计。** **在样本数量比较小时，极大似然估计与贝叶斯估计互有优劣，但随着样本数量的增加， 极大似然估计和贝叶斯估计是相似的。**