最大似然估计
公司年会的抽奖环节需要一个抽奖程序，部门A承担了这个抽奖程序的编写。抽奖结果出来后，发现100个奖项中的35个都被部门A的成员抽到了，而其它部门的中奖数量与与部门的人数占比是一致的；但是部门A的成员人数只占公司总人数的10%。此时大家肯定会认为抽奖程序是不公平的，并且偏向了部门A。
这里涉及了概率与统计的两种应用，假设检验和最大似然估计。
我们想当然地认为抽奖程序是公平的，那么每个人抽到奖项的概率是一样的，最可能出现的抽奖结果是每个部门抽到的奖项与部门的人数占比是一致的，至少不会差异太。这就是假设检验的零假设。零假设给出了各种情况出现的概率，如果出现概率极低的情况出现了就可以否定零假设，此处就是否定了程序是公平的假设。也就是说认为程序是不公平的就是假设检验的应用。
在知道A部门中奖概率为35%的情况下，对抽象程序的真实情况进行了推测——偏向部门A。因为直觉上我们认为只有程序偏向部门A才能说得通，才最可能导致这个结果。这就是极大似然估计的一种应用。
似然函数与概率函数
“最大似然估计”中的最大似然的意思是，在给定的分布模型下这个结果出现的概率最大，估计的意思就是求得此时分布模型的参数。可见似然也是概率，之所以叫做似然只是一种约定。通常我们说概率的时候，表示的是不同的结果在分布模型下的取值。此时结果已经出现了，如果仍然采用在结果出现之前给定的参数，这个结果的概率就是确定的。通过假设检验知道了之前给定的参数是不对的，需要估计新的参数，也就是将参数当作未知的。对于不同的参数，结果将会取得不同的概率值。可见似然函数与概率函数的模型是一样的，区别只在于将谁看作变量，将谁看作参数。
如果将概率函数记作\(p(x;\theta)\)，那么似然函数应当记作\(p(\theta;x)\)，此种函数记法中分号左侧的是变量右侧是参数。对于概率函数更一般的记法为\(p(x|\theta)\)，此时似然函数也记作\(p(x|\theta)\)，这实际上是条件概率的记法，将参数\(\theta\)视作已经出现的条件。
\[ p(x|\theta) = \begin{cases} p(x;\theta), & \text {概率函数} \\ p(\theta;x), & \text{似然函数} \end{cases} \]
为什么参数可以视作先验条件这件事先不管，先说说为什么概率函数和似然函数用条件概率记法时都是\(p(x|\theta)\)？
因为概率函数和似然函数都表示的是在参数条件\(\theta\)成立的条件下，结果\(x\)出现的概率。而\(p(\theta|x)\)表示将参数视作先验条件时，条件\(\theta\)的后验概率，意义为在\(x\)出现后对\(\theta\)概率的再评估。根据贝叶斯公式，有：
\[p(\theta|x) = \frac{p(\theta) \times p(x|\theta)}{p(x)}\]
显然\(p(\theta|x)\)和\(p(x|\theta)\)的函数模型是不一样的。（在贝叶斯公式中\(p(x)\)表示在所有需要考察的不同\(\theta\)出现时，\(x\)出现的概率总和。\(p(\theta)\)表示不同\(\theta\)出现的概率。如果需要考察的\(\theta\)就是全部可能的\(\theta\)，那么\(\int{p(x)}d_x=1\)，仍然有\(p(x)\le1\)。通常套用贝叶斯公式时，都有\(p(x)<1\)。如果\(p(x)=1\)，说明\(x\)只有一个取值，这个取值就是\(x\)的全体，\(\theta\)属于\(x\)的子事件时\(p(x|\theta)=1\)，否则\(p(x|\theta)=0\)。）
回到文章开头给出的例子，此时似然函数或者概率函数应该是什么样的呢？如果不考察部门，别人中奖了，你没中奖，并不能说明抽奖程序对你不公平；因此单考察每个人的中奖情况时没有任何可用的信息对程序的公平性做出判断。将人归属的部门纳入考虑时，我们做出了程序在部门之间的公平性判断，所基于的是部门A实际中奖人数与期望的中奖人数不一致。因而我们的似然函数或者概率函数应该是部门A的中奖人数在给定的部门中奖概率\(\theta\)之下的概率。
\[p(x|\theta) = \begin{pmatrix} 100 \\ x \\ \end{pmatrix}\theta^x(1-\theta)^{100-x}\]
如公式所示，这是一个典型的二项分布。对应到经典的箱子和球时，可以设定部门A为红球，其它部门为白球，已经知道红球和白球的比率为\(\theta\)，即每次抽奖抽到红球的概率为\(\theta\)。从箱子里一次取100个球出来，取到\(x\)个红球的概率就是上式。
当然实际抽奖时是一个一个的抽的，一般都相当于无放回的从箱子里取球。无放回的取球和有放回的取球对于考察某一个球这样的个体来说是有区别的，此时考察的是某一类球这个群体的个数而不是个体，因而是否有放回是无关紧要的。可以换一种思路，无放回的从箱子里一个一个的取100个球，与用一个大铲子一下子铲出来100个球，效果是一样的。但是如果考虑了群体中个体的组合关系，也就是还要考察某种球是第几次抽到的，就不是二项分布了。为了简单起见，就假装是一次把100个奖全抽了吧，这相当于只进行了一次二项分布实验。
似然函数的最大值
似然函数仍然是函数，约定了先验条件\(\theta\)的不同取值与\(x\)的某个确定结果的概率之间的映射关系。为了解出\(\theta\)的取值，必须要有一个概率。我们确实也有了一个概率，将35和0.1带入上式求出来的，再将这个概率带入求出来的\(\theta\)仍然等于0.1，这样做是没有任何意义的。就是因为认为将35和0.1带入上式求出来的概率不可信，太低了，才要重新求解参数\(\theta\)的取值。
既然认为根据先验参数解出来的概率太低了，那么就应该用更大的概率带入反算\(\theta\)。到底应该是多大呢？就让这个概率取到它能取到的最大值吧。对应到我们人类的直觉就是，如果某个事件发生了，那么这个事件应该是所有可能发生事件中概率最大的那个。但是这个最大值是未知的，但是通过微积分的知识可以知道，如果在某个区间内函数图像是凸的，那么最大值和极大值取在同一处，且如果可导的话该处的导数为零。因而我们可以对似然函数求导，令其等于零，然后解出\(\theta\)。这也是最大似然估计也被称作极大似然估计的原因。
\[\frac{{\rm d}p(x|\theta)}{{\rm d}\theta} = 0\]
要能够利用上面的公式，似然函数必须要满足几个条件。
- 似然函数在\(\theta\)的取值区间内函数图形是凸的（凹函数）；
- 似然函数在\(\theta\)的取值区间内处处可导；
- 令似然函数等于零有解析解。
如果似然函数的函数图像是凸的，也处处可导，但是就是没有解析解，可以设目标函数为\(-p(x|\theta)\)，然后用梯度下降算法找到似然概率取到最大值时的\(\theta\)。比如似然概率函数的图像是一个三角尖。
如果似然函数的函数图像是一条水平线，即\(x\)是均匀分布的，此时是不能使用该方法求解参数的，因为不管\(\theta\)取任何值，\(x\)的概率都是一样的。对于均匀分布的假设，只能通过假设检验否定这个假设，而不能用最大似然估计求解分布模型的参数。对于本文给出的例子，实际上隐含了一个均匀分布，那就是每个人中间的概率是一样的。一般遇到这种情况，可以将原始的分布进行分段，然后转换成伯努利分布，进而转换成二项分布。此处，我们就是这么做的。最终得到的结果是分段的不均匀情况。
多样本与乘法公式
在前面讲述中，我们是通过实验或者观察取得了\(x\)某一个结果\(x_i\)（样本），将其带入\(p(x|\theta)\)得到似然函数。仅通过一个样本对分布模型的参数进行估计未免过于武断。更多情况是，进行多次实验或者观察多次，此时将会得到多个样本。此时该如何操作呢，显然不能让每个样本都取到分布能出现的最大值。此时需要换一种看待问题的方式，将多个样本看作一个事件\(X\)，这个事件的意义是观察到样本都发生。那么\(X\)的概率函数可以套用乘法公式。
如果\(x\)是离散型的，\(p(x|\theta)\)就是\(x\)的发生概率，当抽得的样本为\(x_1,x_2,...,x_n\)时，样本总体\(X\)的发生概率为：
\[P(X|\theta)=P(x_1,x_2,...,x_n|\theta)=\prod_{i=1}^n p(x_i|\theta)\]
因为\(x_1,x_2,...,x_n\)都是已知的，将\(\theta\)当做未知的变量，所以这个也就是似然函数。
如果\(x\)是连续型的，\(p(x|\theta)\)不是\(x\)的发生概率，给样本值\(x_1,x_2,...,x_n\)加上一个邻域\(dx_1,dx_2,...,dx_n\)，其中\(dx_i>0\)，又足够小以至于所有的\(\left[x_i,x_i+dx_i \right]\)没有重叠。将抽取到样本\(x_i\)看作“样本的取值落在范围\(\left[x_i,x_i+dx_i \right]\)中”这样一个事件\(X_i\)发生，那么对应的概率为
\[P(X_i|\theta) = \int_{x_i}^{x_i+dx_i}p(x|\theta) dx \approx p(x_i|\theta) \times dx_i\]
进而\(X_i\)的总体发生概率为
\[P(X|\theta)=P(X_1,X_2,...,X_n|\theta)\approx\prod_{i=1}^n p(x_i|\theta) \times \prod_{i=1}^n dx_i\]
又由于\(dx_i\)是给定的微小正常数，不改变“令导数为零”的解，因而连续型的似然函数也被写做（实际上也不改变函数的单调性，使用梯度下降算法时似然函数也可以使用同一个）
\[P(X|\theta)=P(X_1,X_2,...,X_n|\theta)=\prod_{i=1}^n p(x_i|\theta)\]
知道了似然函数，接下来就可以直接套用公式解了。让我们再次回到本文开头给出的例子，继续转换对问题的看待方式。就按照实际抽奖情况来看待，相当于无放回的取球。假设公司总员工数量为100吧，部门A有10个人，此时第一次抽到部门A的概率为0.1，那么第\(n\)次的概率是多少？不妨先看看第2次的。第二次的概率由两种情况构成，第一次抽到部门A和第一次没抽到部门A。
第一次抽到部门A时，第二次还抽到部门A的概率为
\[p_1=\frac{\begin{pmatrix} 10 \\ 1 \\ \end{pmatrix}}{\begin{pmatrix} 100 \\ 1 \\ \end{pmatrix}} \times \frac{\begin{pmatrix} 9 \\ 1 \\ \end{pmatrix}}{\begin{pmatrix} 99 \\ 1 \\ \end{pmatrix}}\]
第一次抽到部门B时，第二次抽到部门A的概率为
\[p_2=\frac{\begin{pmatrix} 90 \\ 1 \\ \end{pmatrix}}{\begin{pmatrix} 100 \\ 1 \\ \end{pmatrix}} \times \frac{\begin{pmatrix} 10 \\ 1 \\ \end{pmatrix}}{\begin{pmatrix} 99 \\ 1 \\ \end{pmatrix}}\]
把他们加到一起，并提取公因式得到
\[p_2=\frac{\begin{pmatrix} 10 \\ 1 \\ \end{pmatrix} \times\left( \begin{pmatrix} 90 \\ 1 \\ \end{pmatrix} + \begin{pmatrix} 9 \\ 1 \\ \end{pmatrix} \right)}{\begin{pmatrix} 100 \\ 1 \\ \end{pmatrix} \times \begin{pmatrix} 99 \\ 1 \\ \end{pmatrix}}=\frac{\begin{pmatrix} 10 \\ 1 \\ \end{pmatrix}}{\begin{pmatrix} 100 \\ 1 \\ \end{pmatrix}}\]
第三次抽取到时，前两次的抽取有\({\begin{pmatrix} 2 \\ 1 \\ \end{pmatrix}}^2\)种情况；第\(n(n\le10)\)次抽取时有\({\begin{pmatrix} 2 \\ 1 \\ \end{pmatrix}}^{n-1}\)中情况，穷举显然是不合适的，不过观察式子的结构，可以相信第\(n\)次抽取时抽到部门A的概率仍然为0.1，不相信的可以自己动手多算算，算多了就信了。
设每次抽到部门A这个事件为1，其概率为\(\theta\)，没抽到部门A这个事件为0，其概率就是\(1-\theta\)，概率函数如下
\[ p(x|\theta) = \theta^x(1-\theta)^{1-x} = \begin{cases} \theta, & \text {if x= 1} \\ 1-\theta, & \text{if x=0} \end{cases}\]
可以看出每一次抽奖都是一个典型的伯努利分布。可能有人会对这个公式表示异议，因为某次抽奖的结果出来之后，剩下的人抽奖结果肯定都会增加的啊。确实是这样的，在已经知道抽奖结果的情况下，再判断剩余人的抽奖概率已经不属于概率（不准确）了，是后验概率。我们现在建立的是抽奖之前对中奖情况预估的概率模型，抽奖程序不是在抽奖过程中编写的不是，而是在抽奖之前编写的吧。最终检验的也是抽奖之前对中奖情况预估的概率模型，最终估计的也是抽奖之前对中奖情况预估的概率模型的先验参数。
将抽奖过程看作多次伯努利实验时，就不能仅仅考虑群体个数了，还需要考虑群体中个体的组合关系。按照时间先后的顺序抽奖了100次，那么个体的组合关系也只有一种，也就是100个伯努利实验结果发生了，那么他们的联合概率函数应该应用乘法公式。
\[P(X|\theta)=P(x_1,x_2,...,x_n|\theta)=\prod_{x_i=1} \theta^1(1-\theta)^{1-1}\times \prod_{x_i=0} \theta^0(1-\theta)^{1-0}\]
设抽到部门A的个数为\(n\)，那么100次伯努利实验结果的联合概率函数就是
\[P(X|\theta) = \theta^n(1-\theta)^{100-n}\]
此时用最大似然估计将会估计到和二项分布一样的参数，但是有100次实验结果，是不是能把我们对\(\theta\)估计结果的信息提升100备呢！
PS: 二项分布是不考虑群体中个体的组合关系的多重伯努利分布。