【机器学习】作业6-EM算法

EM算法

标签（空格分隔）： 机器学习

---

感谢彭先生，感谢彭先生，感谢彭先生

EM算法和朴素贝叶斯

一般机器学习算法都有一个前提，样本的所有属性都被观测到，即样本是完整的。但是在现实环境中，会有很多不完整数据。

未观测变量学名“隐变量”，令$X$为以观测变量集，$Z$表示隐变量，$\Theta$ 表示模型参数，如果要对$\Theta$做极大似然估计，即求MAX：

$$LL(\Theta | X, Z) = lnP(X, Z | \Theta)$$

$Z$为隐变量，上式无法直接求解.此时我们可以通过计算$Z$的期望，来最大化已观测数据的对数的“边际似然”

$$LL(\Theta | X) = lnP(X | \Theta) = ln\sum_Z P(X, Z | \Theta).$$

EM算法是常用的估计参数的隐变量的利器，它是一种迭代式的方法，其基本思想是：若参数 $\Theta$ 已知，则可根据训练数据推断出最优隐变量$Z$的值（E步）；反之。若$Z$的值已知，则可以方便的对参数$\Theta$做极大似然估计（M步）.

于是，以$\Theta ^ 0$ 为起点，上面的式子可以迭代以下步骤直至收敛：

• 基于$\Theta ^ t$ 推断隐变量 $Z$ 的期望， 记为 $Z ^ t$
• 基于已观测变量 $X$ 和 $Z^t$ 对参数$\Theta$做极大似然估计，记为 $\Theta ^ {t+1}$

这就是EM算法的原型

进一步，如果我们取的不是$Z$的期望，而是基于 $\Theta ^ t$ 计算隐变量Z的概率分布 $P(Z | X, \Theta ^ t)$， 则EM算法的两个步骤是：

• E步（Expectation）： 以当前参数$\Theta ^ t$ 推断隐变量分布 $P(Z | X, \Theta ^ t)$, 并计算对数似然 $LL(\Theta | X, Z)$ 关于$Z$ 的期望

$$Q(\Theta | \Theta ^ t) = E_{Z|X, \Theta ^ t} LL(\Theta | X, Z)$$

• M步（Maximization）：寻找参数最大化期望似然，即：

$$\Theta ^{t + 1}, arg_\Theta max Q(\Theta | \Theta ^ t)$$

简要来说，EM算法使用两个步骤：E步：利用当前的估计的参数值来计算对数似然的期望值；M步：寻找能使E步产生的似然期望最大的参数值。反复迭代出解

非梯度优化方法

实现和output

本次实验使用UCI的Iris数据集，数据维度为4，设前面3个维度数据正常，第4个维度存在数据缺失（50%）。
首先对数据进行预处理，然后构造低配版的朴素贝叶斯分类器。
在对最后一维数据进行处理时，仅使用其中一半的数据，然后使用EM算法估算其均值和方差。

# -*- coding: utf-8 -*-
import numpy as np


def EM(data, valid_cnt, total_cnt, eps = 1e-4):
    # :param data: 输入的一维数组
    # :param valid_cnt: 有效样本数
    # :param total_cnt: 样本总数
    # :param eps: 收敛所需精度
    # :return: avg: 隐变量的均值，theta: 隐变量的方差
    valid_data = data[0:valid_cnt]
    avg = np.sum(valid_data) / total_cnt
    theta = np.sum(np.square(valid_data)) / total_cnt - avg
    while True:
        s1 = np.sum(valid_data) + avg * (total_cnt - valid_cnt)
        s2 = np.sum(np.square(valid_data)) + (avg * avg + theta) * (total_cnt - valid_cnt)
        new_avg = s1 / total_cnt
        new_theta = s2 / total_cnt - new_avg * new_avg
        if new_avg - avg <= eps and new_theta - theta <= eps:
            break
        else:
            avg, theta = new_avg, new_theta
    return avg, theta


def elderly_man(dtype1, dtype2, latent_idx):
    # build NAIVE bayesian
    avg, var = [], []
    for idx in range(latent_idx):
        # 对隐变量之前的数据，正常计算其均值和方差
        # dim_type1 和 dim_type2 表示多维数据中的一维
        dim_type1, dim_type2 = dtype1[:, idx], dtype2[:, idx]
        avg.append([np.average(dim_type1), np.average(dim_type2)])
        var.append([np.var(dim_type1), np.var(dim_type2)])
    # 假设维度 3 的数据为隐变量，只有一半的数据是可观测的
    # 使用EM算法估计其均值和方差
    em_avg_type1, em_var_type1 = EM(data_type1[:40, latent_idx], 20, 40)
    em_avg_type2, em_var_type2 = EM(data_type2[:40, latent_idx], 20, 40)
    # 将估计得到的均值和方差加入到数组中，并返回
    avg.append([em_avg_type1, em_avg_type2])
    var.append([em_var_type1, em_var_type2])
    return avg, var


def gauss(x, avg, var):
    t = 1.0 / np.sqrt(2 * np.pi * var)
    return t * np.exp(-np.square(x - avg) / (2.0 * var))


if __name__ == '__main__':
    data_str = open('data/iris/iris.data').readlines()
    # print(data_str)
    data_type1 = np.ndarray([50, 4], np.float32)
    data_type2 = np.ndarray([50, 4], np.float32)
    for idx in range(50):
        data_type1[idx] = data_str[idx].strip('\n').split(',')[0:4]
    for idx in range(50, 100):
        data_type2[idx - 50] = data_str[idx].strip('\n').split(',')[0:4]
    a, v = elderly_man(data_type1[:40], data_type2[:40], 3)
    # 构造测试数据集，correct_times 表示测试结果准确的数据条数
    data_test = np.concatenate((data_type1[40:], data_type2[40:]))
    correct_times = 0
    for data_idx in range(len(data_test)):
        data = data_test[data_idx]
        # 数据集两类数据相同，因此先验概率均为0.5
        val_type1, val_type2 = 0.5, 0.5
        for idx in range(4):
            # 朴素贝叶斯计算
            val_type1 *= gauss(data[idx], a[idx][0], v[idx][0])
            val_type2 *= gauss(data[idx], a[idx][1], v[idx][1])
        # 前10条数据为类型1，后10条数据为类型2
        if val_type1 > val_type2 and data_idx < 10:
            correct_times += 1
        elif val_type1 < val_type2 and data_idx >= 10:
            correct_times += 1
        print("Num: %2d, Type1: %f, Type2: %f"
              % (data_idx + 1, val_type1, val_type2))
    print("Accuracy: %.1f%%" % (correct_times * 5))

output

C:\ProgramData\Anaconda3\python.exe D:/ECNU2017/J机器学习/MachineLearning/EM_Algorithm.py
Num:  1, Type1: 3.068372, Type2: 0.000000
Num:  2, Type1: 0.006664, Type2: 0.000000
Num:  3, Type1: 0.614423, Type2: 0.000000
Num:  4, Type1: 0.001331, Type2: 0.000000
Num:  5, Type1: 0.026249, Type2: 0.000000
Num:  6, Type1: 1.754030, Type2: 0.000000
Num:  7, Type1: 2.557216, Type2: 0.000000
Num:  8, Type1: 2.103980, Type2: 0.000000
Num:  9, Type1: 3.413123, Type2: 0.000000
Num: 10, Type1: 5.123086, Type2: 0.000000
Num: 11, Type1: 0.000000, Type2: 0.363503
Num: 12, Type1: 0.000000, Type2: 0.513148
Num: 13, Type1: 0.000000, Type2: 0.429917
Num: 14, Type1: 0.000000, Type2: 0.000804
Num: 15, Type1: 0.000000, Type2: 0.582627
Num: 16, Type1: 0.000000, Type2: 0.450959
Num: 17, Type1: 0.000000, Type2: 0.642538
Num: 18, Type1: 0.000000, Type2: 0.743852
Num: 19, Type1: 0.000000, Type2: 0.000821
Num: 20, Type1: 0.000000, Type2: 0.630034
Accuracy: 100.0%

Process finished with exit code 0

100%过分了

2017.11.29