On this page
参数估计问题的回顾
比如有一枚不均匀的硬币，它抛出之后正面朝上和反面朝上的概率并不相等，你想去探索硬币正面朝上的概率，那很简单，你把它抛 $N$ 次，其中正面朝上的次数为$n_1$，那么我们就用 $\frac{n_1}{N}$ 的值作为正面朝上概率的估计。
又比如某个地区 18 岁成年男子的身高服从一个正态分布，那么如果想知道这个身高正态分布的均值和方差，那也不难，我们可以随机挑选 $N$个该地区的成年男子作为样本，分别统计出他们的身高 $x_1,x_2,x_3,…,x_N$，通过极大似然估计法得到的样本均值和样本方差：
\[\mu_{mle}=\frac{1}{N}\sum_{i=1}^{N}x_i\] \[\sigma^2_{mle}=\frac{1}{N}\sum_{i=1}^n(x_i-\mu_{mle})^2\]如果想要估计值具备无偏性，按照之前讲过的知识点，做一点简单的修正即可：
\[\sigma^2=\frac{1}{N-1}\sum_{i=1}^n(x_i-\mu_{mle})^2\]新情况：场景中含有隐变量
如果在实际场景中，情形发生一点变化，你就会发现似乎前面非常有效的极大似然估计方法就会变得使不上劲了。
比如，现在手上有两枚硬币，硬币 A 和硬币 B，它们都是不均匀的，但是它们正面朝上的概率 $p_A$ 和 $p_B$并不相等，每次你从这两枚硬币随机摸出一枚硬币进行投掷（但是你并不知道取的是哪一枚），你依然记录每次抛掷的结果是正面还是反面，正面出现的次数 $N_1$，反面出现的次数 $N_2$，但是你能够直接依据这个数据估计出硬币 A 和硬币 B 各自正面朝上的概率吗？恐怕是不行的吧。
再比如，还是统计某地区 18 岁青年的身高，但是这次的样本中，既有男生，也有女生，那么当你获取了 $N$ 个青年的样本身高之后，能够直接去求得男生的身高均值和方差以及女生的身高均值和方差吗？显然是不行的，因为男生和女生的身高各自服从不同参数下的高斯分布。
为什么这两个情况下，直接去用极大似然估计的方法去估计参数就不行了呢？是因为这两个场景下，模型都变成了混合模型，不光有观测变量，更出现了一个隐含变量。
一般而言，第 $i$ 个样本的观测变量记作 $x_i$，隐含变量记作 $z_i$。
比如抛硬币的例子中，每次抛硬币的正反是观测变量，但是到底这次抛的是硬币 A 还是硬币 B，我们不知道，这是一个隐含变量。
又比如量身高的例子，我们每次看到的身高值是观测变量，而这个样本到底是属于男生的还是女生的，就是一个隐含变量。
迭代法：解决含有隐变量情形的抛硬币问题
那么这种混合模型，不能直接用极大似然估计法去估计，那我们就考虑采用迭代法去慢慢地试，这就是我们贯穿后续好几讲内容的 EM 方法。
说实话整个方法的公式和计算技巧有点复杂，这一讲当中，我们暂时不引入公式，我们针对抛硬币的例子，来做一个简单的计算，目的是让大家感性的认识一下，在有隐变量的情况下，迭代探索参数究竟是一个什么样的过程。我们在这一讲当中不涉及具体的公式。
好，我们下面具体来描述一下抛硬币的背景。
我们有不均匀的硬币 A 和硬币 B，我们重复做 5 组试验，每一次试验的过程是这样的，我们任取其中一枚硬币，连续抛掷 10 次，记录正反的次数，5 组试验的结果如下：
| 组别 | 正面次数 | 反面次数 |
|---|---|---|
| 第一组 | 5 | 5 |
| 第二组 | 9 | 1 |
| 第三组 | 8 | 2 |
| 第四组 | 4 | 6 |
| 第五组 | 7 | 3 |
那么，硬币 A 和硬币 B 正面向上的概率 $\theta_A$ 和 $\theta_B$ 各是多少？
目前问题的症结就在于我们不知道每一组到底扔的是哪一枚硬币，我们先插一句嘴，如果我现在告诉你，第一组和第四组试验我们用的是硬币 B，第二、三、五组试验用的是硬币 A，那你会做吗？
那就好办了，对于硬币 A 的参数估计，我们直接拿第二、三、五组试验数据进行即可：
\[\theta_A=\frac{9+8+7}{30}=0.8\]然后那第一、四组试验去估计硬币 B：
\[\theta_B=\frac{4+5}{20}=0.45\]而这里，一步肯定是估计不出参数的，我们使用迭代法，也就是准备用多轮操作逐步接近真实的 $\theta_A$ 和 $\theta_B$。
我们先给一个初值 $\theta_A^{(0)}$ 和 $\theta_B^{(0)}$，我们先随意给一个值：
\[\theta_A^{(0)}=0.6，\theta_B^{(0)}=0.5\]那么基于这个预设的参数，在第一组试验中我们设硬币 A 出现的概率为 $P_{A1}$，硬币 B 出现的概率为 $P_{B1}$。
如果在试验 1 中，我们完全用硬币 A 来抛，抛出 5 正 5 反的概率为：
\[\theta_{A}^{(0)5}(1-\theta_{A}^{(0)})^5=0.6^5\cdot0.4^5=0.0007962624\]同理在试验 1 中，我们完全用硬币 B 来抛，抛出 5 正 5 反的概率为：
\[\theta_{B}^{(0)5}(1-\theta_{B}^{(0)})^5=0.5^5\cdot0.5^5=0.0009765625\]那么第一组试验中，是硬币 A 和硬币 B 的概率分别为：
\[P_{A1}=\frac{0.6^5\cdot0.4^5}{0.6^5\cdot0.4^5+0.5^5\cdot0.5^5}=0.45\] \[P_{B1}=\frac{0.5^5\cdot0.5^5}{0.6^5\cdot0.4^5+0.5^5\cdot0.5^5}=0.55\]按照同样的方式，我们可以分别计算出另外四组试验中，硬币 A 和硬币 B 各自出现的概率：
| 组别 | 硬币 A | 硬币 B |
|---|---|---|
| 第一组 | 0.45 | 0.55 |
| 第二组 | 0.80 | 0.20 |
| 第三组 | 0.73 | 0.27 |
| 第四组 | 0.35 | 0.65 |
| 第五组 | 0.65 | 0.35 |
我们还是用第一组试验中算出来的 $P_{A1}=0.45$ 和 $P_{B1}=0.55$ 来做做文章。
第一组试验的结果是 5 正 5 负，那么由硬币 A 扔出来的结果应该是：0.45*5=2.25 正，0.45*5=2.25 负；由硬币 B 扔出来的结果应该是：0.55*5=2.75 正，0.55*5=2.75 负，由此可以得到各组试验中硬币 A 和硬币 B 各自抛掷出来的结果：
| 组别 | 硬币 A | 硬币 B |
|---|---|---|
| 1 | 正：2.25，负：2.25 | 正：2.75，负：2.75 |
| 2 | 正：7.2，负：0.8 | 正：1.8，负：0.2 |
| 3 | 正：5.84，负：1.46 | 正：2.16，负：0.54 |
| 4 | 正：1.4，负：2.1 | 正：2.6，负：3.9 |
| 5 | 正：4.55，负：1.95 | 正：2.45，负：1.05 |
| 共 | 正：21.24，负：8.56 | 正：11.76，负：8.44 |
有了这个表，我们就能做够重新估计硬币 A 和硬币 B 正面朝上的概率了，因为此时数据是完备的，估计出的值作为新一轮的估计值：
\[\theta_A^{(1)}=\frac{21.24}{21.24+8.56}=0.71\] \[\theta_B^{(1)}=\frac{11.76}{11.76+8.44}=0.58\]就好像当初的 $\theta_A^{(0)}$ 和 $\theta_B^{(0)}$ 那样，此时我们有了 $\theta_A^{(1)}$ 和$\theta_B^{(1)}$，就可以把上面两张表的计算过程不停地重复，不断地去计算和更新出
\[\theta_A^{(2)},\theta_B^{(2)};\theta_A^{(3)},\theta_B^{(3)};\theta_A^{(4)},\theta_B^{(4)};...\theta_A^{(N)},\theta_B^{(N)}\]直到参数收敛于真实的 $\theta_A,\theta_B$。
代码试验
我们用代码来演示一下上述过程，观察一下最终参数的收敛情况。
代码片段：
import numpy as np
from scipy.stats import binom
#一轮迭代处理
def single_iter(theta_priors, observations):
"""
param observations：五组试验的观测值
param theta_priors：上一轮迭代形成的 theta_A 和 theta_B
"""
counts = {'A': {'H': 0, 'T': 0}, 'B': {'H': 0, 'T': 0}}
theta_A = theta_priors[0]
theta_B = theta_priors[1]
#迭代计算每组试验的数据
for observation in observations:
len_observation = len(observation)
num_heads = observation.sum()
num_tails = len_observation-num_heads
P_A = binom.pmf(num_heads, len_observation, theta_A)
P_B = binom.pmf(num_heads, len_observation, theta_B)
# 计算出硬币 A 和硬币 B 各自出现的概率
weight_A = P_A / (P_A + P_B)
weight_B = P_B / (P_A + P_B)
# 更新在当前硬币 A 和硬币 B 各自出现的概率下，硬币 A 和硬币 B 各自的正反面次数
counts['A']['H'] += weight_A * num_heads
counts['A']['T'] += weight_A * num_tails
counts['B']['H'] += weight_B * num_heads
counts['B']['T'] += weight_B * num_tails
#重新估计新一轮的硬币 A 和硬币 B 正面向上的概率
new_theta_A = counts['A']['H'] / (counts['A']['H'] + counts['A']['T'])
new_theta_B = counts['B']['H'] / (counts['B']['H'] + counts['B']['T'])
return [new_theta_A,new_theta_B]
observations = np.array([[1,0,0,0,1,1,0,1,0,1],
[1,1,1,1,0,1,1,1,0,1],
[1,0,1,1,1,1,1,0,1,1],
[1,0,1,0,0,0,1,1,0,0],
[0,1,1,1,0,1,1,1,0,1]])
prior = [0.6, 0.5] #设定初始的参数值
iteration = 0
iterations = 10000 #最多的迭代次数
tol = 1e-6 #判断参数收敛的阈值
while iteration < iterations:
new_prior = single_iter(prior,observations)
print(new_prior)
delta_change = np.abs(prior[0] - new_prior[0])
if delta_change < tol:
break
else:
prior = new_prior
iteration += 1
print('迭代结束，总共迭代轮数{}'.format(iteration))
print('最终估计的参数{}'.format(new_prior))
运行结果：
[0.683267379440028, 0.5794860533160178]
[0.6946902774724677, 0.5830445317245763]
[0.6995256100328493, 0.57971899072605]
[0.7033758035611858, 0.5759155969304315]
[0.7067629289547314, 0.5723818307416716]
[0.7097151528294361, 0.5692497378243309]
[0.7122309725095934, 0.5665503627275396]
[0.7143282388567959, 0.5642771950286414]
[0.7160432311533907, 0.5624008414937195]
[0.7174231039725665, 0.5608782115958729]
[0.7185187176686263, 0.5596601027467699]
[0.7193794191585418, 0.5586969396036339]
[0.7200499095520126, 0.5579425231081682]
[0.7205688014147885, 0.5573560461279302]
[0.720968331704508, 0.5569028257103324]
[0.7212747545382558, 0.5565542085397206]
[0.7215090638835501, 0.5562870189502978]
[0.7216878211235749, 0.5560828082205115]
[0.7218239594521363, 0.5559270662698881]
[0.7219275026279905, 0.5558084843857315]
[0.7220061755631375, 0.5557183096081342]
[0.7220659062448442, 0.5556498025455018]
[0.7221112291752515, 0.5555977947465911]
[0.7221456045736407, 0.5555583344445941]
[0.7221716680847203, 0.555528407038558]
[0.7221914245178972, 0.5555057168235281]
[0.722206397253179, 0.555488517848247]
[0.7222177429417977, 0.5554754835858867]
[0.7222263392688235, 0.5554656069366276]
[0.7222328519363013, 0.5554581237464508]
[0.7222377856908556, 0.5554524544514625]
[0.722241523141501, 0.5554481596321429]
[0.7222443542577712, 0.5554449062077694]
[0.7222464987677085, 0.5554424417501024]
[0.7222481231543918, 0.5554405749813793]
[0.7222493535477021, 0.5554391609763221]
[0.72225028549926, 0.5554380899384829]
迭代结束，总共迭代轮数 36
最终估计的参数 [0.72225028549926, 0.5554380899384829]
我们从程序运行的结果来看，经过 36 轮的迭代，最后收敛出 $\theta_A=0.722$，$\theta_B=0.555$。
从这个例子中，我们初步领略了如何用迭代求解的方法，来探索包含隐含变量情况下的参数估计问题，说实话这个例子非常简单也非常直观，所以我们很清爽地处理了，实际上没有看到太多的数学表达式。
但是实际上，这种极其简单的情况并不多，大部分含有隐变量的问题计算都较为复杂，我们需要归纳出一套公式、一个形式化的统一办法，来描述参数$\theta^{(t)}$ 和 $\theta^{(t+1)}$的迭代关系，那么在这几讲当中，我们就详细介绍一套非常常用的混合模型参数迭代的理论和方法——EM方法，深入剖析这套方法如何使用、为何管用，最终基于它再次回到实践，用 EM 方法解决混合高斯模型中的参数估计问题。