引入
设袋子里有红白两种球，比例为\(p:1-p\)。现在不知道比例值是多少，因此可以随机抽样：\(x_1,x_2...x_n\)，每一个\(x_i\)都是随机变量的随机抽样结果，并进行如下估计：\(p=\frac{|\{x_i=红\}|}{n}\)
如果有两个袋子，各自都有红色的球和白色的球，抽取的时候先随机挑选一个袋子，设选第一个袋子的概率为\(w\)，第一个袋子中\(红球:白球=p:1-p\)，第二个为\(红球:白球=q:1-q\)。
如果此时进行抽取，抽取到红球的概率为：\(P\{x=红\}=w\times p+(1-w)\times q\)，抽到白球的概率同理。
如果想要知道从第一个袋子中抽到红球的概率：\(P\{第一个|x_i=红\}\)，可以把概率乘以权重放到矩阵中：
$$\begin{pmatrix}
wp&w(1-p)\
(1-w)q&(1-w)(1-q)\
\end{pmatrix}$$
取矩阵的第一列分析，不难得出\(P\{第一个|x_i=红\}=\frac{wp}{wp+(1-w)q}\)。
此时模型的\(w,p,q\)三个参数均未知，可以使用极大似然法对其进行估计。令1代表红色，0代表白色，每个\(x_i\)的密度函数为\(wp^{x_i}(1-p)^{1-x_i}+(1-w)q^{x_i}(1-q)^{1-x_i}\)，如果\(x_i=1\)(红色)，带入上式可以发现结果就是\(wp+(1-w)q\)，也就是矩阵第一列的和。
根据极大似然估计，需要把所有样本的密度函数做乘积，似然函数就是：
\(\Pi_{i=1}^n [wp^{x_i}(1-p)^{1-x_i}+(1-w)q^{x_i}(1-q)^{1-x_i}]\)，需要对其求极大值。如果按照极大似然的思路，会发现这个问题很复杂（中括号里的部分不易处理，不信你试试）。不过这个函数可以转换为如下的形式：
\(\Pi_{i=1}^n [(wp+(1-w)q)^{x_i}(w(1-p)+(1-w)(1-q))^{1-x_i}]\)，这直观上也很容易理解，因为\(x_i\)只有0和1两种形式（带入\(x_i=1 \ or\ 0\)试试）。对其处理的话，可以令\(P=wp+(1-w)q\)，\(Q=w(1-p)+(1-w)(1-q)\)，函数继续变形为\(\Pi_{i=1}^nP^{x_i}Q^{1-x_i}=\Pi_{i=1}^nP^{x_i}(1-P)^{1-x_i}\)，对其取对数后得\((\Sigma_{i=1}^nx_i)lnP+(n-\Sigma_{i=1}^nx_i)ln(1-P)\)，对P求导并令导数为0，可得\(P=\frac{\Sigma_{i=1}^nx_i}{n}=\frac{k}{n}\)，\(k\)为样本中1（红球）的个数，且有\(w(1-p)+(1-w)(1-q)=\frac{n-k}{n}\)（本质上等于上面的方程）。
总结一下，刚才设计的模型（先随机选袋子，再随机选球）从极大似然估计的角度，\(wp+(1-w)q=\frac{k}{n}\)就是最佳估计，但是一个方程里含有三个未知数，所以只要三个参数满足这个方程就是一个合理的估计，即给出两个参数就可以把另一个参数解出来。
从另一个角度看问题
现在给出三个参数\(p^{(0)}, q^{(0)}, w^{(0)}\)，对于\(x_i\)，要根据它的颜色判断它来自于哪个袋子（利用上面的矩阵来计算）：
红球来自第一个袋子的概率：\(\frac{wp}{wp+(1-w)q}=a_i\)，白球来自第一个袋子的概率：\(\frac{w(1-p)}{w(1-p)+(1-w)(1-q)}=b_i\)...同理设红球和白球来自第二个袋子的概率为\(c_i,d_i\)，那么对于每个球，这四个概率都可以求出来，列成一张表：
| \(x_1\) | \(x_2\) | .. | ||
|---|---|---|---|---|
| 第一个-红 | \(\frac{wp}{wp+(1-w)q}=a_1\) | |||
| 第一个-白 | \(\frac{w(1-p)}{wp+(1-w)(1-q)}=b_1\) | |||
| 第二个-红 | \(\frac{(1-w)p}{wp+(1-w)q}=c_1\) | |||
| 第二个-白 | \(\frac{(1-w)(1-p)}{wp+(1-w)(1-q)}=d_1\) |
每一列的总和都是一样的（都是1），因为这里考虑的是后验概率，即每个球已知是红色还是白色，所以对于每一列只是把和已知颜色对应的行加起来（而非把所有的行加起来）。
因此来自于第一个袋子的概率\(w\)可以给出一个重新的估计：\(w^{(1)}=\frac{\Sigma_{x_i=1}a_i+\Sigma_{x_i=1}b_i}{n}\)，即表中前两行相加再除以n。这个过程就是：选定\(w^{(0)}\)为初值，根据初值计算出一个个概率值，再根据上表的概率值计算出新的\(w^{(1)}\)。
同理\(p^{(1)}=\frac{\Sigma_{x_i=1}a_i}{\Sigma_{x_i=1}a_i+\Sigma_{x_i=0}b_i}\)，\(q^{(1)}=\frac{\Sigma_{x_i=1}c_i}{\Sigma_{x_i=1}c_i+\Sigma_{x_i=0}d_i}\)。如此可以继续迭代下去得到新的参数值，能够证明这个算法是收敛的，至少可以收敛到局部极大值（不一定是全局）。
整个过程就叫做EM算法。
推广
对于高斯分布\(N(\mu,\sigma^2)\)，设有\(K\)个高斯分布，从中随机抽取到的概率分别为\(w_1,w_2...w_k\)，且满足\(\Sigma_{i=1}^k=1\)。
\(X\)的密度函数\(w_1f_1(x)+....+w_kf_k(x)\)，\(f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\)
现在给出\(x_1,...x_n\)，如何估计出所有的\(\mu_j,\sigma_j^2,w_i\)？
按照极大似然的过程，\(max\ \Pi_{i=1}^n[\Sigma_{j=1}^Kw_jf_j(x_i)]\)等价于\(max\ \Sigma_{i=1}^nlog(\Sigma_{j=1}^Kw_jf_j(x_i))\)。但是这样并没有对问题进行化简，所以根据前面的结果，给出一组参数\(\mu_j^{(0)},\sigma_j^{2(0)},w_i^{(0)}\)，考虑建一个表：
| 高斯分布 | \(x_1\) | \(x_2\) | .. | |
|---|---|---|---|---|
| 1 | \(\frac{w_1f_1(x_1)}{w_1f_1(x_1)+...+w_kf_k(x_1)}=w_{11}\) | \(\frac{w_1f_1(x_2)}{w_1f_1(x_2)+...}=w_{21}\) | .. | |
| 2 | \(\frac{w_2f_2(x_1)}{w_1f_1(x_1)+...}=w_{12}\) | .. | .. | |
| .. | .. | .. | .. | |
| k | .. | .. | .. |
这样就又可以进行估计了。同样，每一列的总和是1，则对于\(w_1\)有一个很自然的估计：\(w_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}}{n}\)，其余的w值同理。现在看均值和方差：
\(\mu_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}x_i}{\Sigma_{i=1}^nw_{i1}}\)，注意这里分布是对每一行的\(w_{i1}\)求和，并不是1。\(\sigma_1^{(1)}=\frac{\Sigma_{i=1}^nw_{i1}(x_i-\mu_1^{(1)})^2}{\Sigma_{i=1}^nw_{i1}}\)
迭代若干次以后就会得到比较合适的高斯分布了。
对于二维高斯分布，其参数也都可以用EM算法进行估计。
需要注意的是，这里的\(K\)是超参数，需要人为给定。