2. 最大似然估计¶
在统计学中,把需要调查或者研究的某一现象或者事物的全部数据称为统计总体,或简称 总体(population)。 比如,我们要研究中国人的身高分布,那么全国14亿人的身高数据就是总体(population), 这14亿身高数据所属的数据分布称为 总体分布 (population distribution), 其中每一个人的身高数据,即单个数据称为个体(individual)。 然而在实际中,我们不可能得到14亿的全部数据,也就是 总体数据通常是无法得知的 。 这时,可以选择抽样(sampling),即从总体当中随机抽取出部分个体,然后得到这部分抽样个体的数据, 一次抽样的结果称为一份样本(sample)。 比如,从14亿的人群中随机抽取出1万的个体,然后去测量这1万人的身高数据, 这样就得到了一份包含1万个数据的样本,样本的容量(sample size),或者说样本的大小,是1万。 注意样本(sample)和个体(individual)的区别,样本(sample)是一次抽样的结果,包含多个个体(individual)数据, 一份样本中包含的个体数据的数量称为本容量(sample size)。 通常我们会假设总体分布服从某种已知的概率分布,但是分布的某些参数是不确定的, 比如全国身高数据服从正态分布,但是期望和方差不知道, 这时我们期望能通过样本推断(估计)出总体正态分布的期望和方差参数。
推断统计学(或称统计推断,英语:statistical inference), 指统计学中,研究如何根据样本(sample)数据去推断总体(population)特征(或者参数)的方法, 比如根据样本的平均值去估计总体的均值参数。 它是在对样本数据进行描述的基础上,对统计总体的未知数量特征做出以概率形式表述的推断。 更概括地说,是在一段有限的时间内,通过对一个随机过程的观察来进行推断的。 在统计学中,利用样本推断(估计)总体分布参数方法有很多, 比如矩估计、最大似然估计、贝叶斯估计等等, 本章我们讨论其中应用最为广泛的最大似然估计算法。
2.1. 最大似然估计¶
最大似然估计(Maximum Likelihood Estimation,MLE),又叫极大似然估计,是统计学中应用最广泛的一种未知参数估计方法。 它可以在已知随机变量属于哪种概率分布的前提下, 利用随机变量的一些观测值估计出分布的一些参数值。 所谓观测值,就是随机变量的采样值,也就是这个随机变量试验的真实结果值, 因为是我们能”看到”的值,所以称为观测值。
假设有一个离散随机变量 \(X\) ,其概率质量函数是 \(P(X;\theta)\) ,其中 \(\theta\) 是这个概率分布的参数,其值是未知的。 函数 \(P(X;\theta)\) 本身是已知的,也就是我们知道 \(X\) 所属何种概率分布,比如是高斯分布等等。
现在假设我们有一些变量 \(X\) 的观测值,这些观测值集合用符号 \(\mathcal{D}=\{x^{(1)},x^{(2)},\ldots,x^{(N)}\}\) 表示。这些观测值都是从同一个概率分布 \(P(X;\theta)\) 得到的, 并且这些样本是独立获取的,即每条样本值不依赖其它样本值, 我们可以称这些样本是 独立同分布 的。
我们知道其中任意一条样本 \(x_i\) 的发生概率是 \(P(x_i;\theta)\) , 那么所有样本发生的联合概率是 \(P(\mathcal{D};\theta)=P(x^{(1)},\ldots,x^{(N)};\theta)\) , 又由于所有样本是满足独立同分布的(i.i.d)的,根据联合概率分布的分解法则有
假设 \(\theta\) 的可能取值空间为 \(\Theta\) ,记作 \(\theta \in \Theta\) 。 不论 \(\theta\) 取何值,都有一定的可能(概率)产生出这个样本集 \(\mathcal{D}\) ,但显然 \(\theta\) 的值会影响着这个样本的产生概率 \(P(\mathcal{D};\theta)\) 。换句话说就是,不同的 \(\theta\) 值会得到不同的样本联合概率 \(P(\mathcal{D};\theta)\) 。
现在我们思考 \(\theta\) 真实值是什么。 事实上,我们根本无从得知参数 \(\theta\) 的真实值。 但我们可以换个思路,我们可以从 \(\theta\) 的取值空间 \(\Theta\) 中 挑一个最好的 的出来。 那么什么是最好的,这个最好的标准是什么?
在概率统计中,把观测样本的联合概率称为 似然(likelihood), 一般用符号 \(L(\theta;\mathcal{D})=P(\mathcal{D};\theta)\) 表示, 有时也称为似然函数(likelihood function)。
观测样本集的似然函数就是样本集的联合概率
最优的 \(\theta\) 值是令观测样本发生概率最大的值, 也就是令似然函数取得最大。 参数 \(\theta\) 的最大似然估计值可以写为
仔细观察后发现,似然函数是每条样本概率 \(P(x_i;\theta)\) 的连乘, 而概率值都是在 \([0,1]\) 之间的,一系列小于 \(1\) 的数字连乘会趋近于 \(0\)。 而计算机在处理浮点数时存在精度问题,太小的值是无法表示的。 所以一般我们会为似然函数加上一个对数操作来解决计算机的精度问题, 我们把加了对数的似然函数称为 对数似然函数(log-likelihood function) , 一般用符号 \(\ell\) 表示。
通过极大化对数似然函数 \(\ell(\theta;\mathcal{D})\) 得到 \(\hat{\theta}\) 和极大化似然函数 \(L(\theta;\mathcal{D})\) 是等价的,这里不再证明,有兴趣的读者可以参考其他资料。
虽然这里我们是以离散随机变量为例,但最大似然估计同样可以应用于连续值随机变量的参数估计。 连续值随机变量用的是概率密度函数函数表示其每个状态的概率大小情况, 概率密度函数表示是每一个点的”密度”,而不是概率值, 但每个点的密度是和它的概率呈正比的。 假设连续值随机变量 \(X\) 的概率密度函数是 \(f(x;\theta)\) ,则有
最大似然估计是通过极大化对数似然函数求解,对于连续值随机变量用概率密度函数 \(f(X=x;\theta)\) 替换 \(P(X=x;\theta)\) ,对极大化求解没有任何影响。 因此在使用最大似然估计概率模型的分布时,如果是离散随机变量就用概率质量函数, 如果是连续值随机变量就是概率密度函数。
那么如何进行极大化求解呢?通常有如下三种方法:
解析法(Analytic),又叫直接求解法。我们知道一个函数在取得极值时其一阶导数是为 \(0\) 的, 因此可以通过令对数似然函数的一阶导数为 \(0\) 得到一个方程等式,然后解这个方程得到 \(\hat{\theta}_{ML}\) 。 这种方法得到的解称为解析解。
(2.1.7)¶\[\frac{\partial \ell}{\partial \theta} = 0\]函数的一阶导数为 \(0\) 的点称为“驻点”(stationary point),可能为(局部)极大或者极小值点,也可能为鞍点(saddle point), 可以通过极值点的二阶导数判断是极大值点还是极小值点。 并不是所有情况都能得到解析解的,很多时候是无法直接求得的,在后面的章节中我们会详细讨论。
网格搜索法(Grid Search)。如果我们知道 \(\hat{\theta}\) 的值在空间 \(\Theta\) 中, 可以对这个空间进行搜索来得到使得似然函数最大的参数值。 换句话说,就是尝试这个空间中的每个值,找到令似然函数取得最大的参数值。 网格搜索方法是一种很好的方法,它表明可以通过重复逼近和迭代来找到似然函数的最大值。 但是,它在大多数情况下不切实际,并且当参数数量变多时变得更加困难。
数值法(Numerical)。这是现在最常用的算法。本质上就是先为 \(\theta\) 赋予一个初始值, 然后利用爬山法找到最优解。梯度下降(上升)法(Gradient descent),牛顿法(Newton-Raphson),BHHH,DFP等等都属于这类。 关于这类算法,读者可以先参考其他资料,我暂时没有精力写。
本章我们只使用解析法求解,在正式讲广义线性模型时再介绍最大似然估计的数值求解法, 下面介绍几个应用最大似然估计的具体例子。
2.2. 伯努利分布¶
假设一个随机变量 \(X\) 服从伯努利分布(Bernoulli distribution), 即只有两种可能的取值 \(X \in \{0,1\}\), 设其取值为 \(1\) 的概率为 \(P(X=1)=\theta\) , 其概率质量函数为
其中 \(\theta\) 是未知的参数,需要使用最大似然估计得到。 假设变量 \(X\) 的独立同分布的观测样本集为 \(\mathcal{D}=\{x^{(1)},\ldots,x^{(N)}\}\) , 样本集的规模为 \(|\mathcal{D}|=N\) 。
现在我们利用最大似然估计法估计出参数 \(\theta\) , 首先写出观测样本的对数似然函数。
为方便表述,我们定义几个统计值,在样本集 \(\mathcal{D}\) 中, \(n_0\) 表示随机变量 \(X=0\) 在观测样本中的次数, \(\hat{p}_{\mathcal{D}}(0)=\frac{n_0}{N}\) 表示 \(X=0\) 在观测集中出现的相对频率(经验分数); \(n_1\) 表示 \(X=1\) 在样本中的次数, \(\hat{p}_{\mathcal{D}}(1)=\frac{n_1}{N}\) 表示 \(X=1\) 在样本中出现的相对频率(经验分数)。 把 \(n_0,n_1\) 代入到对数似然函数中可得
我们知道当对数似然函数的导数为 \(0\) 时,函数取得极值。 现在对对数似然函数求导,并令导数为 \(0\) 。
又由于 \(n_0 = N - n_1\) ,代入上式可得:
化简可得 \(\theta\) 的估计值
假设投硬币 \(N\) 次,正面向上的次数是 \(n_1\) , 在你朴素的认知里,你认为正面向上的概率是多少?仔细观察这个结果,是不是和你的经验认知是一样的。
2.3. 类别分布¶
现在假设随机变量 \(X\) 是类别变量,其取值空间是 \(\mathcal{X}=\{x_1,\dots,x_K\}\) , 类别分布的概率质量函数可以写成下面的形式
其中 \(\theta=[\theta_1,\dots,\theta_K]\) 为分布的参数。 公式(2.1.21) 的含义是,变量 \(X\) 取值为类别 \(x_k\) 的概率是 \(\theta_k\) , 即 \(P(X=x_k;\theta)=\theta_k\) 。其中参数向量 \(\theta\) 需要满足约束:
现在利用最大似然估计出参数向量 \(\theta\) 的值, 我们继续用符号 \(\mathcal{D}\) 表示观测样本集, \(|\mathcal{D}|=N\) , 则似然函数为
其中 \(n_k\) 表示类别 \(x_k\) 在样本中出现的次数,那么有 \(\sum_{k=1}^K n_k=N\) 。 我们为似然函数加上对数操作,以便把连乘符号转换成加法。
为了找到 \(\theta_k\) 的最大似然解,我们需要最大化对数似然函数 \(\ell(\theta;\mathcal{D} )\) , 并且限制 \(\theta_k\) 的和必须等于1。这样 带有约束的优化问题需要借用拉格朗日乘数 \(\lambda\) 实现,即需要最大化如下方程。
对上述公式求偏导,并令偏导等于 \(0\)。
把 \(\theta_k = - \frac{n_k}{\lambda}\) 代入到约束条件 \(\sum_{k=1}^K \theta_k=1\) 中, 解得 \(\lambda = -N\) ,最终可求得 \(\theta_k\) 的值:
我们发现,伯努利分布与类别分布的参数最大似然估计值具有相同的形式,并且最大似然估计值是可以通过样本统计的得到, 这是最大似然估计的一个特性。 样本的统计值被称为统计量(statistic), 统计量(statistic)是 样本的一个函数,其代表着从样本中提取的一些”信息”,比如样本的均值(mean),样本的总和(sum)等等。 很多时候这些信息可以用于确定这个分布的未知参数, 如果仅需要一个统计量就能确定这个分布的未知参数,而不再需要其它的额外”信息”,那么这个统计量就称为这个分布(或者分布族) 的 充分统计量(sufficient statistic) ,在后面的章节中我们会详细讨论充分统计量。
2.4. 高斯分布¶
现在我们假设有一个高斯随机变量 \(X \sim N(\mu,\sigma^2)\) ,其参数有两个,均值 \(\mu\) 和方差 \(\sigma^2\) 。观测样本集为 \(\mathcal{D}=\{x^{(1)},x^{(2)},\ldots,x^{(N)}\}\) ,我们利用最大似然估计出参数 \(\mu,\sigma^2\) ,首先写出高斯分布的概率密度函数。
似然函数为
高斯分布的概率密度函数中有自然常数 \(e\) 的指数, 因此其对数似然函数,我们选择以常数 \(e\) 为底数的对数。
然后对参数求偏导数,并令偏导数为0。
由第一个等式可以解的:
其中 \(\bar{x}\) 表示样本的平均值,然后代入第二个等式解得:
至此,我们就得到了高斯分布的均值参数和方差参数的最大似然估计值。 仔细观察可以发现,参数的估计值只依赖两个观测样本的统计量。
这两个量被称为高斯分布的充分统计量(sufficient statistics), 如果你对统计量和充分统计量不理解,没关系,下一章我们会详细讨论。
2.5. 总结¶
最后我们来总结下用最大似然方法估计参数的一般模式。
令 \(X\) 表示随机变量, 令 \(f(X;\theta)\) 为随机变量 \(X\) 的概率质量(密度)函数的参数化表示, 其中 \(\theta\) 是未知参数向量 (注意 \(\theta\) 表示是多个参数的集合,不一定只有一个未知参数)。 我们把 \(\theta\) 看做是一个 非随机变量 ,是一系列数值变量,但是其具体取值却是未知的。
我们给出这个变量的独立同分布观测样本集,\(\mathcal{D}=\{x^{(1)}, \ldots, x^{(N)}\}\) ,上标代表样本编号。 目标是利用样本估计出参数 \(\theta\) 的值。 一个在给定观测时估计未知参数的方法是使用最大似然估计,也叫极大似然估计。
\(f(\cdot;\theta)\) 是一个参数为 \(\theta\) 的概率分布的概率质量函数(pdf)或者概率密度函数(pmf)。
\(L(\cdot ; \mathcal{D}) \triangleq P(\mathcal{D} ; \cdot)\) 是随机变量 \(X\) 的观测数据集 \(\mathcal{D}\) 的似然函数(likelihood function) 。
- 最大似然估计(maximum likelihood estimation,MLE)¶
最直接的理解是:当使得给定观测数据的似然函数取得最大值时,此时似然函数中未知参数的取值是最优的。
似然函数的概率解释是:这些观测值(观测值是可观测到的随机变量的取值)”同时(不是时间上的同时,是联合发生)”发生的概率。 最大似然就是:既然这些样本事件已经发生了,那么我就假设他们的发生的概率是最大的,就认为使得这些样本具有最大发生概率时参数的值是最优取值。 注意,最大似然估计的解不一定存在,也不一定唯一,这取决于似然函数是否有极值点,以及有几个极值点。
观测样本集的对数似然函数是
然后通过极大化对数似然函数的方式求解参数的值。当然有些时候可以通过令偏导数为 \(0\) 直接求得解析解, 然而有时候却无法得到解析解,需要用梯度迭代的方法求解。
虽然最大似然估计使用十分广泛,但是它不是完美的,在样本较少时,或者样本有偏时,得到的估计值偏差较大。 例如投掷一个普通的硬币3次,每次都是正面朝上,这时最大似然估计正面朝上的概率时结论会是1, 表示所有未来的投掷结果都是正面向上。这明显是有问题的,当数据集(观测样本集)较少时非常容易出现错误的结果, 当然也不是没有解决办法,比如可以使用贝叶斯估计,贝叶斯估计可以算是最大似然估计的升级版, 通过增加先验的方式解决这种极端的场景。有兴趣的读者可以参考其他资料了解贝叶斯估计。 通过这个小例子让我们意识到,使用最大似然估计的得到的参数估计值未必是符合我们心意的, 这个估计值到底行不行,我们需要一种手段来评价参数估计值的好坏, 下一章我们讨论如何评价最大似然估计值的好坏。
最后我们讨论一下几个细节的地方,一个是关于独立同分布的,一个是关于对数的。 首先讨论独立同分布的问题。前面在讨论用最大似然估计概率分布的未知参数时,我们假设观测样本是独立同分布的, 这里”独立”和”同分布”分开讲。”同分布”顾名思义,所有观测样本都是同一个概率分布的,这个显而易见的。 我们要估计某个概率分布的未知参数,必然需要这个概率分布的观测样本,其他概率分布的观测样本也没用啊。 再看”独立”性,这里指所有样本之间互相独立,没有依赖关系。 如果所有样本之间独立,根据第一章讲的链式法则,那么它们的联合概率就可以拆解成每个样本独立发生概率的连乘, 这样联合概率的计算就变得十分简单。
如果样本之间”不独立”,根据链式法则,联合概率就需要分解成一系列条件概率的乘积,显然这样变复杂很多。 实际上样本独立性并不是必须得,与你的概率模型相关,比如马尔科夫链,其中既有独立也有不独立,这里暂时不过多讨论了。
接下来讨论对数问题。首先是为什么加对数,或者说加对数的好处的是什么? 前文已经提到了一点,加对数后连乘变成连加,避免了计算机无法处理极小浮点数的问题。除此外,还有一个好处是求导变简单了。 似然函数求极值的过程需要求导,乘法求导是十分困难的,加法就简单了很多。 最后一个问题,对数的底数有什么特别要求没有?答案是没有。理论上,底数是多少都行,但如果你的概率分布是指数族, 即概率质量(密度)函数是自然常数 \(e\) 的指数形式,此时用以 \(e\) 为底的对数是十分合适的, 正好把 \(e\) 给抵消掉。有关指数族分布是本书的重点,后面的章节会详细讨论。