SC-复习提纲
1 Methods for Generating Random Variables
- The Inverse Transform Method
- The Acceptance-Rejection Method
- 根据一个参考的分布生成；
- 注意接受的概率徐取决于概率分布的比值 c，因此选取的参考分布的好坏决定了算法的效率
- Transformation Methods
- Sums and Mixtures：随机变量的和或者混合是一种特殊类型的变换
- 理解：sum 是直接把多个随机变量加起来；而 mixture 则是把分布函数「加起来」；具体的实现的话，就是按照不同 组分（参数不同的随机变量）的权重，从各自的分布中去抽样，组合起来即可。
- 一个例子：A Poisson–Gamma Mixture Is Negative-Binomially Distributed，参见 https://gregorygundersen.com/blog/2019/09/16/poisson-gamma-nb/
2 Methods for Generating Random Vectors
Multivariate Normal Distribution
思想是利用 MND 的性质，通过一元标准正态分布变换得到
\[ X=ZQ+J\mu' \]
所以关键是如何得到矩阵 Z 满足 \(Q'Q=\Sigma\) ，可以通过 谱分解、奇异值分解、Choleski 分解 得到
- 比较了不同采样方式的效率问题，主要的计算量在协方差矩阵的分解上。
Mixtures of Multivariate Normals
Wishart Distribution
一个好复杂的分布 \[ M=X'X \] 其中 \(X\) 是从 \(N_d(0,\Sigma)\) 中抽取的 \(n\times d\) 随机矩阵，则称 \(M\sim W_d(\Sigma, n)\)
直接通过正态分布需要生成 \(nd\) 个随机数，来决定 \(d(d+1)/2\) 个元素值；更为高效的方法是基于 Bartlett 分解，只需要生成 \(d(d+1)/2\) 个随机数（从标准正态和卡方分布中采样）；和生成 MND 一样最后需要做一些变换。
Uniform Distribution on the \(d\)−Sphere
- 即生成在 d 维球面上的均匀分布；思想是生成一个随机向量，做归一化处理；而所谓「随机方向」的向量，可以从 d 维的标准正态分布中得到（直接生成 d 个一维标准正态即可）。
3 Monte Carlo Integration
Simple Monte Carlo estimator
针对的是 \[ \theta=\int_a^bg(x)dx \] 这类积分问题，将其看作是 \((b-a)E_{X\sim U[a,b](g(X))}\) 的期望形式，于是根据 LLN 进行采样计算
由于这里引入了一个 \([a,b]\) 上的均匀分布，因此需要加上其分布函数，也就是要注意在最后的结果中乘以这个系数，即 \[ \hat\theta=(b-a)\overline{g_m(X)}=(b-a){1\over m}\sum g(X_i)\tag{3.1} \] 其中的 \(X_i\) 根据 \([a,b]\) 上的均匀分布得到
Variance and Efficiency
针对标准正态的 CDF \[ \Phi(x)=\int_{-\infty}^x{1\over\sqrt{2\pi}}e^{-{t^2\over 2}}dx \] 的 MC 估计问题，给出了两种方法：
一种是利用之前的 SMC 去估计 \(\Phi(x)-\Phi(0)\) （区间有界）；
（Hit-or-miss Monte Carlo）另一种是利用 \[ \Phi(x)=P(Z\le x)=EI(Z\le x) \] 将此概率转化为示性函数的期望估计问题，也可以用 MC 方法求解
两种方法均为无偏估计，那么就需要比较那一个的方差比较小（哪一个更有效率，仅需要较少的采样即可使得方差较小）
事实上，上述两种方法的方差分别为 \[ {(b-a)^2\over m}Var(g(X))\\ {F(x)(1-F(x))\over m} \]
Variance Reduction
下面介绍了集中方差缩减的方法
Antithetic Variables
对偶变量法
想法是利用两个随机变量之和的方差：若这两个随机变量是负相关的，则他们之和的方差要小于各自的方差之和（从而达到方差缩减的效果）
利用到 MC 积分上来。我们可以把 \([a,b]\) 上的随机分布通过变量替换全部转换到 \([0,1]\) 上的均匀分布；注意到 \(U\) 和 \(1-U\) 同分布，如果我们把全部的采样数 m 中的一半生成标准均匀分布，另一半用 \(1-u_j\) ，如果我们能过说明最终生成的样本对 \(g(U_i), g(1-U_i)\) 之间是负相关的，那么就可以达到 VR 的效果
事实上我们有更一般的结论：定理，设 \(X=(X_1, ...,X_n)\) 各分量相互独立，\(f,h\) 是同方向的单调函数，则 \[ E[f(X)h(X)]\ge E[f(X)]E(h(X)) \] 即 \(f(X),h(X)\) 正相关。
应用到这里，我们取 \(f(X)=g(U_i), h(X)=-g(1-U_i)\) （第 i 个样本仅用所生成的第 i 个标准均匀分布），则 \(g(U_i), g(1-U_i)\) 负相关。
5 Control Variates
这里的方法是引入控制变量，假设要估计的量为 \(\theta=E[g(X)]\)，而另一个函数满足 \(E[f(X)]=\mu\)，则 \[ \hat \theta_c=g(X)+c(f(X)-\mu) \] 仍为无偏估计，其方差为 \[ Var(\hat\theta_c)=Var(g(X))+c^2Var(f(X))+2cCov(g(X),f(X)) \] 上式最小是在 \[ c^*=-{Cov(g(X),g(X))\over Var(f(X))} \] 处取得，最小值为 \[ Var(\hat\theta_{c^*})=Var(g(X))-{[Cov(g(X),f(X))]^2\over Var(f(X))} \] 显然，方差缩减为 \[ {Var(\hat\theta)-Var(\hat\theta_{c^*})\over Var(\hat\theta)}=[Corr(g(X),f(X))]^2 \]
可知，当 \(f,g\) 强相关时，方差缩减是有效的，否则不会导致方差减少。
Antithetic variate as control variate
对偶变量的方法可以看成是控制变量方差的特例
一般地，若 \(\hat\theta_1, \hat\theta_2\) 为两个无偏估计，则 \[ \hat\theta_c=c\hat\theta_1+(1-c)\hat\theta_2 \] 也是无偏估计，若 \(\hat\theta_1,\hat\theta_2\) 同分布，计算可得最优值为 \(c^*={1\over2}\) ，即 \[ \hat\theta_{c^*}={\hat\theta_1+\hat\theta_2\over2} \]
显然这也是（这种特定选择下）对偶变量方法下的估计量
Several control variates
- 可以将控制变量的方法扩展到采用多个控制变量的场合
- 此时，控制变量方法下的估计量 \(\hat\theta_{\hat c}\) 以及最优的 \(c_i\) 可以通过回归模型来估计
Importance Sampling
这部分讲义上有点乱，参考 https://statweb.stanford.edu/~owen/mc/Ch-var-is.pdf；主要的问题，或许是，之前的 SMC 仅能解决 \(\int g(x)dx\) 这样的积分问题，而下面的 IS 虽然是由这一问题引入的，但是其能解决更为普遍的 \(E[g(X)]=\int f(x)g(x)dx\) 的形式，因此出现了混乱。
在前面内容中,我们将有限区间上的积分视为是此区间上的均匀分布随机变量的某个函数的期望值. 这种方法的缺点是不能直接用于无穷积分的估计, 以及当被积函数在积分区间上不是很均匀的时候 效率很低. 那么很自然地可以考虑均匀分布以外的其他加权函数. 这就导致“重要性抽样方法”.
假设有一个分布 \(f(x)\) ，其在 \(g(x)>0\) 的区间上值也 \(>0\) （其他地方密度函数为 0），我们记 \(Y=g(X)/f(X)\) ，则简单的积分可化为 \[ \int g(x)dx=\int{g(x)\over f(x)}dx=EY \] 的期望形式，于是可以通过 SMC 方法来估计 \(EY\) \[ {1\over m}\sum Y_i={1\over m}\sum {g(X_i)\over f(X_i)} \] 其中 \(X_i\) 为从 \(f\) 中抽取的样本，将 \(f\) 称为 重要函数。
显然估计的方差取决于 \[ Y={g(X)\over f(X)}\tag{6.1} \] 的方差，也就是说，我们所选取的重要函数 \(f\) 约接近被积函数 \(g\) ，估计的方差越小。
更为一般的，我们的目标是求 \[ \theta=E[g(X)]= \int_Ag(x)p(x)dx \] 的形式，注意在上面我们是求一个积分，而这里是求随机变量 \(X\) 的函数的期望（当然也是积分，不过积分的形式改变了）；不同于之前的 \(f\) 是我们选定的「重要函数」，这里的 \(p(x)\) 是随机变量 \(X\) 的 pdf【教案中仍然用 f 表示 X 的分布，这里采用上面链接中的符号记其分布为 p】
\(p(x)\) 不好 sampling，引入 重要性抽样函数/包络函数 envelope \(q(x)\)，则 \[ \theta=\int_Ag(x){p(x)\over q(x)}q(x)dx \] 于是我们可以借由 \(q(x)\) 来采样，\(\theta\) 的估计量为 \[ \hat\theta={1\over n}\sum g(X_i){p(X_i)\over q(X_i)} \]
和（6.1）式类似，我们记 \(Y=g(X){p(X)\over q(X)}\) ；类似的，当我们所选取的重要性抽样函数 \(q(x)\) 越接近 \(g(x)p(x)\) ，估计的方差就越小。
可以证明 \(E[Y]=\theta\) ，则其方差为
\[ Var(Y)=E[Y^2]-E[Y]^2=\int {(g(x)p(x))^2\over q(x)}dx-\theta^2 \]
7 Stratified Sampling and Stratified Importance Sampling
Stratified Sampling
相较于之前的 SMC 积分中，通过变量替换转换为在 \([0,1]\) 上积分，之后通过采样标准均匀分布的方法，我们可以利用积分的线性性，将整个积分转化为多个小区间上的积分之和，也即，将积分区间 \([0,1]\) 分成若干个自区间，分别在子区间上采样并估计，最终的估计为 \[ \hat\theta^S={1\over k}\sum_{i=1}^k\hat\theta_i \] 这里前面的系数又是一个符号问题：按照积分的线性性本来是简单的求和，事实上这里的 \({1\over k}\) 正是（3.1）式中的系数 \(b-a\)
可以证明，基于分层抽样法得到的估计的方差要小于 SMC
Stratified Importance Sampling
- 分层抽样的思想也可以结合到重要性采样的方法中
8 Monte Carlo Methods for Estimation
这一部分介绍了使用 MC 方法对于一个统计量的 标准差（标准误差）、MSE，以及 CI 的估计；也即，MC 方法在统计检验中的作用。
Monte Carlo Methods for Estimation
Monte Carlo Estimation and Standard Error
我们采用「plug-in」方法。
例如，对于均值的估计 \(\bar X\)，我们要估计其标准差 \[ se(\theta)=\sqrt{Var(X)/m}\\ \] 而对于方差的「plug-in」估计为 \[ \hat\sigma^2=\hat Var(X)={1\over m}\sum(x_i-\bar x)^2 \] 于是标准差的估计量为 \[ \hat{se}(\bar X)={1\over m}\sqrt{\sum(x_i-\bar x)^2}\tag{8.1} \]
当然，也可以采用无偏估计量，即将前面的系数改为 \({1\over\sqrt{m(m-1)}}\)
Estimation of MSE
事实上这里很迷惑，想了半天感觉和前面的 SE 的估计是完全不一样的东西，不知道为什么放在了一起讲。
在上面估计统计量的 SE 的过程我们在普通的统计课程中就有讲到，我们可以直接通过「plug-in」的思想，利用所观测到的数据计算得到；相较而言，MSE 的估计（或者 bias）则更为困难，下面介绍了对于 MSE 的 MC 估计方法，注意到这里需要反复抽取 n 组大小为 m 的样本（而上面我们则直接用了一次取样的 m 个样本）。
这里原因在于，在 \(SE(\hat\theta)\) 的估计中，我们将其转化到了对于 \(Var(X)\) 的估计，而 \(MSE(\hat\theta)\) 中，我们无法化简掉未知的参数 \(\theta\) ，因此在这里 \(\theta\) 是已知的。这是和上面的对于 SE 的估计完全不同的地方，事实上若我们不知道参数的真实值我们是完全无法进行估计的。因此，为了估计这一期望，我们从原始分布中取了多组观测（因此事实上我们模拟的是 \(\hat\theta\) 的分布），然后计算这一期望值。
具体来说，均方误差为 \[ MSE(\hat\theta)=E[(\hat\theta-\theta)^2] \] 我们从 总体 （注意 \(\theta\) 已知）中抽取 n 组大小为 m 的样本，分别用这些样本计算估计量 \(\hat\theta^{(j)}\) ，于是 \(\hat\theta\) 的 MSE 的 MC 估计为 \[ \hat{MSE(\hat\theta)}={1\over m}\sum(\hat\theta^{(j)}-\theta)^2\tag{8.2} \]
这样子看下来好像 MSE 的估计没什么意义？下面给出了一个 例子：当样本存在异常值的时候，我们会采用截尾分布（即删去数据中最大和最小的那些数据点）；通过模拟实验，我们发现
均值的稳健估计(截尾均值估计)在总体分布被污染时能降低 MSE
Estimating a confidence level
在应用统计中经常遇到的一个问题是估计总体的分布. 比如, 许多常用的统计推断 方法和工具都是基于正态性假设下的. 而在实际中, 总体分布非正态是经常的, 估计 量的分布可能不知道或者没有显示表示, 此时, Monte Carlo 方法则可以用来进行 统计推断.
这里其实和前面的 MSE 一样，我们都是已知样本点的分布的，这里做的意义，不是为了进行统计检验，而是通过这些实验，来判断一个统计检验的好坏（？）
例子：我们假定数据分布为正态的情况下，对于其方差可以借由样本方差的分布得到检验统计量，最终得到一个区间估计；然而，若假设不满足，那么这样做的结果可能就会造成很大的误差；通过 MC 方法，我们可以度量这一可能的误差。对于一个分布 \(X\sim F_X\) （不满足正态分布），我们进行 m 次大小为 n 的采样，对于这 m 组样本仍然采用正态假设下的 CI 计算公式，统计总体的真实方差落在计算到的 CI 中的概率，即得到 经验置信水平。
前面的例子中考虑的问题都是在总体分布已知的情形下, 对参数进行Monte Carlo估计. 因此这种 情况下的Monte Carlo方法也称为是参数Bootstrop 方法. 有关于Bootstrap 方法我们将在后面学习.
9 Monte Carlo Methods for Hypothesis Testing
Empirical Type I error rate
和之前计算 CI 的 经验置信水平 的方法一样，我们同样可以估计一个 HT 的 经验一类错误率
给了一个 例子 说明这种估计的意义：正态分布的偏度检验 （Wald 检验）。利用 CLT 可知偏度系数的估计量 \(b_1\) 的渐进分布 \[ \sqrt{n}(b_1-\beta_1)\to N(0,6) \] 于是可以基于此构建偏度系数是否为 0 的 HT；
然而，当 n 比较小的时候，这种近似可能存在问题，我们计算了不同 n 下的 经验一类错误率，发现 n 较小时经验一类错误率达不到预期的显著性水平如 \(\alpha=0.05\)
事实上，可以求出 \[ Var(b_1)={6(n-2)\over(n+1)(n+3)} \] 将这样的正态近似可能取得更好的 经验一类错误率。
Power of a Test
- 回忆检验的 功效 的定义：当 \(\theta\in\Theta_1\) ，即备择假设成立的情况下，拒绝原假设的概率；
- 和前面估计一类错误率不同之处在于，这里的备择假设是一个区间而非单个点；因此，事实上我们可以画出 经验功效函数曲线
- 具体来说，对于一个特定的 \(\theta_1\in\Theta_1\) ，利用 MC 方法估计在这种情况下的检验的功效，把点连起来构成 经验功效函数曲线。
Power Comparisons
- 终于有一个比较好理解的应用了：基于 MC 方法估计得到的经验功效函数曲线，我们可以比较不同检验方法的优劣。
- 例如，在书中给出了三种不同的 正态性检验 （Shapiro-Wilk检验, Energy检验 和偏度检验）的功效曲线，从而给出哪种检验可以保证较高的 power 的经验估计。
10 The Bootstrap
Efron 在1979,1981和1982年的工作中引入和进一步发展了Bootstarp方法, 此后发表了大量的关于 此方法的研究.
Bootstrap方法是一类非参数Monte Carlo方法, 其通过再抽样对总体分布进行估计. 再抽样方法将 观测到的样本视为一个有限总体, 从中进行随机(再)抽样来估计总体的特征以及对抽样总体作出统计推断. 当目标总体分布没有指定时, Bootstrap方法经常被使用, 此时, 样本是唯一已有的信息.
Bootstrap 一词可以指非参数Bootstrap, 也可以指参数Bootstrap(上一讲中). 参数Bootstrap是指总体分布 完全已知, 利用Monte Carlo方法从此总体中抽样进行统计推断; 而非参数Bootstrap是指总体分布完全未知, 利用再抽样 方法从样本中(再)抽样进行统计推断.
可以视样本所表示的有限总体的分布为一个“伪”总体, 其具有和真实总体类似的特征. 通过从此“伪”总体中 重复(再)抽样, 可以据此估计统计量的抽样分布. 统计量的一些性质, 如偏差, 标准差等也可以通过再抽样来估计.
Bootstrap Estimation of Standard Error
Bootstrap Estimation of Bias
这一部分在之前的文章 https://www.cnblogs.com/easonshi/p/12122470.html 中写得比较详细了，参看。这里的关键，在于理解 \(\theta=t(F)\)，以及 「plug-in」 方法的操作/推导过程。
11 Bootstrap Confidence Intervals
- The Standard Normal Bootstrap Confidence Interval
- The Percentile Bootstrap Confidence Interval
- The Basic Bootstrap Confidence Interval
- The Bootstrap t interval
12 The Jackknife
偏差的 Jackknife 估计
标准差的 Jackknife 估计
同样在上面的那篇笔记中有总结了
另外，需要指出的是，Jackknife 方法仅在估计量是「平滑」的时候才有效（数据的小幅变化对应于估计量的小幅变化）
例如，中位数就不是平滑的，给了一个例子：用 Jackknife 方法估计从 \(1,2,...,100\) 中随机抽取的10个数的 中位数 的标准差，这个例子中 Jackknife 和 Bootstrap 估计相差很远，就是因为这一问题。
13 Introduction to Markov Chain Monte Carlo Methods
Integration problems in Bayesian inference
介绍了一个重要的应用背景：在 Bayesian 推断中，不同于我们之前讨论的 MC 积分的问题 \(E[g(X)]\) ，我们需要做的积分是后验分布 \[ E[g(\theta|x)]=\int g(\theta)f_{\theta|x}(\theta|x)d\theta={\int g(\theta)f_{x|\theta}(x|\theta)\pi(\theta)d\theta\over \int f_{x|\theta}(x|\theta)\pi(\theta)d\theta} \]
可以写成更为普遍/简单的形式 \[ E[g(Y)]={\int g(t)p(t)dt\over \int p(t)dt}\tag{13.1} \] 其中的 \(p(t)\) 为一个密度或者似然（积分不需要=1）；显然在贝叶斯统计中 \(p(t)\) 即为后验分布（可以忽略常数项）
这里的好处是，后验分布我们不需要求出前面的归一化系数，只需要知道含有参数的主要部分即可；
式（13.1）在形式上和之间的积分问题忽略掉分母的部分是一致的（当然分母很容易处理），区别之处在于，这里的后验分布可能很复杂，很难从中抽取样本。
于是就有了 MCMC，前一个 MC，马尔科夫链就是要生成这种分布，后一个 MC 则是正常的蒙特卡洛积分。
Markov Chain Monte Carlo Integration
- 先给出了一些定理
- 说明了
- 不可约+非周期+正常返的马氏链存在唯一的平稳分布
- LLN：马氏链的实值函数的遍历均值以概率 1 收敛到极限分布下的均值
- CLT：遍历均值作合适的变化依分布收敛到标准正态分布
14 The Metropolis–Hastings Algorithm
MH(Metropolis-Hastings) 算法是一类常用的构造马氏链的方法, 其包括了: Metropolis抽样, Gibbs抽样, 独立抽样, 随机游动抽样等；
这里的主要部分就是要寻找一个合适的对于马氏链的更新方式，也即从 \(X_t\) 更新到 \(X_{t+1}\) ，算法表示为
- 构造合适的 proposal 分布 \(g(.|X_t)\)
- 从 proposal 中生成 \(Y\)
- 若 \(Y\) 被接受，则 \(X_{t+1}=Y\) ，否则 \(X_{t+1}=X_t\)
其中 proposal 的选择要使得产生的马氏链的平稳分布为目标抽样分布 f，需满足的正则化条件包括 不可约、非周期、正常返 。一个具有和目标分布相同支撑集的 proposal 一般会满足这些正则化条件
Metropolis-Hastings Sampler
https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm#Formal_derivation 天这里讲得太清楚了，为什么之前的文章从来没有过从这个角度来写的
下面的这两篇可以拓展阅读
简单来说，我们生成的马氏链的过程中，转移概率由两部分组成，从 proposal 中得到，经过拒绝接受的过程，也即 \[ P(x'|x)=g(x'|x)A(x',x)\tag{14.1} \] 而为了使得生成的马氏链的平稳分布为目标分布 \(P(x)\)，我们需要满足 细致平稳 detailed balance 条件 \[ P(X'|x)P(x)=P(x|x')P(x') \] 即任意两个状态之间的转移概率是一样的，也即 \[ {P(x'|x)\over P(x|x')}={P(x')\over P(x)} \] 将式（14.1）代入，即接受概率需要满足 \[ {A(x',x)\over A(x,x')}={P(x')g(x|x')\over P(x)g(x'|x)}\tag{14.2} \] 为了满足上式，一个常用的接受概率为 \[ A(x',x)=\min\{1,{P(x')g(x|x')\over P(x)g(x'|x)}\}\tag{14.3} \] 显然 \(A(x',x), A(x,x')\) 中肯定有一个为 1，即满足式（14.2）
The Metropolis Sampler
- Metropolis sample 是 MH 算法的一种特殊形式，这里的 proposal 是对称的 \(g(x'|x)=g(x|x')\)，因此接受概率为 \(\min\{1,{P(x')\over P(x)}\}\)
Random Walk Metropolis
- 随机游动Metropolis抽样方法是Metropolis方法的一个例子。假设候选点 Y 从一个对称的 proposal \(g(Y|X_t)=g(|X_t-Y|)\) 中产生；则在每次的迭代中，仅需要从 proposal 中产生一个增量 Z，然后 \(Y=X_t+Z\)；比如增量从正态分布中产生，此时候选点
\[ Y|X_t\sim N(X_t,\sigma^2) \]
The Independence Sampler
- 是 MH 算法的另一个特例，这里的 proposal 不依赖于前一步的状态值，即 \(g(Y|X_t)=g(Y)\) ，此时的接受概率为 \[ A(Y,X_t)=\min\{1,{P(Y)g(X_t)\over P(X_t)g(Y)}\} \]
15 The Gibbs Sampler
- Gibbs 抽样是MH算法的一个特例, 其经常用于目标分布是多元分布的场合. 假设所有的一元条件分布 (每个分量对其他分量的条件分布)都是可以确定的, Gibbs抽样使用这些一元条件分布进行抽样.
- Gibbs 采样主要解决的是多元分布的问题，这里有一个很强的假设是所有的边际分布都是已知的；这样的话，也就不需要 proposal 了，并且所有生成的点都是接受的。
16 Monitoring Convergence
检查产生的马氏链是否收敛，最直接简单的就是图示了
- Convergence diagnostics plots
- Monte Carlo Error
- The Gelman-Rubin Method
17 EM Algorithm
期望-最大化(EM)算法是一种在观测到数据后, 估计未知参数的迭代优化方法. 其能非常简单地执行并且能够通过稳定,上升的步骤非常可靠地找到全局最优值. 对EM方法详细介绍请参考非常好的教材 McLachlan and Krishnan, 1997, EM algorithm and extensions, Wiley.
https://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm#Proof_of_correctness wiki 真的太赞了，符号表达要比之前的笔记 https://www.cnblogs.com/easonshi/p/12113808.html 清楚无数倍，有机会重写一些
基本的想法就是，在 MLE 中，由于存在一些 删失数据/隐变量，我们无法得到具体的 marginal likelihood 的形式，因此需要通过 complete likelihood 积分将其消去 \[ L(\theta;X)=p(X|\theta)=\int p(X,Z|\theta)dZ\tag{17.1} \] 然而，我们无法得到 Z 的分布，EM 的思想就是利用现有的的对于参数的估计来近似 Z 的分布，即利用了 \(Z|X,\theta^{(t)}\)。
Expectation step (E step)：在当前的估计 \(\theta^{(t)}\) 下，得到（17.1）中的 marginal likelihood 的估计的形式，即 \[ Q(\theta|\theta^{(t)})=E_{Z|X,\theta^{(t)}}\log L(\theta;X,Z)\tag{17.2} \] 注意这里我们需要得到的是 Q 的形式，其是未知参数 \(\theta\) 的似然函数的估计
Maximization step (M step)：通过上述的公式更新我们对于 \(\theta\) 的估计 \[ \theta^{(t+1)}=\arg\max_\theta Q(\theta|\theta^{(t)})\tag{17.3} \]
迭代求解
这是算法的基本框架，上面的 wiki 链接中给出了算法有效性的证明，其核心的结论就是：我们通过最大化 Q 得到的更新值 \(\theta^{(t+1)}\) ，其所对应的 marginal 似然也是更好的，即 \[ \log p(X|\theta^{(t+1)})\ge \log p(X|\theta^{(t)}) \]
