在统计学习中，我们去训练一个模型的时候，可以通过求解似然函数的方法来进行推断。而在贝叶斯方法中，我们可以计算后验概率，然后通过最大后验概率（MAP）等方法进行推断。
EM Algorithm
我们知道在最小二乘问题中，对似然函数的求解比较方便，因为其对于我们要求的参数是二次的，所以可以方便求得解析解。然而，在很多情况下，对似然函数的直接求解相当困难，而EM是简化求解似然函数的一大利器，让我们进行一些简单的讨论。
Jensen's inequality
如果自变量x是一维的，我们说是凸函数当且仅当
。如果x是多维的，是多维的对应于hessian矩阵， H是半正定的，及 。而严格的凸函数则分别对应于 或 。
Jesen inequalilty : 如果 是一个凸函数，X是一个随机变量，那么 。如果 是一个严格的凸函数，那么当且仅当。下面是一个示意图，其中X是一个随机变量满足
当然，对于凹函数的情况则完全类似，只是把上面的等号或者不等号换了方向而已。
The EM algorithm
假设训练数据集表示为 ，那么似然函数可以表示为 。其中 为我们要估计的参数， 表示隐变量。因为直接通过 估计参数往往是非常困难的，这时EM借助隐变量采取两步走的策略就有效的多。在E-step中，我们估计似然函数 的一个lower-bound,然后在M-step中再去优化这个lower-bound。在第i个样本点处，假设其隐变量的分布为 ，那么
最后一步中不等式的成立使用了Jensen's inequality。首先 是一个严格的凹函数。 为 关于 的分布 的期望。根据Jensen's inequality,
因此，对于隐变量z任何形式的分布 ，公式（6.1）给出了一个估计似然函数 的下界。的选择多种多样，假定我们当前的参数 是固定的，当然我们希望选择的 使得这个下界 和似然越接近越好。根据Jensen's inequality，要使和相等，那么 。（6.2）
当 时，其中c为某个不依赖与 的常数，明显（6.2）成立。因此只要我们选择 ，使得 ，就可以使和相等。又因为 ，所以
因此，这里 可选择为给定了样本和参数 的后验概率。现在给定了似然函数的下限，我们在M-step中要做的就是选择新的参数最大化，然后不断地迭代直至收敛。EM算法可表示为
Repeat until convergence
{
E Step ： 对于每个样本点i，计算
M Step :
}
假设和分别为第t次和第(t+1)次迭代中的参数，我们来证明。在第t次迭代中，，因此根据公式（6.1）我们有
而的选择是要使上面式子右边的项最大化，因此
上面第一个不等式成立是因为对于任意的，都成立。因此在每次EM的迭代过程中，似然函数都会递增直至收敛。
Mixture of Gaussians
假如我们要进行密度估计，而我们对数据进行统计，比如画直方图，发现数据的密度直方图存在不只一个波峰，这时假设数据是高斯分布的就不太合适，下面是一个示意图。横轴表示X的取值，纵轴表示密度。
这时，我们可以使用混合模型，比如假设数据是由M个相互独立的高斯分布组成的
其中表示第m个高斯分布，表示属于第m个高斯分布的概率。因此似然函数可以表示为
显然，直接对（6.3）求最大值是很困难的，这里我们可以用EM的来求解。
在E-step中，隐变量表示第i个样本点属于哪一个高斯函数，其属于第m个高斯函数的概率为
其中，。
在M-step中，这时lower bound 可表示为
对上面的式子求导可解的
所以，在高斯混合过程的求解中，首先对参数做一个初始的估计，然后迭代的重复E-step和M-step直至收敛。
Sampling
前面提到在贝叶斯的推断和决策中，我们需要根据后验概率和损失函数进行推断和决策。在很多时候，后验分布的形式可能非常复杂，这个时候寻找其中的最大后验估计或者在后验概率进行积分等计算往往非常困难，这时我们可以通过采样的方法来求解。
Monte Carlo Integration
Monte Carlo 积分是对积分运算的一种近似求解。假设我们要计算一个积分 ，我们知道计算机能够处理的都是离散的数，要进行真正的积分运算是不可能的。这时我们需要斗转星移，想方设法把积分转化为离散的求和运算。我们可以把 分解为两项 ，其中 是某种形式的函数，而重要的是 可以表示为某种概率分布。这时对 求积分就可以表示为求 在密度函数为 的某种分布下的期望。可表示如下
当然，这个分布函数的形式可能很复杂，无法直接求解，这时我们可以通过采样的方法去近似满足这个分布。那么，积分可表示如下
在贝叶斯分析中，Monte Carlo 积分运算常常被用来在对后验概率进行积分时做近似估计。比如我们要计算
则可以通过下面的近似运算。
其中每一个 表示从x的分布中采样出来的点。通过Monte Carlo积分计算出来的Monte Carlo标准差可表示如下
Basic Sampling
可以看到，在用Monte Carlo积分近似求解的时候，最重要的一点是从特定的分布中去采样。先来看看什么是采样。在采样的时候，我们对特定分布的了解一般有两种形式。第一种是知道了分布的具体形式，比如其密度函数或者概率函数；第二种是我们不知道密度函数 ，但是我们知道 ，其中 ，其中 是一个常数，但是这个常数可能非常难求解，这时其实 ，采样的时候不受影响。这二种情况常常发生在贝叶斯推断中，比如 ，其中 容易求解，而积分运算 则可能非常难以计算。这时 就相当于，就相当于。Anyway，对于变量x取值空间中的任意一点 ，我们都可以轻易的计算 （第一情况）或者 （第二种情况）。
知道了分布的形式，在Basic Sampling中我们可以通过一些转化的方法来采样。首先我们知道在计算机中早已实现了产生随机数的方法，其实这个随机数并不是真正的随机数，它也有一定的周期的，anyway，在实际应用中这个影响不大，我们就天真的认为这是真正的随机数。我们利用其产生（0,1）上面的随机数，显而易见，这就是对 上面均匀分布的采样。在采样中，这种（0，1）上面的随机数是最最基本的采样，我们来看看对一些比较简单的分布采样是怎样从这个最最基本的采样转为过来的。
假设我们有两个随机变量 ，并且满足，其反函数为 。假设 的密度函数为 ， 的密度函数为 ，因为 和 是相对应的，所以我们可以得到 。其中 分别表示 和 对应的累积概率函数。那么 ，对两边对 求导可以得到
因为如果 是均匀分布，那么 ，所以
假定我们现在知道了某个特定 的分布，给定了上面的关系，我们就可以从（0,1）上产生的随机数转化为符合这个特定分布的采样。假定是满足分布产生的随机数，那么 就是满足特定分布 的采样。举个例子，假设我们要对指数分布进行采样 。那么通过上面的式子我们可以计算得到 。那么 。
多元的情况也是类似的，和 的关系可表示为下面的式子
Basic Sampling的方法虽然很有效，但是其实用范围比较小，只适合于一些简单的分布，因为如果分布比较复杂就难以实施了。因为从上面看到，在计算过程中我们需要求积分和反函数，想想我们采样的目的往往就是把积分运算离散化，所以这种方法很有局限性。
Rejection sampling
相比于Basic Sampling，这种方法允许我们对更加复杂的分布进行采样。为了实施Rejection sampling,我们需要定义一个称之为Proposal Distribution的分布 ，这个分布 一般都具有较简单的形式，我们可以用前面的Basic Sampling等方法对其进行采样，下面我们来看看怎样把对 的采样和我们的目标联系起来。假设这里我们知道 ，前面提到过对 或 采样是无差别的。我们需要引进一个常数 使得对于任意的，其满足下面的性质。这时每当我们从中采样一个点 ，我们计算 ，然后从均与分布 中采样一个点 ，如果 ，那么这个样本点 就被拒绝，否则就被保留。下面是一个Rejection Sampling的示意图。图中蓝线表示 ，红线表示 ，灰色区域表示拒绝样本点 的区域。
从上面的描述中可以看出，我们最终接受的样本点其分布正好满足 。下面来看看所有从 中采样的样本被接受的概率为
可以看到，我们样本的接受率直接受的影响，所以一方面我们要保证比较大，使得，另一方面我们又要使得尽量小，因为很大的使得接受率非常低，这样我们采样的效率就非常低，耗时又耗力。因此，选择合适的 和 对rejection sampling来说就非常重要。
Importance sampling
前面提到，我们采样的目的很多时候都是为了近似积分运算，前面的采样方法都是先对分布进行采样，然后我们再用采样的结果近似计算积分。Importance sampling则两步并做一步，直接的近似计算积分。
要直接近似的计算积分，一种很自然的想法就是把的取值空间均匀地划分为很多的大小相等的小格子，假设划分的小格子共有L个，然后计算在每个小格子中随机选一个样本点 或者格子最中心的点 ，然后通过下面的式子计算近似值
这种方法在低维空间中还比较使用，但是在高维空间中却并不奏效。因为我们知道在高维空间中样本点是非常稀疏的，这时很多的格子中对应的 其实都是0，这时我们就需要浪费很多的精力在这些对计算结果没有什么帮助的地方。因此相对于这种均匀的选取样本点，我们更加倾向于去选择那是使得比较大的那些样本点。同前面的rejection sampling一样，我们需要一个能够轻易采样的proposal distribution ，这样每次从 采样，我们可以近似计算积分如下
其中 被称之为importance weights，其纠正了从中采样时对最终计算结果带来的偏差。同rejection Sampling不同的是，这里所有采样的样本点都被保留了下来。
当我们只知道 的时候情况也是类似的，这时
其中 ，
因此，
下面是一个importance sampling的示意图，其中蓝线表示 ，红线表示 ，绿线表示 。
同rejection sampling类似，选择一个合适的 也非常重要。最好能够保证通过 采样的那些样本点能够集中到 比较大的那些取值空间，不然的话也可能像前面提到的均匀选取样本点那样效率非常低下，而且这时得到的结果误差也可能很大。
Markov Chain Monte Carlo
MCMC绝对是sampling方法中大杀器，其最早是使物理学家们提出来的一种方法，后来貌似在90年代以后逐渐收到统计学家和做machine learning的热捧，成为sampling中当之无愧的大明星。MCMC构造markov chain，使其稳态分布等于我们要采样的分布 ，这样我们就可以通过markov chain来进行采样。
MARKOV CHAINS
首先来看看什么是markov chain及它的一些性质。先假设我们用一些跟时间相关的随机变量 ，每一个 对应于在时间时刻的随机变量。这些随机变量的取值空间都是一样的，先假设它们取有限个离散的值，这时的取值空间称之为state space, 里面的每一个取值称之为state 。我们说这些随机变量是一阶的Markov process，如果他们满足下面的式子
即下一个时刻的状态只跟当前的状态和转移概率有关系，跟之前的状态都没有关系。我们说一条一阶的markov chain就是一阶的markov process产生的随机变量序列 。我们这里讨论的markov chain都是一阶的，所以后面就不再提及。如果转移概率对任意时刻的 都是一样的，我们可以说这条markov chain是homogeneous的。我们用下面的式子来表示每一步中从状态到状态的转移概率。这里的一步是指都时间过渡到下一时刻。
Markov chain在时刻t处于状态的概率可以表示为 。
我们用向量 来表示在时刻各个状态的概率向量。在初始时刻，我们需要选取一个特定的，通常情况下我们可以使向量中一个元素为1，其他元素均为0，即从某个特定的状态开始。随着时间的推移，状态会通过转移矩阵，向各个状态扩散。
chain在时刻处于状态的概率 可以通过的状态概率和转移概率来求的，可通过下面的称之为Chapman-Kolomogrov equation来得到
用转移矩阵写成矩阵的形式如下
其中转移矩阵 表示 。因此 。
定义 为 ，可以看出，表示矩阵 中第个元素。
下面来看看markov chain的一些性质。
我们说一条markov chain是irreducibile当其满足下面的条件：对于任意的两个状态 都存在 使得 。及从任意一个状态 通过有限步之后都可以到达另外任意一个状态 。比如下面的这个图从状态D永远到不了状态A和B。所以这个markov chain不是irreducibile。
Figture (6.1)
而我们可以看到下面的图是irreducibile的，因为从任意一个状态都可以到达另外任意一个状态。
我们说一条markov chain是aperiodic当从一个状态 到另一个状态 所需要的步长的最小公约数为1。下图是一个不满足aperiodic的例子，因为从状态A到状态E所需的步长只能是
Figture(6.2)
一条markov chain有一个平稳分布 ，是指给定任意一个初始分布 ，经过有限步之后，最终都会收敛到平稳分布。平稳分布具有性质 。换言之，是转移矩阵特征值 所对应的左特征向量。
一条markov chain如果是irreducible和aperiodic的，那么这条markov chain就具有平稳分布。想一想，如果这条markov chain不满足irreducible，比如下面的图中，如果初始状态为 ,那么最终会一直保持状态。如果初始状态为 ，那么最终会在 中徘徊，所以明显不满足平稳分布的条件。
另外一种不满足irreducible情况是在图（6.1）中，无论从哪个状态开始，最终都会在徘徊，这样的平稳分布是不好的，其中 ，因为我们希望平稳分布中每一个状态的概率都大于0。
如果如果这条markov chain不满足aperiodic，比如在图(6.2)中，我们设置初始状态为 ，那么状态 只能在步数为奇数的时候被访问， 只能在步数为偶数的时候被访问，那么明显最终不能达到稳态分布。可以证明，如果转移矩阵 没有-1的特征值，那么这条markov chain就是aperiodic的。
一条markov chain拥有平稳分布的一个充分条件是对于任意两个状态和， 其满足detailed distance ： 。
可以看到，此时 。
所以，明显满足平稳分布的条件。
如果一条markov chain满足detailed distance，我们就说它是reversible的。
在markov chain中，对于随机变量的取值是连续的情况，上面的这些定义和性质都是类似的。比如这时转移概率为 ，那么。
Chapman-Kologronvo equation这时可以表示为 。
平稳分布的性质表示为。
Metropolis-Hastings algorithm
先来看看Metropolis algorithm，其是MCMC的一个初级版本。现在假设我们要从 中采样，Metropolis algorithm可表示如下
如果 ，即 ，说明这个地方密度比较大，可以多在这儿采样，这时我们就必然接受这个候选样本点。如果 ，即，说明这个地方密度相对来说不怎么大，我们就以的概率接受这个候选点。
我们来看看Metropolis algorithm的性质，我们可以看出这个markov chain依次产生 。从算法中可以轻易看出， 只依赖于 ，所以是一阶的markov chain。经过有限的步之后，这条chain就会处于平稳分布 ，所以我们从 中采样相当于从 中采样。
Metropolis-Hastings algorithm
Hastings对Metropolis algorithm稍作改进，得到了所谓的Metropolis-Hastings algorithm。在Metropolis algorithm中，我们要求proposal distribution 必须是对称的，这里我们对其放宽了条件， 不满足对称性也是可以的。这时对候选样本点的接受概率就可表示为
可以看到，当 是对称的时候，Metropolis algorithm和Metropolis-Hastings algorithm是一样的。
举个例子，我们要对分布 进行采样，如果使用，即每次 是从均匀分布 采样的一个样本点，下面是示意图，这里 。横轴表示步长，纵轴表示采样点的值。
Figure (6.3)
注意到图中有很多阶段曲线是非常平稳的，这时我们说这条chain是poorly mixing。相对的如果使用作为 proposal distribution，那么得到示意图如下。这时我们说这条chain是well mixing.
Figure （6.4）
Metropolis-Hasting Sampling as a Markov Chain
在Metropolis-Hasting Sampling中，我们看到转移概率可以表示为
为了要证明Metropolis-Hasting Sampling中，我们产生的markov chain存在平稳分布并且其平稳分布就是我们要采样的分布 。我们先来证明其满足detailed balance, 因为这是其平稳分布的一个充分条件。
对于任意的一对
所以，对已任意的 我们都有，因此markov chain的平稳的，并且其平稳分布为 。
一般在实施Metropolis-Hasting Sampling或者其他的MCMC sampling的时候，前面的1000至5000个采样点都会被扔弃，然后一些测试收敛的test会被用来观测是否已达到平稳分布。因为如果初始样本点选择在中一个很尖的波峰，那么后面很长一段时间内的候选采样点都很可能被拒绝，之后才慢慢进入到平稳分布中。在达到平稳分布之前的时间被称之为burn-in period。一个糟糕的初始采样值 或者一个糟糕的proposal distribution都可能使burn-in period的时间很长，严重影响采样的效率。其中对初始采样点 的建议是尽量选取那些在分布中心的样本点。
一条chain被认为是poorly mixing的，如果它在某一段时间内都停留在取值空间的某个小区域，反之则认为是well mixing。可对照图（6.3）和（6.4）。当一个分布存在多峰的情形时，很可能会发生poorly mixing的情况，这时一般由两种方法去处理。第一种是使用多条链来采样，每条链采用不同的初始值。另外一种是使用Simulated Annealing的方法，下面简单介绍一下这种方法。
Simulated Annealing
如果我们的目标是在一个复杂的函数上求最大值，一种常用的贪婪式算法称之为爬山算法，也是就所谓的梯度方法。但是这种方法有一个很大的缺点，当函数表示的曲线或者曲面存在多个波峰的时候，这时找到的局部最大值很可能不是全局最大值，而且差距还可能非常之大。
Simulated Annealing来自冶金学的专有名词退火。退火是将材料加热后再经特定速率冷却，目的是增大晶粒的体积，并且减少晶格中的缺陷。材料中的原子原来会停留在使内能有局部最小值的位置，加热使能量变大，原子会离开原来位置，而随机在其他位置中移动。退火冷却时速度较慢，使得原子有较多可能可以找到内能比原先更低的位置。Simulated Annealing的原理也和金属退火的原理近似：将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有"能量"，以表示该点对命题的合适程度，在我们采样的时候就表示为对该候选样本点的接受程度。算法先以搜寻空间内一个任意点作起始：每一步先选择一个"邻居"，然后再计算从现有位置到达"邻居"的概率。可以证明，Simulated Annealing算法所得解依概率收敛到全局最优解。
利用Simulated Annealing，我们只需要把Metropolis sampling中的 修改如下
Choosing a Jumping (Proposal) Distribution
在Metropolis Sampling中， 必然满足对称性原则，而Hasting-Metropolis中则不必。选取proposal distribution一般有两种方法，第一钟方法是Random Walk。新的候选采样的 是当前采样点 加上一个随机变量。即 。所以 ，其中 是随机变量的密度函数。如果 ，明显这时是对称的，可用于Metropolis Sampling。第二种方法是independent chain，即 。这时候选采样点的选择独立于当前采样点的值。注意此时明显 不是对称的，所以此时只能用于metropolis Hasting sampling。
我们可以通过调节 的variance或者covariance matrix的特征值来调节采样的接受率。如果variance变小，就会提高接受率，反之亦然。如果variance很大，那么采样在取值空间中移动的比较大，但这个时候的接受率很低，所以容易出现poorly mixing。但是如果variance很小，那么在取值空间移动的就很小，这时接受程度又比较高，也容易出现poorly mixing。所以要掌握好variance的tradeoff。
Convergence Diagonistics
我们可以看到，在MCMC采样的过程中，序列中的每个样本都是跟其前面一个样本强烈相关联的。因为我们希望我们采样得到的样本点是独立同分布的，这时我们可以对采样的序列进行处理，比如在burn-in period之后，我们每隔个样本点选取一个样本点最为我们最终的采样。
可以进行一些测试来验证采样是否已达到平稳。其中一种叫做Geweke test的方法是在去除 burn-in period之前的那些采样点后，把这些采样点的前10%和后50%拿出来作为两个数据集，然后用假设检验的方法比较这两个数据集的均值是否相等，如果相等的话就说明确实达到了平稳。
在使用MCMC采样的时候，可以使用一条很长的chain或者几条不同初始值的chain进行采样，在并行环境下，使用多条chain可能更有优势，但是使用一条长的chain的效果很可能会更好，因为burn-in period的时间可能会非常之久。
The Gibbs Sampler
Gibbs Sampling是Metropolis Hastings Sampling的一个特例。在Gibbs Sampling中， ，即每个候选采样点都会被保留。Gibbs Sampling的主要特性是其针对条件概率采样，如果我们的取值空间是多维的，在联合密度函数下进行采样通常是比较困难的，而Gibbs Sampling每次采样只针对一个随机变量时，其他随机变量都已经取好了某个特定的值。Gibbs Sampling的过程可以表示为
Sampling 和EM的关系
当EM中隐变量的取值空间是连续的时候，在M-step中，目标函数可表示如下
这时计算这个积分我们可以用前面提到的采样的方法。
Bagging
Bagging方法其实就是用bootsrap方法来做推断。后面我们会看到用Boosting的方法往往会更加有效。对于用bootsrap选取的B个bootsrap数据集， ，用Bagging方法得到估计如下
Model Averaging and Stacking
先从贝叶斯的角度来看看Model Average。假设我们用 来表示训练数据集，我们拥有共M个子模型。是我们感兴趣的某个东东，比如说是在自变量X =x处的预测值等等。那么
所以在有很多子模型时，一个贝叶斯的预测值就是子模型预测的一个平滑，每个子模型的贡献正比于与其在训练集下的后验概率。
在Model average的策略中，最简单的就是Committee Methods，其假设每个子模型的后验概率都相等，这时的预测值就是每个子模型预测值的平均。当然稍微复杂一点，我们可以使用前面讲到过的BIC来计算每个子模型的后验概率，然后最后的预测值是它们的加权平均。
我们再从经典统计的角度来看看，假设M个子模型分别表示为 如果我们使用squared error，每个模型的权重表示为 。那么
这里预测变量x是固定的，假设训练集（包括隐变量Y）服从分布真实的联合分布P。
因此设，
那么
可以证明
所以从这一点上看，使用Model average的效果貌似肯定不会比使用单个子模型的效果差。但是我们的训练数据集的联合分布跟真实的联合分布可能并不一致，所以上面的结论并不总是成立。
Stacking
前面在做Model Average的时候并没有考虑到子模型的复杂度，如果某个子模型的复杂度很高，其对training data的拟合效果就会很好，这时在average的时候就会偏向于给这个子模型较大的权重。但其实我们知道复杂度很高的子模型往往具有很高的variance，其在新数据的预测效果上并不好。Stacking方法就通过一点改进尽量去避免这种情况。 令 表示用训练数据集中除第i个样本点外的所有样本点拟合出的第m个子模型。子模型的权重向量可以表示为
因为在计算第i个点的误差时，用的是不包含着这个点的训练集拟合的子模型，所以会尽量减小子模型overfitting的情况。最终的预测值可表示为 。如果在增加一些限制，比如权重向量中每个元素都大于等于0，并且和加起来为1，这时权重向量的每个元素可以看做是对应子模型的后验概率，只是这里的求解过程是一个quadratic programming问题。可以看出，stacking和leave-one-out cross validation有很大的联系，如果限制权重向量中只有一个元素为1，其余均为0，那么这时stacking就相当于用leave-one-out cross validation做model selection。
Stochastic Search: Bumping
Bumping其实并不是一种Model Average的方法，其倒更像是一种model selection的方法。我们使用bootstrap方法从训练集中抽出B个bootstrap数据集。对于每个Bootstrap数据集拟合出一个预测模型。我们要选择的Bootstrap数据集满足下面的条件
然后我们使用来做预测。这里我们对每个bootstrap数据集使用的子模型都是一样的，只是使用不同的训练集去训练，这样可以看出用bumping拟合出来的模型就是这个特定的子模型能够在数据集的子区域上表现最好的状态。
联系客服
