EM算法最通俗理解

期望最大化算法（EM）是机器学习领域非常重要的算法之一，但作为一个工科生，每次看其推导过程，总会怀疑自己的智商是不是不够用，为什么每一步推导都能看懂，但放到一起就崩了呢？可能还是因为数学基础太差了，总不能真正理解，尤其是不理解为什么要这样，相信工科的同学都有这种疑问：到底是先有了实际意义再将其公式化，还是先有了公式再赋予其意义呢？前者应该是工科生的思维，而后者可能只有数学大牛才能搞定吧。

今天我们就以工科生的思维来理解一下EM算法，抛开数学公式推导，看能不能从实际意义上出发，最终得到类似的结果呢？从此刻开始，请忘记所有公式推导，忘记EM算法，回归问题本身。
让我们从经典的男女身高问题开始：这里有一笔学生身高统计数据（假设100个），但是统计的时候忘了统计性别，现在要分别估计男生和女生的身高分布。
你可能会说，这数据是男是女都不知道，根本没办法统计啊！但是现在告诉你，如果你做不出结果，老板就会把你开除，女朋友就要跟你分手。啥意思呢，就是无论如何你要做出一个结果！(言外之意，结果是否100%正确并不重要，统计本来就是一个概率的问题，只要你的结果"差不离儿"就可以了)
好了，为了不被开除/被分手，只能硬着头皮上了。这里先统一一下基本的常识，也就是根据经验来总结几条：
1) 男生身高一般来说是高于女生的
2) 从统计学来看，身高应该是服从正态分布的(也就是我们要估计的其实是4个数，女均值/女方差、男均值/男方差)
基于上面两条经验，我们可以得到一个最暴力的版本：
版本一：快刀斩乱麻
把所有身高从小到达排序一下，前50个作为女生身高，后50个作为男生身高，然后根据概率论里学到求均值和方差的公式，就可以求出男生和女生的分布了。
问题就这么简单地解决了，但这个做法显然过于暴力了，能不能做得细致一些呢？这种直接一刀切的方式，最大的问题是切的点可能不对，因为没有人说男生女生的数量是一样的。那怎么切能更准确一点儿呢？我们先来看一眼统计数据，假如一个身高是155，那你自己判断下它是属于男生还是女生呢？这还用说吗，这个八成是女生的啊。那你是基于什么判断的呢？显然这里有个经验问题，即女生的身高基本是160左右，男生基本是180左右（不接受任何反驳）。所以靠近160的统计数据，大概率是属于女生的，而靠近180的数据，大概率是属于男生的。有了这个经验我们就可以做出第二个版本了：
版本二：根据经验的一刀切
样本数据小于170((160+180)/2)的按照女生处理，大于170的按照男生处理。然后同样根据均值方差公式，得到男生女生的分布。
这个版本连先验知识（也就是没遇到这个问题之前你就知道的知识），感觉已经够极限了。但是还不够，老板跟女朋友还是逼问你：能不能做得更好呢？你可能有点儿发疯了，我这人生经验都用上了，这还能怎么搞啊？那我们换个问法：如果告诉你的结果不够准确，你觉得最有可能是哪里出了问题呢？显然，你就做了一个规定那就是按照170划分，那肯定就是划分错了呗，这里确实有点儿不太"优雅"，因为169跟171其实就差了2厘米，却划分给了不同的组。甚至如果正好有个170的，该划分给男生还是划分给女生呢？你说那怎么办呢？难道还要把它劈开吗？哎，对了，为啥不能劈开呢？人当然是不能劈开的，数据是可以劈开的啊，也就是男女占比就按照(0.5/0.5)，这样是不是更好呢？然后再看一下169的，前面也说了，它应该是属于女生的概率稍微大一些，那我们是不是可以认为它的男女占比是(0.4/0.6)。既然有了这样的思路，那190的其实也可以认为是男女占比是(0.9/0.1)嘛。
把一个样本按照概率分给男女两个分布去统计，这是个关键点。但前面也说了，这也是没办法的办法，可以说完全是基于直觉的。
于是我们就有了一个新的问题，那就是如何计算这个比例（为了避免混淆，后文以"分配比例"说明）。其实这个分配比例代表的是一个样本属于男还是属于女的概率的比。那我们就求出这两个概率不就可以了吗，用什么求概率呢，用概率密度公式啊，概率密度公式从哪儿来，从变量的分布来啊，那分布从哪儿来......哎，我们要求的不就是分布吗，怎么转了一圈又回到了原点？
回到了原点是不是说明无解了呢？稍等一下，谁说男女分布不知道呢？我们在版本一版本二中，不都得到了一个男女身高的分布吗？假如我们把这个分布拿过来不就算出了这个分配比例吗，那用分配比例又可以计算一个新的分布了.....事情有点儿变得有意思了，我们来捋一下：
我们拿一个预先估计的分布，计算出来一个样本的分配比例，然后用这个比例计算出一个新的分布。由于第一次预估的分布是强行划分得到的，第二次则用了一个分配比例，那第二次得到的分布直觉上应该更准确一些，更关键的是，我们得到了一个自循环，因为我们又可以用新的(应该是更准确的)分布去计算一个更准确的分配比例，这样迭代下去，问题就变成了估算分布->计算分配比例->计算分布->计算分配比例...
版本三：按比例划分样本
1)我们预先估计一个初始概率分布，可以用一些最简单的统计方式（比如方式一），甚至直接根据经验给出（比如前面说的男生均值就是180，女生就是160这样）
2)用估计的概率分布把样本按照比例划分，再根据新的比例去计算一个新分布
3)不断重复执行过程2
那问题来了，这个过程会一直进行下去吗？应该是不会，因为如果我们的估算越来越准确，那我们迟早会跟真实值非常接近，这个时候两次的计算结果应该就非常接近了。也就是这个循环肯定会有一个极限值存在的，不过极值一般是对一个函数而言的。哎对了，这个地方好像能隐约能感觉到我们是在求一个函数的极值有不有？虽然这个函数我们不知道是啥。如果你愿意，那你可以试着定义一个这样的函数，如果你定义了这样一个函数，那我愿意替你的函数起个名字：就叫似然函数吧。既然模糊感觉到的函数都有名字了，那我们用了半天的这个分配比例是不是也该有个名分呢？分配比例在概率论中的专业名词叫分布（注意不要跟我们估计的分布混了）。既然分布都有了，那总得定义一个变量吧，有变量才能有分布啊，那这个用来表示样本分配比例的变量叫啥呢？就叫它隐变量吧。

以上就是对男女身高问题的分析，既然有了分析，肯定要实验一下，用代码产生两组服从正态分布的数据，然后混合在一起，看用我们的方法能不能得到正确的结果。

N = 1000 #样本数量
# 模拟产生男女身高N个
garray = np.random.normal(loc = 165, scale = 5.0, size = N)
barray = np.random.normal(loc = 175, scale = 5.0, size = N)

# 混合数据
merged_array = []
for i in range(len(male_array)):
    merged_array.append(male_array[i])
    merged_array.append(fmale_array[i])
    
# 根据经验给出初始参数
gmu, gsigma = 160, 10 #女均值、标准差
bmu, bsigma = 180, 10 #男均值、标准差

# 计算该参数下的权重分配
gweight = []
bweight = []
for h in height:
    gw = stats.norm.pdf(h, gmu, gsigma)
    bw = stats.norm.pdf(h, bmu, bsigma)
    gweight.append(gw / (gw + bw))
    bweight.append(bw / (gw + bw))
    
# 用新权重更新参数(女生的)
wsum = 0.0
gmuSum = 0.0
for h,w in zip(height, gweight):
    gmuSum += h * w
    wsum += w

gmu = gmuSum / wsum

这里要注意一下，有个比例归一的问题，假设没有比例分配问题，那女生的平均值就是gmuSum/N，但这里的gmuSum都乘了一个比例，假如所有比例都是(0.1/0.9)，gmuSum实际是乘了个0.1的。所以，后面应该除上一个0.1。那如果比例不相等呢，那就除一个比例的平均值嘛。将所有比例求和再除N,即 wsum/N 然后用 gmuSum/N 除这个数，两个N消除了，就得到gmu = gmuSum / wsum

# 再算下标准差（女生的）
gsigmaSum = 0.0
for h,w in zip(height, gweight):
    gsigmaSum += (h - gmu) * (h - gmu) * w
gsigma = gsigmaSum / wsum
gsigma = math.sqrt(gsigma)
 

同样的方式，计算男生的

wsum = 0.0
bmuSum = 0.0
for h,w in zip(height, bweight):
    bmuSum += h * w
    wsum += w

bmu = bmuSum / wsum

bsigmaSum = 0.0
for h,w in zip(height, bweight):
    bsigmaSum += (h - bmu) * (h - bmu) * w
bsigma = bsigmaSum / wsum
bsigma = math.sqrt(bsigma)

这样，一次运算就完成了，可以用新的参数去估算新的比例了，程序里加个循环就OK啦。后面会给出完整的程序代码，大家可以直接跑的，但在此之前，我们还有三个问题要考虑，这三个问题如果考虑明白了，那基本就算是理解了
1) 如果我们在一定范围内改变给定的初始值，结果会怎么样呢？
2) 假如我不小心把男女均值写反了，那结果会是怎么样的？
3) 假如我给定的男女均值是一样的(比如都是170)，结果会是怎样呢？
你可以先考虑下这三个问题的结果，然后用下面的程序验证一下，或者也可以直接在程序上实验，先看结果再分析
....
PS：如果你看过别人或自己实现过EM算法，你会发现我的代码跟EM算法是完全一样的流程，说白了，这其实就是EM算法(为了回归主题，我们还是叫回它原来的名字吧)
假如你已经做完了实验，我们来尝试回答一下刚才的三个问题：
1) 如果我们在一定范围内改变给定的初始值，结果会怎么样呢？
如果你做了实验就会发现，在一定范围内改变初值，最后的结果是一样的。当然如果你改的太多了，比如改成了身高均值是10000，那可能会得到一些报错，这是因为程序的计算精度是有限的，如果中间的结果太小了，程序里就按照0处理了，所以可能除0异常。但在理论上，如果我们的计算精度是无限的，那无论初值是什么，我们最终都可以得到正确结果。
2) 假如我不小心把男女均值写反了，那结果会是怎么样的？
这是个有趣的问题，如果你把男女身高均值搞反了，那最终的结果也是反的，这里面好像有点儿说法，我们先看最后一个问题，再一起分析。
3) 假如我给定的男女均值是一样的(比如都是170)，结果会是怎样呢？
不知道这个结果有没有让你很惊讶，如果我们给的初值一样，虽然它相当合理，相当接近真实值，但算法就是得不到正确结果。无论迭代多少次，男女的分布参数都是完全一样的。但只要我们稍微改变一点点，比如是170/170.1，那随着迭代进行，估计的均值和标准差都会慢慢的接近真实值，只是速度会很慢，可以预想到，如果迭代次数足够多，我们最终还是会得到正确结果。
如果你真的思考了，那你肯定已经有一些想法了，看看跟我们想的是不是一致呢？
EM算法其实是一种聚类的算法，它之所有能跑起来，并能得到一个不错的结果，主要依据是样本的“扎堆”性，因为样本分布其实是分了两堆的，我们可以把它想象成一排铁钉，然后我们选择的初始化的值就像两个磁铁一样，会被铁钉所吸引。然后磁铁就会向铁钉分布的中心移动，随着磁铁的移动，每个铁钉对它的引力也会产生变化，最终经过一段时间的移动，磁铁最终会到达铁钉堆的中心。而之所以两个磁铁会分开，就是因为它们的初始位置不一样，我们又强迫每个铁钉对磁铁的吸引力总和是一样的，这样就产生了蝴蝶效应，靠近下面的磁铁收到下面铁钉的引力就大于上面的铁钉，这样慢慢就将两个磁铁分开了。而磁铁本身是没有区分的，也就是两个磁铁换个位置，那结果就正好反过来了。
最后再提一个问题吧，如果我们统计的不是男生女生的身高分布，而是按照年龄段统计。比如年龄段在 5-10、10-15、15-20 的分别统计其分布，这时候该怎么处理呢？
给出一点提示：建立模型的时候，不是样本分了几堆，而是你认为它应该分几堆。借用投资大师芒格的一句话，“如果你有一把锤子，你看什么都像钉子”。