贝叶斯分类器（二）贝叶斯网和EM算法

vlambda
2022-04-19

贝叶斯分类器（二）贝叶斯网和EM算法

一、贝叶斯网

贝叶斯网（Bayesian network）亦称“信念网”（belief network），它借助于有向无环图（Directed Acyclic Graph, DAG）来刻画属性之间的依赖关系，并使用条件概率表（Conditional Probability Table, CPT）来描述属性的联合概率分布。

一个贝叶斯网 $B$ 由结构 $G$ 和参数 $\Theta$ 两部分构成，即 $B=\langle G,\Theta\rangle$ 。网络结构 $G$ 是一个有向无环图，其每个结点对应于一个属性，若两个属性有直接依赖属性，则它们由一条边连接起来；参数 $\Theta$ 定量描述这种依赖关系，假设属性 $x_i$ 在 $G$ 中的父结点为 $\pi_i$ ，则 $\Theta$ 包含了每个属性的条件概率表 $\theta_{x_i|\pi_i}=P_B(x_i|\pi_i)$ 。

贝叶斯分类器（二）贝叶斯网和EM算法

贝叶斯网有效地表达了属性间的条件独立性。给定父结点集，贝叶斯网假设每个属性与它的非后裔属性独立，于是 $B=\langle G,\Theta\rangle$ ，将属性 $x_1,x_2,...,x_d$ 的联合概率分布定义为：

P_B(x_1,x_2,...,x_d)=\prod^d_{i=1}P_B(x_i|\pi_i)=\prod^d_{i=1}\theta_{x_i}|\pi_i

以上图为例，联合概率分布定义为：

P(x_1,x_2,...,x_5)=P(x_1)P(x_2)P(x_3|x_1)P(x_4|x_1,x_2)P(x_5|x_2)

$P(x_1,x_2)=\sum_{x_4}P(x_1,x_2,x_4)=\sum_{x_4}P(x_4|x_1,x_2)P(x_1)P(x_2)=P(x_1,x_2)$ 称为“边际独立性”， $x_1$ 独立于 $x_2$ 。

贝叶斯网中三个变量之间的典型依赖关系：

贝叶斯分类器（二）贝叶斯网和EM算法

在“同父”结构中，给定父结点 $x_1$ 的取值，则 $x_3$ 和 $x_4$ 条件独立。在“顺序”结构中，给定 $x$ 的值，则 $y$ 与 $z$ 条件独立。 $V$ 型结构亦称“冲撞”结构，给定子结点 $x_4$ 的取值， $x_1$ 和 $x_2$ 必不独立；奇妙的是，若 $x_4$ 的取值完全未知，则 $V$ 型结构下 $x_1$ 和 $x_2$ 却是相互独立的。

事实上，一个变量取值的确定与否，能对另两个变量间的独立性发生影响，这个现象并非 $V$ 型结构所特有。如在同父结构中，条件独立性 $x_3\perp x_4|x_1$ 成立，但若 $x_1$ 的取值未知，则 $x_3$ 和 $x_4$ 就不独立；在顺序结构中， $y\perp z|x$ ，但 $y$ 与 $z$ 不独立。

“有向分离”（D-separation），分析有向图中变量之间的条件独立性。先把有向图转变为一个无向图：

找出有向图中的所有 $V$ 型结构，在 $V$ 型结构的两个父结点之间加上一条无向边；
将所有有向边改为无向边。

由此产生的无向图称为“道德图”（moral graph），令父结点相连的过程称为“道德化”（moralization）。

贝叶斯分类器（二）贝叶斯网和EM算法

贝叶斯网学习的首要任务就是根据训练数据集找出结构最“恰当”的贝叶斯网。“评分搜索”是求解这一问题的常用方法。具体来说，先定义一个评分函数，评估贝叶斯网域训练数据的契合程度，然后基于评分函数寻找结构最优的贝叶斯网。

每个贝叶斯网描述了一个训练数据上的概率分布，自有一套编码机制能使那些经常出现的样本有更短的编码。“最小长度描述”（Minimal Description Length， MDL）准则：选择综合编码长度（包括描述网络和编码数据）最短的贝叶斯网。

给定训练集 $D=\{x_1,x_2,...,x_m\}$ ，贝叶斯网 $B=\langle G,\Theta\rangle$ 在 $D$ 上的评分函数可写为：

s(B|D)=f(\theta)|B|-LL(B|D)

其中 $|B|$ 是贝叶斯网的参数个数； $f(\theta)$ 表示描述每个参数 $\theta$ 所需的字节数；而 $LL(B|D)=\sum^m_{i=1}logP_B(x_i)$ 是贝叶斯网 $B$ 的对数似然。

显然， $s(B|D)=f(\theta)|B|-LL(B|D)$ 中的第一项是计算编码贝叶斯网 $B$ 所需的字节数，第二项是计算 $B$ 所对应的概率分布 $P_B$ 需多少字节来描述 $D$ 。于是，学习任务转化为一个优化任务，即寻找一个贝叶斯网 $B$ 使评分函数 $s(B|D)$ 最小。

若 $f(\theta)=1$ ，即每个参数用1字节描述，则得到AIC(Akaike Information Criterion)评分函数：

$AIC(B|D)=|B|-LL(B|D)$
若 $f(\theta)=\frac{1}{2}logm$ ，则得到BIC(Bayesian Information Criterion)评分函数：

$BIC(B|D)=\frac{1}{2}logm|B|-LL(B|D)$
若 $f(\theta)=0$ ，即不计算对网络进行编码的长度，则评分函数退化为负对数似然，相应的，学习任务退化为极大似然估计。

若贝叶斯网 $B=\langle G,\Theta\rangle$ 的网络结构 $G$ 固定，则评分函数 $s(B|D)$ 的第一项为常数。此时，最小化 $s(B|D)$ 等价于对参数 $\Theta$ 的极大似然估计。参数 $\theta_{x_i|\pi_i}$ 能直接在训练数据 $G$ 上通过经验估计获得，即：

\theta_{x_i|\pi_i}=\hat P_D(x_i|\pi_i)

其中 $\hat P_D(·)$ 是 $D$ 上的经验分布。最小化评分函数 $s(B|D)$ ，只需要对网络结构进行搜索，而候选结构的最优参数可直接在训练集上计算得到。

通过已知变量观测值来推测待查询变量的过程称为“推断”（inference），已知变量观测值称为“证据”（evidence）。

现实应用中，贝叶斯网的近似推断常使用吉布斯采样（Gibbs sampling）来完成，这是一种随机采样方法。

令 $Q=\{Q_1,Q_2,...,Q_n\}$ 表示待查询变量， $E=\{E_1,E_2,..,E_k\}$ 为证据变量，已知其取值为 $e=\{e_1,e_2,...e_k\}$ 。目标是计算后验概率 $P(Q=q|E=e)$ ，其中 $q=\{q_1,q_2,...,q_n\}$ 是待查询的一组取值。如待查询变量为 $Q=\{好瓜,甜度\}$ ，证据变量 $E=\{色泽,敲声,根蒂\}$ 且已知其取值为 $e=\{青绿,浊响,蜷缩\}$ ，查询的目标值是 $q=\{是,高\}$ ，即这是好瓜且甜度高的概率有多大。

下图为吉布斯采样算法：

贝叶斯分类器（二）贝叶斯网和EM算法

实质上，吉布斯采用是在贝叶斯网所有变量的联合状态空间与证据 $E=e$ 一致的子空间进行“随机漫步”（random walk）。每一步仅依赖于前一步的状态，这是一个“马尔科夫链”（Markov chain）。在一定条件下，无论从什么初始状态开始，马尔科夫链第 $t$ 步的状态分布在 $t\rightarrow \infty$ 时必收敛于一个平稳分布；对于吉布斯采样来说，这个分布恰好是 $P(Q|E=e)$ 。因此，在 $T$ 很大时，吉布斯采样相当于根据 $P(Q|E=e)$ 采样，从而保证了式 $P(Q=q|E=e)\simeq \frac{n_q}{T}$ 收敛于 $P(Q=q|E=e)$ 。

需要注意的是，由于马尔科夫链通常需要很长时间才能趋于平稳分布，因此，吉布斯采样算法的收敛速度较慢。此外，若贝叶斯网中存在极端概率“0”或“1”，则不能保证马尔科夫链存在平稳分布，此时吉布斯采样会给出错误的估计结果。

二、EM算法

概率模型有时含有观测变量（observable variable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测数据，那么给定数据，可以直接用极大似然估计或贝叶斯估计法估计模型的参数。但当模型含有隐变量时，就不能简单地使用这些估计方法。

EM算法是一种迭代方法，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。

EM算法的每次迭代由两步组成：E步，求期望（expectation）；M步，求极大（maximization），故称为期望极大算法（expectation maximization algorithm），简称EM算法。

例（三硬币模型）：假设有3枚硬币，分别记作A,B,C。这些硬币正面出现的概率分别是 $\pi,p和q$ 。进行如下实验：先掷硬币 A，根据其结果选出硬币B或硬币C，正面选硬币B，反面选硬币C；然后掷选出的硬币，掷硬币的结果，出现正面记作1，反面记作0；独立的重复 $n$ 次实验（n=10），观测结果如下：1,1,0,1,0,0,1,0,1,1。假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问，如何估计三个硬币正面出现的概率，即三硬币模型的参数。

贝叶斯分类器（二）贝叶斯网和EM算法

EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

一般地，用 $Y$ 表示观测随机变量的数据， $Z$ 表示隐随机变量的数据。 $Y$ 和 $Z$ 连在一起称为完全数据（complete-data），观测数据 $Y$ 又称为不完全数据（incomplete-data）。

EM算法通过迭代求 $L(\theta)=logP(Y|\theta)$ 的极大似然估计。每次迭代包含两步：E步，求期望；M步，求极大化。

import numpy as np
def expectation_maximization(obs_value, theta, tol, max_iter): """ EM算法 obs_value: 观测值 theta: 模型参数 tol: 精度要求 max_iter: 最大迭代次数 """ m, n = obs_value.shape for it in range(max_iter): theta_old = np.copy(theta) # 用于精度判断，退出迭代 for i in range(m): mu_array = np.zeros(n) # 存储每个观测值yi的mu值 for j in range(n): # 计算模型参数pi^i，p^i，q^i下观测数据y_i来自掷硬币B的概率 yi = obs_value[i, j] p_bcoin = theta[0] * theta[1] ** yi * (1 - theta[1]) ** (1 - yi) p_ccoin = (1 - theta[0]) * theta[2] ** yi * (1 - theta[2]) ** (1 - yi) mu_array[j] = p_bcoin / (p_bcoin + p_ccoin) # 更新参数值 theta[0] = np.mean(mu_array) theta[1] = mu_array.dot(obs_value[i]) / mu_array.sum() theta[2] = (1 - mu_array).dot(obs_value[i]) / (1 - mu_array).sum() # 满足精度，迭代退出 print("iter: ", it + 1, " model args: ", theta) if abs(theta - theta_old).max() <= tol: break def random_sample_observations(m, n): """ 随机产生m轮，每轮n次实验 """ random_obs = np.zeros((m, n)) for i in range(m): for j in range(n): random_obs[i, j] = np.random.randint(0, 2) return random_obs

observations = np.array([[1, 1, 0, 1, 0, 0, 1, 0, 1, 1]])theta = np.array([0.5, 0.5, 0.5])expectation_maximization(observations, theta, 1e-15, 10)

EM算法的描述：

输入：观测变量数据 $Y$ ，隐变量数据 $Z$ ，联合分布 $P(Y,Z|\theta)$ ，条件分布 $P(Z|Y,\theta)$ ；

输出：模型参数 $\theta$ 。

（1）选择模型的初值 $\theta^0$ ，开始迭代；

（2）E步：记 $\theta^{(i)}$ 为第 $i$ 次迭代参数 $\theta$ 的估计值，在第 $i+1$ 步迭代的E步，计算

Q(\theta,\theta^{(i)})=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}]=\sum_ZlogP(Y,Z|\theta)P(Z|Y,\theta^{(i)})

这里 $P(Z|Y,\theta^{(i)})$ 是给定观测数据 $Y$ 和当前的参数估计 $\theta^{(i)}$ 下隐变量数据 $Z$ 的条件概率分布。

（3）M步：求使 $Q(\theta,\theta^{(i)})$ 极大化的 $\theta$ ，确定第 $i+1$ 次迭代的参数的估计值 $\theta^{(i+1)}$

\theta^{(i+1)}=\underset{\theta}{argmax}Q(\theta,\theta^{(i)})

（4）重复第（2）和第（3）步，直到收敛。

函数 $Q(\theta,\theta^{(i)})$ 是EM算法的核心，称为Q函数。

步骤（1），参数的初值可以任意选择，但需要注意EM算法对初值是敏感的。
步骤（2），E步求 $Q(\theta,\theta^{(i)})$ 。Q函数式中 $Z$ 是未观测数据， $Y$ 是观测数据。注意， $Q(\theta,\theta^{(i)})$ 的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求Q函数及其极大。
步骤（3），M步求 $Q(\theta,\theta^{(i)})$ 的极大化，得到 $\theta^{(i+1)}$ ，完成一次迭代 $\theta^{(i)}\rightarrow \theta^{(i+1)}$ 。
步骤（4）给出停止迭代的条件，一般是对较小的正数 $\epsilon_1,\epsilon_2$ ，若满足

||\theta^{(i+1)}-\theta^{(i)}||<\epsilon_1

或

||Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})||<\epsilon_2

则停止迭代。

在现实应用中往往会遇到“不完整”的训练样本，这种存在“未观测”变量的情形下，是否仍能对模型参数进行估计呢？

未观测变量的学名是“隐变量”（latent variable）。令 $X$ 表示已观测变量集， $Z$ 表示隐变量集， $\Theta$ 表示模型参数。若欲对 $\Theta$ 做极大似然估计，则应最大化对数似然： $LL(\Theta|X,Z)=lnP(X,Z|\Theta)$ 。然而由于 $Z$ 是隐变量，无法直接求解。此时，可通过对 $Z$ 计算期望，来最大化已观测数据的对数“边际似然”（marginal likelihood）：

LL(\Theta|X)=lnP(X|\Theta)=ln\sum_ZP(X,Z|\Theta)

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

基于 $\Theta^t$ 推断隐变量的期望，记为 $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）：寻找参数最大化期望似然，即：