memetic算法代码

文章目录

• 浅析：从最大似然估计（MLE）、最大后验估计（MAP）到EM算法

• 1. 最大似然估计算法（Maximum Likelihood Estimation, MLE）
• 2. 最大后验估计（Maximum A Posteriori, MAP）
• 3. 举例区分MLE和MAP
• 4. 期望最大算法（Expectation Maximum, EM）

浅析：从最大似然估计（MLE）、最大后验估计（MAP）到EM算法

最近阅读推荐系统领域的一些文章后，感觉最深的是在采样是对概率论方面理论知识的要求颇高。且建立优化目标时，广泛应用MLE、MAP、EM等优化方法。本文在参考一些优秀的博客后简单给出自己对于三种优化方法的粗鄙理解。

1. 最大似然估计算法（Maximum Likelihood Estimation, MLE）

在我看来MLE本质上就是一种用频率估计概率的思想。

先来举个简单的例子，我们小学时候都学过Buffon的投针实验，当时是为了估算圆周率 $\pi$。其基本思想可以简化为有边长为1的正方形内切一个单位元，随机投针在正方形内，落在圆内的概率为p，现投针有如下结果：
投针共n次，其中投在圆内m次，那么现在要由此来估算p。根据概率p，我们可以知道得到如此投针结果的概率为：

$p=C_n^mp^m(1-p)^{n-m}$

现在的问题是已知结果如何估计概率，MLE给出的思路是既然出现了这种结果，那就让这种结果出现的概率最大化，即最大似然。

现在问题变成

$\argmax_p \,C_n^mp^m(1-p)^{n-m}$

将该式对 p 进行求导，导数为0，得到等式：

$C_n^mmp^{m-1}(1-p)^{n-m}-C_n^mp^m(n-m)(1-p)^{n-m-1}=0$

容易解得： $p=m/n$

这就是简单的使用MLE方法对分布中参数的估计，下面给出比较一般解释。

现有$X=(x_1, x_2,\cdots,x_n)$从某一分布独立采样得到，分布中参数为 $\theta$，那么概率可以表示为：
$P = P(X\vert \theta)$
根据MLE最大似然的思想，我们需要优化的问题变成了：
$\argmax P(X\vert \theta) = \argmax \prod_{i=1}^{n}P(x_i\vert \theta)$

在实际问题中连乘的形式一般计算起来比较复杂，所以一般我们采用负对数求最小的形式：
$\argmin -\log P(X\vert \theta) = \argmin -\sum_{i=1}^n P(x_i\vert \theta)$

若分布类型已知，如二项分布，高斯分布，则可得到 $p(x;\theta)$ 的表达式，上式对 $\theta$ 求导使其为0即可解出 $\theta$.

2. 最大后验估计（Maximum A Posteriori, MAP）

如果说MLE是频率学派用来求解分布中参数的方法，那么MAP则是贝叶斯学派对应的解决办法。

简单讲下频率学派和贝叶斯学派的思想的区别：对于从某一个分布中随机获得的样本，在频率学派眼里，分布的参数是不变的；而在贝叶斯学派看来恰恰相反的，他们认为现实是定的，也就是获取的样本是固定的，而随机的是分布的参数。

首先来看一看贝叶斯理论中的几个基本概念。

• 先验概率：根据以往经验和分析得到的概率。即在本次实验结果出来之前的概率。
• 后验概率：指依据得到"结果"信息所计算出的最有可能是那种事件发生。即在已知结果的前提下的概率。

下面来展示贝叶斯理论中最核心的公式——贝叶斯公式：
$P(B\vert A)=P(A\vert B)\frac{P(B)}{P(A)}$

实际上该公式可由全概率公式推导。

在该公式中，$P(B)$即先验概率，$P(B|A)$为后验概率。

举个简单的例子方便大家理解：
假设我们根据一些经验得知男生喜欢玩游戏的概率为60%，是否喜爱用$\theta$表示，喜爱则$\theta=1$，否则为0，此处我们得到了 $P(\theta=1)=0.6$。为先验概率。
再用另一个变量$x$表示男生是否有对象，若有对象则为1，否则为0。
现在知道喜欢玩游戏的男生中30%有对象，那么即 $P(x=1\vert \theta=1)=0.3, P(x=0\vert theta=1)=0.7$；不喜欢玩游戏的男生有对象的概率是50%

• 那么首先可以得出任一男生有对象的概率是 $P(x=1)=0.6\times 0.3+0.4\times 0.5=0.38$
• 现在问一个男生有对象，那么他喜欢玩游戏的概率有多大，即
  $P(\theta=1\vert x=1)=P(x=1 \vert \theta=1) \frac{P(\theta=1)}{P(x=1)}=0.3\times \frac{0.6}{0.38}$

此处要求的概率即为后验概率。

现有样本 $X=(x_1, x_2,\cdots, x_n)$ 从某一分布中独立采样得到。我们设该分布的参数为 $\theta$。那么现在先验概率为 $P(theta)$，后验概率为$P(\theta \vert X)$。 根据第一部分对MLE的分析，我们知道MLE最大化的是 $P(X\vert theta)$，而MAP最大化的是后验概率$P(\theta \vert X)$：

$\begin{aligned} \hat\theta_{map} &= \argmax_\theta P(\theta|X)\\ &= \argmax_\theta P(X\vert \theta) \frac{P(\theta)}{P(X)}\\ &=\argmin_\theta -\log P(X\vert \theta)-\log P(\theta) +\log P(X)\\ &=\argmin_\theta -\log P(X\vert \theta)-\log P(\theta) \end{aligned}$

最后两个等式去掉 P(X) 项是因为X是确定的，该项为常数不影响结果。

注意观察到，该式最后结果和MLE的优化目标只差一个 $-\log P(\theta)$

回到我们前面说的贝叶斯学派和频率学派的区别，我们说了贝叶斯学派将样本视为固定的，而参数是随机的，也就是说 $\theta$ 是一个随机的变量，它服从一定的分布。那么这个 $P(\theta)$

3. 举例区分MLE和MAP

现有样本 $X=(x_1, x_2,\cdots,x_n)$ 是从二项分布$B(m,\theta)$中独立采样得到。现在要估计参数 $\theta$。

• MLE方法
  $\begin{aligned} \hat\theta_{MLE} &=\argmin_{\theta} -\sum_{i=0}^n\log C_m^{x_i} \theta^{x_i}(1-\theta)^{m-x_i}\\ &=\argmin_{\theta} -\sum_{i=0}^n\log C_m^{x_i}+ {x_i}\log(\theta)+(m-x_i)\log(1-\theta) \end{aligned}$
  求导易得出 $\theta$.
• MAP方法
  需要对$\theta$给出先验分布，令 $\theta\sim N(0, \delta^2)$，那么有
  $P(\theta) = C\exp(-\frac{\theta^2}{2\delta^2})$
  $\begin{aligned} \hat\theta_{MLE} &=\argmin_{\theta} \sum_{i=0}^n-\log C_m^{x_i} \theta^{x_i}(1-\theta)^{m-x_i}-\log C\exp(-\frac{\theta^2}{2\delta^2})\\ &=\argmin_{\theta} \sum_{i=0}^n -\log C_m^{x_i}- {x_i}\log(\theta)-(m-x_i)\log(1-\theta)+\frac{\theta^2}{2\delta^2} \end{aligned}$
  相比较MLE而言，在优化目标中多了一个关于参数$\theta$的二次项。熟悉机器学习的朋友都知道，这就是对MLE优化目标后面添加L2正则的结果。

4. 期望最大算法（Expectation Maximum, EM）

上面简单的描述了MLE和MAP算法过后，我们来介绍一下重头戏，机器学习领域极其出名的算法——EM算法。

首先让我们来搞清楚这个算法解决的是一个什么样的问题。EM算法用于含有隐变量的概率参数模型的最大似然估计或最大后验估计。举个例子来说明这个解释：

现有 A，B两种硬币，其各自抛掷后正面朝上的概率分别为 $\theta_A,\theta_B$。现在每轮选择一枚硬币进行抛掷10次，统计正面朝上的次数。如此反复进行1000轮。要求估计$\theta_A,\theta_B$。

• 如果给定每轮抛掷的是那种硬币，那么可以对两种硬币分别估计，统计证明朝上的次数后直接使用MLE或者MAP即可；

但是现在每一轮选择的是哪枚硬币是未知的，可以视为一个变量z，z取A或B。此时z即为一个隐藏的变量。这时候我们面临一个问题：如果我们不能知道每轮选择的是A还是B，就不能估计$\theta_A,\theta_B$；而在不知道$\theta_A,\theta_B$的情况下，我们又无法估计每轮选择的硬币是A或者是B的概率。这就成为了一个先有鸡还是先有蛋的问题。

EM算法是怎么解决这一问题的呢？
EM算法核心思想分为两步：Expectation-step和Maximum-step。其中Expectation-step就是通过观测数据估计参数，然后用估计的参数值来计算似然函数的期望；在Maximum-step则是根据寻找似然函数最大化时参数的值。

我们先来给出该方法的理论推导：
需要优化的为：
$\begin{aligned} L(\theta) &= \sum_{i=1}^n \log p(x_i;\theta)\\ &= \sum_{i=1}^n \log\sum_{z}^m p(x_i, z;\theta)\\ &= \sum_{i=1}^n \log\sum_{z}^m Q_i(z)\frac{p(x_i, z;\theta)}{Q_i(z)}\\ &\ge \sum_{i=1}^n \sum_{z}^m Q_i(z) \log\frac{p(x_i, z;\theta)}{Q_i(z)} \end{aligned}$
该式子中，$z$ 为引入的隐变量，$Q_i(z)$ 表示 $x_i$

Jensen不等式：
对任意凸函数 $f(\cdot)$ 有 $f(\mathbb{E}[X])\le \mathbb E[f(X)]$
此处，$f(\cdot)$为 $log(\cdot)$，将 $\frac{p(x_i, z;\theta)}{Q_i(z)}$

记 $J(z, Q)= \sum_{i=1}^n \sum_{z}^m Q_i(z) \log\frac{p(x_i, z;\theta)}{Q_i(z)}$，我们有：
$L(\theta) \ge J(z,Q)$

下面就是EM算法根据此不等式进行操作：

1. 固定 $\theta$，调整 $Q$ 使得 $J(z,Q)$ 增大到与 $L(\theta)$
2. 固定 $Q(z)$，调整 $\theta$ 使得 $J(z,Q)$
3. 再固定 $\theta$，重复前两步，最终收敛至局部最优解。

几个问题：

1. 第一步如何调整使得 $J(z,Q)=L(\theta)$ ？
   根据不等式等号成立的条件，当且仅当满足Jensen不等式的条件时成立，即 $X=\mathbb E[X]$，也就是变量X为常数。
   那么有：
   $\begin{aligned} &\frac{p(x_i, z;\theta)}{Q_i(z)} = c\\ &\Rightarrow p(x_i,z;\theta) = c \,Q_i(z)\\ &\Rightarrow \sum_zp(x_i,z;\theta) = c\\ &\Rightarrow \frac{p(x_i, z;\theta)}{ \sum_zp(x_i,z;\theta)} = Q_i(z) （此处是将第一个式子中c代换成上一步和式）\\ &\Rightarrow Q_i(z) = \frac{p(x_i, z;\theta)}{ p(x_i;\theta)}=p(z\vert x_i;\theta) \end{aligned}$
   也即当 $Q_i(z)$ 取上式所示后验概率的时候，$J(z,Q)=L(\theta)$
2. 为什么算法会收敛到最优解（局部最优，初值选取合理可能为全局最优）？
   因为在更新 $\theta$ 时，一般采用MLE或MAP诸如此类方法，都会极大化 $J(z,Q)$。如此迭代下去一定回收敛，具体的理论性证明本文不作讨论。

EM算法思路很像坐标上升法：

对于含有多个参数时，每次选择一个维度进行优化，经过多次迭代达到收敛。这里面$\theta$ 和 $Q$