机器学习 EM算法

EM算法

复习：Jensen不等式：若f是凸函数

基本Jensen不等式
$$f(\theta x+(1-\theta )y)\le \theta f(x)+(1-\theta )f(y)$$
若${\theta }_1,\dots ,{\theta }_k\ge 0,{\theta }_1+\cdots +{\theta }_k=1$
则$f\left({\theta }_1{x}_1+\cdots +{\theta }_k{x}_k\right)\le {\theta }_{1}f\left(x_1\right)+\cdots +{\theta }_{k}f\left(x_k\right)$

若$p(x)\ge 0$ on $S\subseteq \mathrm{dom}f,{\int }_{S}p(x)dx=1$
则$f\left({\int }_{S}p(x)xdx\right)\le {\int }_{S}f(x)p(x)dx$

$f(\mathbf{E}x)\le \mathbf{E}f(x)$

引子：K-means算法

K-means算法，也被称为k-平均或k-均值，是一种广泛使用的聚类算法，或者成为其他聚类算法的基础。

K-means过程

思考

经典的K-means聚奏方法，能够非常方便的将未标记的样本分成若干簇；但无法给出某个样本属于该练的后验概率。
其他方法可否处理未标记样本呢？

最大似然估计

找出与样本的分布最接近的概率分布模型。
简单的例子
■10次抛硬币的结果是：正正反正正正反反正正
假设p是每次抛硬币结果为正的概率。则：
得到这样的实验结果的概率是：
$\begin{array}{rl}P& =pp(1-p)ppp(1-p)(1-p)pp\\ & =p^7(1-p{)}^3\end{array}$
■最优解是：p=0.7

二项分布的最大似然估计

投硬币试验中，进行N次独立试验，n次朝上，N-n次朝下。
假定朝上的概率为p，使用对数似然函数作为目标函数：$\begin{array}{rl}f(n|p)& =\log \left(p^n(1-p{)}^{N-n}\right)\stackrel{\Delta }{⟶}h(p)\\ \frac{\partial h(p)}{\partial p}& =\frac{n}{p}-\frac{N-n}{1-p}\stackrel{\Delta }{⟶}0\Rightarrow p=\frac{n}{N}\end{array}$
进一步考察
若给定一组样本${\mathrm{X}}_1,{\mathrm{X}}_2\dots {\mathrm{X}}_{\mathrm{n}}$，已知它们来自于高斯分布$\mathrm{N}(\mu ,\sigma )$，试估计参数$\mu ,\sigma$。
按照MLE的过程分析
高斯分布的概率密度函数：
$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$
将Xi的样本值xi带入，得到：
$L(x)={\prod }_{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(x_i-\mu \right)}^2}{2{\sigma }^2}}$
化简对数似然函数
$l(x)=\log {\prod }_i\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(x_i-\mu \right)}^2}{2{\sigma }^2}}$
$={\sum }_i\log \frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(x_i-\mu \right)}^2}{2{\sigma }^2}}$
$=\left({\sum }_i\log \frac{1}{\sqrt{2\pi }\sigma }\right)+\left({\sum }_i-\frac{{\left(x_i-\mu \right)}^2}{2{\sigma }^2}\right)$
$=-\frac{n}{2}\log \left(2\pi {\sigma }^2\right)-\frac{1}{2{\sigma }^2}{\sum }_i{\left(x_i-\mu \right)}^2$
参数估计的结论
目标函数$l(x)=-\frac{n}{2}\log \left(2\pi {\sigma }^2\right)-\frac{1}{2{\sigma }^2}{\sum }_i{\left(x_i-\mu \right)}^2$
将目标函数对参数$\mu ,\sigma$分别求偏导，很容易得到$\mu ,\sigma$的式子：$\{\begin{array}{l}\mu =\frac{1}{n}{\sum }_i{x}_i\\ {\sigma }^2=\frac{1}{n}{\sum }_i{\left(x_i-\mu \right)}^2\end{array}$
符合直观想象
$\{\begin{array}{l}\mu =\frac{1}{n}{\sum }_i{x}_i\\ {\sigma }^2=\frac{1}{n}{\sum }_i{\left(x_i-\mu \right)}^2\end{array}$
上述结论和矩估计的结果是一致的，并且意义非常直观：样本的均值即高斯分布的均值，样本的伪方差即高斯分布的方差。
■该结论将作为下面分析的基础。
问题：随机变量无法直接（完全）观察到

• 随机挑选10000位志愿者，测量他们的身高：若样本中存在男性和女性，身高分别服从N(${\mu }_1,{\sigma }_1$)和N(${\mu }_2,{\sigma }_2$)的分布，试估计${\mu }_1,{\sigma }_1$，${\mu }_2,{\sigma }_2$。

• 给定一幅图像，将图像的前景背景分开

• 无监督分类：聚类/EM

从直观理解猜测GMM的参数估计

随机变量X是有K个高斯分布混合而成，取各个高斯分布的概率为${\pi }_1{\pi }_2\dots {\pi }_{\mathrm{K}}$，第i个高斯分布的均值为${\mu }_i$，方差为${\Sigma }_i$。若观测到随机变量X的一条到样本${\mathrm{X}}_1,{\mathrm{X}}_2,\dots ,{\mathrm{X}}_{\mathrm{n}}$，试估计参数元$\pi$，$\mu ,{\Sigma }_{}$。

建立目标函数

对数似然函数
$l_{\pi ,\mu ,\Sigma }(x)={\sum }_{i=1}^N\log \left({\sum }_{k=1}^K{\pi }_{k}N\left(x_i|{\mu }_k,{\Sigma }_k\right)\right)$
由于在对数函数里面又有加和，无法直接用求导解方程的办法直接求得最大值。为了解决这个问题，我们分成两步。
第一步：估算数据来自哪个组份
估计数据由每个组份生成的概率：对于每个样本x，
它由第k个组份生成的概率为
$\gamma (i,k)=\frac{{\pi }_{k}N\left(x_i|{\mu }_k,{\Sigma }_k\right)}{{\sum }_{j=1}^K{\pi }_{j}N\left(x_i|{\mu }_j,{\Sigma }_j\right)}$

上式中$\mu ,{\Sigma }_{}$也是待估计的值，因此采样选代法：
在计算$\gamma (\mathbf{i},\mathbf{k})$)时假定$\mu 和{\Sigma }_{}$已知；
■需要先验给定$\mu 和{\Sigma }_{}$。
■$\gamma (\mathbf{i},\mathbf{k})$)亦可看成组份k在生成数据${\mathrm{x}}_{\mathrm{i}}$时所做的贡款。
第二步：估计每个组份的参数
对于所有的样本点，对于组份k而言，可看做生成了 { γ ( i , k ) x i ∣ i = 1 , 2 , ⋯ N } \left\{\gamma(i, k) x_{i} | i=1,2, \cdots N\right\} {γ(i,k)xi​∣i=1,2,⋯N}这些点。组份k是一个标准的高斯分布，利用上面的结论：
$u=\frac{1}{n}{\sum }_i{x}_i$
${\tau }^2=\frac{1}{n}{\sum }_i{\left(x_i-\mu \right)}^2$$\{\begin{array}{l}{N}_k={\sum }_{i=1}^N\gamma (i,k)\\ {\mu }_k=\frac{1}{N_k}{\sum }_{i=1}^N\gamma (i,k)x_i\\ {\Sigma }_k=\frac{1}{N_k}{\sum }_{i=1}^N\gamma (i,k)\left(x_i-{\mu }_k\right){\left(x_i-{\mu }_k\right)}^T\\ {\pi }_k=\frac{N_k}{N}=\frac{1}{N}{\sum }_{i=1}^N\gamma (i,k)\end{array}$

EM算法的提出

假定有训练集
$\left\{x^{(1)},x^{(2)},\cdots ,x^{(m)}\right\}$
包含m个独立样本，希望从中找到该组数据的模型p(X，z)的参数。
如：X是观测到的如身高，Z是性别，未观测到的
通过最大似然估计建立目标函数
取对数似然函数
$l(\theta )={\sum }_{i=1}^m\log p(x;\theta )$
$={\sum }_{i=1}^m\log {\sum }_{z}p(x,z;\theta )$
问题的提出$l(\theta )={\sum }_{i=1}^m\log {\sum }_{z}p(x,z;\theta )$
z是隐随机变量，不方便直接找到参数估计。
策略：计算$l(\theta )$下界，求该下界的最大值；
重复该过程，直到收敛到局部最大值。

Jensen不等式

寻找尽量紧的下界
为了使得等号成立$\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}=c$
进一步分析$\begin{array}{rl}{Q}_i\left(z^{(i)}\right)\propto p\left(x^{(i)},z^{(i)};\theta \right)& \sum _z{Q}_i\left(z^{(i)}\right)=1\\ Q_i\left(z^{(i)}\right)& =\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{{\sum }_{z}p\left(x^{(i)},z^{(i)};\theta \right)}\\ & =\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{p\left(x^{(i)};\theta \right)}\\ & =p\left(z^{(i)}|x^{(i)};\theta \right)\end{array}$

EM算法整体框架

坐标上升

总结

EM算法实现

# !/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
from scipy.stats import multivariate_normal
from sklearn.mixture import GaussianMixture
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
import matplotlib.pyplot as plt
from sklearn.metrics.pairwise import pairwise_distances_argmin


mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False


if __name__ == '__main__':
    style = 'myself'

    np.random.seed(0)
    mu1_fact = (0, 0, 0)
    cov1_fact = np.diag((1, 2, 3))
    data1 = np.random.multivariate_normal(mu1_fact, cov1_fact, 400)
    mu2_fact = (2, 2, 1)
    cov2_fact = np.array(((1, 1, 3), (1, 2, 1), (0, 0, 1)))
    data2 = np.random.multivariate_normal(mu2_fact, cov2_fact, 100)
    data = np.vstack((data1, data2))
    y = np.array([True] * 400 + [False] * 100)

    if style == 'sklearn':
        g = GaussianMixture(n_components=2, covariance_type='full', tol=1e-6, max_iter=1000)
        g.fit(data)
        print('类别概率:\t', g.weights_[0])
        print('均值:\n', g.means_, '\n')
        print('方差:\n', g.covariances_, '\n')
        mu1, mu2 = g.means_
        sigma1, sigma2 = g.covariances_
    else:
        num_iter = 100
        n, d = data.shape
        # 随机指定
        # mu1 = np.random.standard_normal(d)
        # print mu1
        # mu2 = np.random.standard_normal(d)
        # print mu2
        mu1 = data.min(axis=0)
        mu2 = data.max(axis=0)
        sigma1 = np.identity(d)
        sigma2 = np.identity(d)
        pi = 0.5
        # EM
        for i in range(num_iter):
            # E Step
            norm1 = multivariate_normal(mu1, sigma1)
            norm2 = multivariate_normal(mu2, sigma2)
            tau1 = pi * norm1.pdf(data)
            tau2 = (1 - pi) * norm2.pdf(data)
            gamma = tau1 / (tau1 + tau2)

            # M Step
            mu1 = np.dot(gamma, data) / np.sum(gamma)
            mu2 = np.dot((1 - gamma), data) / np.sum((1 - gamma))
            sigma1 = np.dot(gamma * (data - mu1).T, data - mu1) / np.sum(gamma)
            sigma2 = np.dot((1 - gamma) * (data - mu2).T, data - mu2) / np.sum(1 - gamma)
            pi = np.sum(gamma) / n
            print(i, ":\t", mu1, mu2)
        print('类别概率:\t', pi)
        print('均值:\t', mu1, mu2)
        print('方差:\n', sigma1, '\n\n', sigma2, '\n')

    # 预测分类
    norm1 = multivariate_normal(mu1, sigma1)
    norm2 = multivariate_normal(mu2, sigma2)
    tau1 = norm1.pdf(data)
    tau2 = norm2.pdf(data)

    fig = plt.figure(figsize=(13, 7), facecolor='w')
    ax = fig.add_subplot(121, projection='3d')
    ax.scatter(data[:, 0], data[:, 1], data[:, 2], c='b', s=30, marker='o', depthshade=True)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('原始数据', fontsize=18)
    ax = fig.add_subplot(122, projection='3d')
    order = pairwise_distances_argmin([mu1_fact, mu2_fact], [mu1, mu2], metric='euclidean')
    print(order)
    if order[0] == 0:
        c1 = tau1 > tau2
    else:
        c1 = tau1 < tau2
    c2 = ~c1
    acc = np.mean(y == c1)
    print('准确率：%.2f%%' % (100*acc))
    ax.scatter(data[c1, 0], data[c1, 1], data[c1, 2], c='r', s=30, marker='o', depthshade=True)
    ax.scatter(data[c2, 0], data[c2, 1], data[c2, 2], c='g', s=30, marker='^', depthshade=True)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('EM算法分类', fontsize=18)
    plt.suptitle('EM算法的实现', fontsize=21)
    plt.subplots_adjust(top=0.90)
    plt.tight_layout()
    plt.show()