试答系列：“西瓜书”-周志华《机器学习》习题试答
系列目录
[第01章：绪论]
[第02章：模型评估与选择]
[第03章：线性模型]
[第04章：决策树]
[第05章：神经网络]
[第06章：支持向量机]
第07章：贝叶斯分类器
第08章：集成学习
第09章：聚类
第10章：降维与度量学习
第11章：特征选择与稀疏学习
第12章：计算学习理论(暂缺)
第13章：半监督学习
第14章：概率图模型
(后续章节更新中...)
目录
- 7.1 试使用极大似然法估算西瓜数据集3.0中前三个属性的类条件概率。
- 7.2 试证明：条件独立性假设不成立时，朴素贝叶斯分类器仍有可能产生最优贝叶斯分类器。
- 7.3 试编程实现拉普拉斯修正的朴素贝叶斯分类器，并以西瓜数据集3.0为训练集，对p151“测1”样本进行判别。
- 7.4 实践中使用式(7.15)决定分类类别时，若数据的维度非常高，则概率连乘$\Pi_{i=1}^dP(x_i|c)$的结果通常会非常接近于0从而导致下溢。试述防止下溢的可能方案。
- 7.5 试证明：二分类任务中两类数据满足高斯分布且方差相同时，线性判别分析产生贝叶斯最优分类器。
- 7.6 试编程实现AODE分类器，并以西瓜数据集3.0为训练集，对p.151的“测1”样本进行判别。
- 7.7 给定d个二值属性的二分类任务，假设对于任何先验概率项的估算至少需30个样例，则在朴素贝叶斯分类器式(7.15)中估算先验概率项$P(c)$需30×2=60个样例。试估计在AODE式(7.23)中估算先验概率$P(c,x_i)$所需的样例数(分别考虑最好和最坏情形)。
- 7.8 考虑图7.3，试证明：在同父结构中，若x1的取值未知，则x3⊥x4不成立；在顺序结构中，y⊥z|x，但y⊥z不成立。
- 7.9 以西瓜数据集2.0为训练集，试基于BIC准则构建一个贝叶斯网。
- 7.10 以西瓜数据集2.0中属性“脐部”为隐变量，试基于EM算法构建一个贝叶斯网。
- 附：编程代码
7.1 试使用极大似然法估算西瓜数据集3.0中前三个属性的类条件概率。
答：只需对西瓜数据集3.0中具体相同类别和属性取值的样本进行计数即可得到前三个属性的类条件概率：
P青绿|是=3/8，P乌黑|是=4/8，P浅白|是=3/8，P青绿|否=3/9，P乌黑|否=2/9，P浅白|否=4/9；
P蜷缩|是=5/8，P稍蜷|是=3/8，P硬挺|是=0/8，P蜷缩|否=3/9，P稍蜷|否=4/9，P硬挺|否=2/9；
P浊响|是=6/8，P沉闷|是=2/8，P清脆|是=0/8，P浊响|否=4/9，P沉闷|否=3/9，P清脆|否=2/9；
7.2 试证明：条件独立性假设不成立时，朴素贝叶斯分类器仍有可能产生最优贝叶斯分类器。
答：可以参考“阅读材料”部分的解释。(待完善)
1、对分类任务来说，只需各类别的条件概率排序正确、无须精准概率值即可导致正确分类结果。
2、若属性间依赖对所有类别影响相同，或依赖关系的影响能相互抵消，则属性条件独立性假设在降低计算开销的同时不会对性能产生负面影响。
7.3 试编程实现拉普拉斯修正的朴素贝叶斯分类器，并以西瓜数据集3.0为训练集，对p151“测1”样本进行判别。
答： 代码附后。
其中编写了函数c=nb(x,X,Y,laplace=True)，可以通过参数laplace选择是否进行拉普拉斯修正。将其设为False时，可以将计算结果与教材计算结果进行对比，对比发现，教材中对于连续型属性-----密度和含糖率的正态分布参数估计中，对于方差的估计采用了无偏估计，亦即\(\sigma_{c,i}^{2}=\frac{1}{|D_c|-1}\sum_{x\in D_c}(x-\mu_{c,i})^2\)。
7.4 实践中使用式(7.15)决定分类类别时，若数据的维度非常高，则概率连乘\(\Pi_{i=1}^dP(x_i|c)\)的结果通常会非常接近于0从而导致下溢。试述防止下溢的可能方案。
答： 通常采用取对数的方法将连乘变为连加：\(\Pi_{i=1}^dP(x_i|c)\rightarrow log[\Pi_{i=1}^dP(x_i|c)]=\Sigma_{i=1}^dlogP(x_i|c)\)
7.5 试证明：二分类任务中两类数据满足高斯分布且方差相同时，线性判别分析产生贝叶斯最优分类器。
答：假设数据满足高斯分布：\(P(x|c)\sim N(\mu_c,\Sigma)\)，模型中需要确定的参数有均值\(\mu_1\)和\(\mu_0\)，以及共同的方差\(\Sigma\)，\(\Sigma\)为对称正定矩阵。数据集的对数似然为：
\[\begin{aligned} LL(\mu_1,\mu_0,\Sigma)&=\sum_ilogP(x_i|c_i)\\ &=\sum_ilog\{\frac{1}{(2\pi)^{d/2}|\Sigma|^{1/2}}exp[-\frac{1}{2}(x_i-\mu_{c_i})^T\Sigma^{-1}(x_i-\mu_{c_i})]\} \end{aligned} \]
通过最大化对数似然可以求得参数估计：
\[\begin{aligned} \nabla_{\mu_c}LL=0 &\Rightarrow \mu_c=\frac{1}{|D_c|}\sum_{x_i \in D_c}x_i\\ \nabla_{\Sigma^{-1}}LL=0 &\Rightarrow \Sigma=\frac{1}{m}\sum_i(x_i-\mu_{c_i})(x_i-\mu_{c_i})^T \end{aligned} \]
上式中求取\(\Sigma^{-1}\)梯度时应用了关系\(\nabla_A|A|=|A|(A^{-1})^T\)。
那么，该贝叶斯分类器的决策函数为：
\[h_{Bayes}(x)=argmax_cP(c)P(x|c) \]
对于二分类任务，这等价于：
\[\begin{aligned} h_{Bayes}(x)&=sign[P(1)P(x|1)-P(0)P(x|0)]\\ &=sign\{exp[-\frac{1}{2}(x-\mu_1)^T\Sigma^{-1}(x-\mu_1)]-exp[-\frac{1}{2}(x-\mu_0)^T\Sigma^{-1}(x-\mu_0)]\}\\ &=sign[\mu_0^T\Sigma^{-1}\mu_0-\mu_1^T\Sigma^{-1}\mu_1+2(\mu_1-\mu_0)^T\Sigma^{-1}x] \end{aligned} \]
上式中第二行采取了同先验假设，亦即P(0)=P(1)=1/2.
在3.4节线性判别分析(LDA)中，关于\(\mu_1,\mu_0\)的定义与上面求得的\(\mu_c\)完全相同，而根据(3.33)式可知，\(S_w=m\Sigma\)。在3.4节中求得最优投影直线方向为\(w=S_w^{-1}(\mu_0-\mu_1)\)， LDA对于新数据的分类是根据投影点距离两个投影中心的距离远近决定的，可以将其表达为:
\[\begin{aligned} h_{LDA}(x) &=sign[(w^Tx-w^T\mu_0)^2-(w^Tx-w^T\mu_1)^2]\\ &=sign\{[2w^Tx-w^T(\mu_1+\mu_0)][w^T(\mu_1-\mu_0)]\} \end{aligned} \]
注意到上式第二行右边项\(w^T(\mu_1-\mu_0)=-(\mu_1-\mu_0)^TS_w^{-1}(\mu_1-\mu_0)\)，而\(S_w\)为对称正定矩阵，因此该项恒为负，因此，上式可以进一步化简为：
\[\begin{aligned}h_{LDA}(x) &=sign[w^T(\mu_1+\mu_0)-2w^Tx]\\ &=sign[(\mu_0-\mu_1)^T\Sigma^{-1}(\mu_1+\mu_0-2x)]\\ &=sign[\mu_0^T\Sigma^{-1}\mu_0-\mu_1^T\Sigma^{-1}\mu_1+2(\mu_1-\mu_0)^T\Sigma^{-1}x] \end{aligned}\]
对比可知，决策函数\(h_{Bayes}(x)\)和\(h_{LDA}(x)\)完全相同，因此可以说LDA产生了最优Bayes分类。
7.6 试编程实现AODE分类器，并以西瓜数据集3.0为训练集，对p.151的“测1”样本进行判别。
答： 编程代码附后。
在实现AODE分类器过程中，需要应用(7.23)~(7.25)式计算各个概率值。由于西瓜数据集3.0中含有连续属性，需要对(7.24)和(7.25)计算式进行相应的拓展。
-
对于(7.24)式，若\(x_i\)为连续属性，则按下式计算：
\[P(c,x_i)=P(c)P(x_i|c)=\frac{|D_c|}{|D|}\cdot N(x_i|\mu_{c,i},\sigma_{c,i}^2) \]其中\(N\)为正态分布的概率密度函数，参数为\(N(x|\mu,\sigma^2)\)。
-
对于(7.25)式，存在多种情况：
(a). 若\(x_i\)为离散属性，\(x_j\)为连续属性，此时：\[P(x_j|c,x_i)=N(x_j|\mu_{c,x_i;j},\sigma_{c,x_i;j}^2) \](b). 若\(x_i\)为连续属性，\(x_j\)为离散属性，此时：
\[\begin{aligned}P(x_j|c,x_i)&=\frac{P(x_i|c,x_j)P(c,x_j)}{P(c,x_i)}\\ &=N(x_i|\mu_{c,x_j;i},\sigma_{c,x_j;i}^2)\cdot\frac{P(c,x_j)}{P(c,x_i)} \end{aligned}\](c.). 若\(x_i\)为连续属性，\(x_j\)为连续属性，此时：
\[\begin{aligned} P(x_j|c,x_i)&=\frac{P(x_i,x_j|c)P(c)}{P(c,x_i)}\\ &=N[(x_i,x_j)|\mu_{c,(i,j)},\Sigma_{c,(i,j)}]\cdot\frac{P(c)}{P(c,x_i)} \end{aligned}\]需要注意的是，对于上面(a),(b)两种情况下，求取正态分布参数时，可能因为子集\(D_{c,x_i}\)中样本数太少而无法估计。在编程中，若样本数为零，则假定为正态分布，亦即\(\mu=0,\sigma=1\)；若仅一个样本，则将平均值设为该唯一取值，方差设为1.
7.7 给定d个二值属性的二分类任务，假设对于任何先验概率项的估算至少需30个样例，则在朴素贝叶斯分类器式(7.15)中估算先验概率项\(P(c)\)需30×2=60个样例。试估计在AODE式(7.23)中估算先验概率\(P(c,x_i)\)所需的样例数(分别考虑最好和最坏情形)。
答：需要计算的先验概率\(P(c,x_i)\)的数目共有2×2×d个，因为c有2个取值，i有d个取值，\(x_i\)有2个取值。
首先考虑属性\(x_1\)，对于参数\(P(c=0,x_1=0)，P(c=0,x_1=1)，P(c=1,x_1=0)，P(c=1,x_1=1)\)，分别需要30个样例来估算，总计需要120个样例。
在最好情形下，对于所有其他属性\(x_i，i=2,3,…,d\)，假设刚好能在c=0的60个样例中找到30个\(x_i\)=0和30个\(x_i\)=1的样例，同样地，对于c=1的情况也如此，则无需更多样例，此时对于所有属性总计需要120个样例；
在最坏情形下，比如对于属性\(x_2\)，假设60个c=0的样例中\(x_2\)取值全部为0，则需要另外补充30个\(c=0,x_2=1\)的样例，同样的，对于c=1的情况，也需要补充30个样例。考虑所有其他属性，总共要补充60×(d-1)个样例，这样总计需要120+60(d-1)=60(d+1)个样例。
因此，对于最好和最坏情形，分别需要样例数120和60(d+1)个样例。
7.8 考虑图7.3，试证明：在同父结构中，若x1的取值未知，则x3⊥x4不成立；在顺序结构中，y⊥z|x，但y⊥z不成立。
证明：
- 在同父结构中。假设x3⊥x4成立，则有：
\(P(x_3,x_4)=P(x_3)P(x_4)=\sum_{x_1}P(x_1)P(x_3|x_1)P(x_4)\)
另一方面有：
\(P(x_3,x_4)=\sum_{x_1}P(x_1,x_3,x_4)=\sum_{x_1}P(x_1)P(x_3|x_1)P(x_4|x_1)\)
两式相减有：
\(\sum_{x_1}P(x_1)P(x_3|x_1)[P(x_4)-P(x_4|x_1)]=0\)
然而上式未必成立，因此x3⊥x4成立未必成立。 - 在顺序结构中。首先有联合概率为：
\(P(x,y,z)=P(z)P(x|z)P(y|x)\)
可以证明：
\(P(y,z|x)=\frac{P(x,y,z)}{P(x)}=\frac{P(z)P(x|z)P(y|x)}{P(x)}=\frac{P(x,z)P(y|x)}{P(x)}=P(z|x)P(y|x)\)
上式即表明 z⊥y | x.
至于z⊥y是否成立，先假设其成立，则有：
\(P(y,z)=P(y)P(z)=P(z)\sum_xP(y|x)P(x)\)
又因为：
\(P(y,z)=\sum_xP(x,y,z)=\sum_xP(z)P(x|z)P(y|x)=P(z)\sum_xP(x|z)P(y|x)\)
两式相减有：
\(\sum_xP(y|x)[P(x)-P(x|z)]=0\)
然而上式未必成立，因此z⊥y未必成立。
7.9 以西瓜数据集2.0为训练集，试基于BIC准则构建一个贝叶斯网。
答：
1、贝叶斯网络结构在计算机中的表示。
首先，我们约定一下如何在计算机中表达一个贝叶斯网络：设有\(n\)个节点，\(x_1,x_2,…,x_n\)，利用\(B\in R^{n×n}\)矩阵的右上角元素来表达连接边，右上角元素\(B_{ij}\)可能取值为0,1,-1，取值0意味着\(x_i\)与\(x_j\)之间没有连接边，等于1意味着\(x_i \rightarrow x_j\)的连接边，等于-1意味着\(x_i \leftarrow x_j\)的连接边。比如，将下图中的贝叶斯网络用矩阵表达为：
2、BIC评分函数。
在贝叶斯网络中，对于某个节点\(x_i\)，设它的可能取值数目为\(r_i\)，它的父节点集为\(\pi_i\),父节点取值的可能组合数目为\(q_i\)。在数据集D中，父节点集\(\pi_i\)取值组合为第\(j\)种组合的样本数目为\(m_{ij}\)，与此同时，节点\(x_i\)取第\(k\)个值的样本数为\(m_{ijk}\)。
那么，对于某个网络结构B，它的BIC评分函数计算公式可以表示[1]：
\[\begin{aligned} BIC(B|D)&=\frac{logm}{2} |B|-LL(B|D) \\ &=\frac {logm}{2} \sum_{i=1}^{n}q_i(r_i-1)-\sum_{i=1}^{n}\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}m_{ijk}log\frac{m_{ijk}}{m_{ij}}\\ \end{aligned}\tag{1}\]
下面试着对上式进行详细说明。
我们用\(\theta_{ijk}=P(x_i=k|\pi_i=j)\)来表示网络参数：对于第i个结点，父节点取值为第j种组合，并且该结点取第k个离散值的概率。\(\theta_{ijk}\)对应于第i格概率表，第j行，第k列的取值，如下图所示：
BIC评分公式右边第一项为结构风险项，|B|表示网络中的参数数目，对于结点xi处的概率表包含\(q_ir_i\)个参数，但是考虑到概率表存在归一化的约束条件:\(\Sigma_k\theta_{ijk}=1\)，因此，独立的参数个数有\(q_i(r_i-1)\)个。
第二项为经验风险项：负对数似然概率。下面试着对其进行详细推导：
\[\begin{aligned} LL(B|D)&=\sum_{s=1}^mlogP(X^{(s)};\theta)\\ &=\sum_{s=1}^m\sum_{i=1}^nlogP(x_i^{(s)}|\pi_i^{(s)})\\ &=\sum_{s=1}^m\sum_{i=1}^nlog\theta_{i,\pi_i^{(s)},x_i^{(s)}} \end{aligned}\tag{2}\]
上式是对很多个概率值对数\(logθ_{ijk}\)的累加求和，其计算结果必然是这样的形式：
\[LL(B|D)=\sum_i\sum_j\sum_km_{ijk}log\theta_{ijk}\tag{3} \]
其中\(m_{ijk}\)为正整数，它的意义恰好代表：对于某个节点\(x_i\)，父节点集\(\pi_i\)取第j种组合，节点\(x_i\)取第k个离散值的样本数。
考虑到网络参数\(θ_{ijk}\)应该满足约束条件：\(\Sigma_k\theta_{ijk}=1\)，利用带等式约束的拉格朗日方法对似然函数求取极大值来求解网络参数：
\[\begin{aligned} \max_\theta L&=LL(\theta|D)+\sum_i\sum_j \lambda_{ij}(1-\sum_k\theta_{ijk})\\ &=\sum_i\sum_j\sum_km_{ijk}log\theta_{ijk}+\sum_i\sum_j\lambda_{ij}(1-\sum_k\theta_{ijk}) \end{aligned}\tag{4}\]
其中\(\lambda_{ij}\)为拉格朗日参数。对上式求导数，并令其等于零：
\[\begin{aligned} &\frac{\partial L}{\partial\theta_{ijk}}=\frac{m_{ijk}}{\theta_{ijk}}-\lambda_{ij}=0\\ &\Rightarrow \theta_{ijk}=\frac{m_{ijk}}{\lambda_{ij}} \end{aligned}\tag{5}\]
考虑约束条件\(\sum_k\theta_{ijk}=1\)，可以得到\(\lambda_{ij}=\sum_k m_{ijk}=m_{ij}\)，这里\(m_{ij}\)恰好代表着对于某个节点\(x_i\)，父节点集\(π_i\)取第j种组合的样本数。因此，网络参数为：
\[\theta_{ijk}=\frac{m_{ijk}}{m_{ij}}\tag{6} \]
这就是为什么在固定网络结构下，网络参数简单地等于对应取值样本的“计数”。
将网络参数(6)带回到前面的(3)式即得到：
\[LL(B|D)=\sum_{i=1}\sum_{j=1}\sum_{k=1}m_{ijk}log\frac{m_{ijk}}{m_{ij}}\tag{7} \]
3、搜索算法。
我们的目标是搜寻BIC评分最小化的贝叶斯网络结构，这里采用下山法进行搜索，每一步随机选取一个边进行调整，调整包括三种操作[1:1]：加边、减边、转向。需要注意的是，加边和转向后可能会引入环，要予以避免。
如果新网络的BIC评分结果更低，则进行更新，否则维持原网络。同时允许一个较小的概率调整为BIC分值更大的网络。算法如下：
初始化网络结构B，
设置参数eta和tao，
for loop in range(step):
随机选取需要调整连接边的两个节点xi和xj，
对Bij值进行修改，
if 新网络有环：
continue
if BIC值减小 or rand()<eta*exp(-loop/tao):
更新网络，输出当前BIC值
else:
continue
输出最终得到的网络结构和参数。
详细的编程代码附后。
说明一点：观察前面的BIC评分函数计算式可以发现，总的BIC分值为各个节点处BIC分值之和。因此，在比较调整前后BIC分值变化时，由于网络仅在节点xi和xj处发生变化，为了减小计算量，可以仅计算这两个节点处的BIC分值进行比较。
4、计算结果。
下面是在西瓜数据集2.0上某一次运行的结果：
初始BIC评分：277.572(结构203.991,经验73.580)
最终BIC评分：119.360(结构 43.915,经验75.445)
多次运行，在不同初始网络的情况下，结果均表现出一个规律：“色泽”是一个独立的属性，是否是“好瓜”主要由“纹理”和“触感”两个属性决定。
7.10 以西瓜数据集2.0中属性“脐部”为隐变量，试基于EM算法构建一个贝叶斯网。
答：
1、EM算法总结
先来总结一下EM算法，大概是这么回事：我们要建立一个概率模型，它具有参数θ，能够告诉我们任意一种变量取值下的概率分布，但是这个模型所涉及的变量中既包含可观测变量X，又包含不可观测变量Z(或者说隐变量)。在利用关于X的数据集D来估计参数θ时只能采用“边际似然”最大化：
\[max LL(\theta|X)=logP(X|\theta)=log\sum_zP(X,Z|\theta)\tag{1} \]
按理说，求上式的极大便可以得到所需估计的参数θ，但是上式往往求解较为困难，于是便有了EM算法。
按照周志华《机器学习》上的介绍，没有过多的推导，完全是一种直观的认识，很容易理解。EM算法是重复以下两个步骤：假如已经知道了参数θ，便可以根据参数θ来估计隐变量Z的分布P(Z|X;θ)或者最佳取值，这便是E步骤(求期望)；假设知道了隐变量Z的分布或者最佳取值，便可以按照完整数据的极大似然的方法来求解θ，这便是M步(求极大)。
但是在其他很多书籍中，介绍的EM算法是这样推导得到的：
\[\begin{aligned} LL(\theta|X)&=\sum_{i=1}^mlog\sum_zP(X^{(i)},z;\theta)\\ &=\sum_{i=1}^mlog\sum_zQ^{(i)}(z)\frac{P(X^{(i)},z;\theta)}{Q^{(i)}(z)}\\ &\geq\sum_{i=1}^m\sum_zQ^{(i)}(z)log(\frac{P(X^{(i)},z;\theta)}{Q^{(i)}(z)}) \end{aligned}\]
上式中Q(Z)是关于Z的一个分布，EM算法过程为：
\[E-Step: Q(z)=P(z|X;\theta)\tag{2} \]
\[M-Step:\theta=argmax_\theta\sum_{i=1}^m\sum_zQ^{(i)}(z)log(\frac{P(X^{(i)},z;\theta)}{Q^{(i)}(z)})\tag{3} \]
将上式与教材中的(7.37)式对比，貌似不同，实际上两者是等价的，说明如下：
教材中的(7.36)式可以详细表达为:
\[\begin{aligned} ELL(\theta|\theta^t)&=E_{Z|X,\theta^t}LL(\theta|X,Z)\\ &=\sum_ZP(Z|X,\theta^t)LL(\theta|X,Z)\\ &=\sum_ZQ(Z)logP(X,Z;\theta) \end{aligned}\]
上式中将教材中所采用的字母Q改为ELL，以防与前文中的Q(Z)相混淆；
将教材中表示参数的Θ换成小写θ表示，以与前文保持一致。
上式是对于单个样本而言的，对于m个样本有：
\[ELL(\theta|\theta^t)=\sum_{i=1}^m\sum_ZQ^{(i)}(Z)logP(X^{(i)},Z;\theta) \]
然后(7.37)式即为：
\[\begin{aligned} \theta^{t+1}&=argmax_\theta ELL(\theta|\theta^t)\\ &=argmax_\theta \sum_{i=1}^m\sum_ZQ^{(i)}(Z)logP(X^{(i)},Z;\theta) \end{aligned}\tag{4}\]
由于在M步骤中将\(Q^{(i)}(Z)\)视为常量，所以(4)式和(3)式所得到的\(θ^{(t+1)}\)完全一样。
2. 本题解答思路。
具体在贝叶斯网络模型中情况又有所不同，除了计算网络参数，还需要搜索网络结构。回顾上一题7.9中，基于BIC准则搜索最佳贝叶斯网络，算法中主要关注于网络结构的搜索，而网络参数的确定则通过简单计数来确定，并隐含在BIC评分函数中。
上题中的算法流程为：
本题解答采用完全一样的算法流程，只不过在计算BIC分值这个步骤需要多做一些工作，此时需要(同样是在固定的网络结构下)通过EM算法来确定网络参数，然后再计算出BIC评分分值。
3. EM算法下计算网络参数和BIC评分函数。
本题中将属性“脐部”作为隐变量，单个隐变量的情况应该较为简单，下面具体分析。
3.1 E-Step： 求取P(Z|X,θ)。也就是说，根据当前的网络参数来确定“脐部”的条件概率分布，“脐部”分别取值凹陷、平坦、稍凹的概率是多少。
\[\begin{aligned} Q(z)&=P(z|X;\theta^t)\\ &=\frac{P(z,X;\theta^t)}{P(X;\theta^t)}\\ &\propto P(z,X;\theta^t)\\ &=\prod_{x_i}P(x_i|\pi_i)\\ &\propto P(z|\pi_z)\prod_{x_i\in son(z)}P(x_i|\pi_i) \end{aligned}\tag{5}\]
上式中多处引入正比例符号∝，这是因为只需要计算z在各种取值下的相对概率大小，然后再进行归一化即可，因此与z无关的部分，或者说对于所有z相同的部分可以不予考虑。
上式中倒数第二行对\(x_i\)的求积是对所有结点的求积，包括隐变量结点z。这些乘积因子中只有z结点及其子节点与z有关，其中son(z)是指结点z的子结点。因此，隐变量的条件概率仅与隐变量结点本身及其子结点处的网络参数(概率表)有关。
比如，在下图的网络结构中，只有Z结点和\(X_5\)结点处的概率表影响Q(z)的计算结果。
3.2 M-Step： 根据上一步所确定的Q(z)，通过最大化似然函数的期望来求取网络参数。
\[\begin{aligned} \max ELL(\theta|\theta^t)&=\sum_{s=1}^m\sum_z Q^{(s)}(z)logP(X^{(s)},z;\theta)\\ &=\sum_{s=1}^m\sum_z Q^{(s)}(z)\sum_i logP(x_i^{(s)}|\pi_i^{(s)})\\ &=\sum_i[\sum_{s=1}^m\sum_z Q^{(s)}(z) logP(x_i^{(s)}|\pi_i^{(s)})]\\ &=\sum_i\sum_{s=1}^m\sum_z Q^{(s)}(z) log\theta_{i,\pi_i^{(s)},x_i^{(s)}} \end{aligned}\]
与前边一样，上式中从第2行开始对i的求和是指对所有结点的求和，这其中包括了隐变量结点Z；而对z的求和是指对z不同取值的求和。
其中网络参数\(θ_{i,j,k}\)与上一题7.9中的含义一样：对于第i个结点，父节点取值为第j种组合，并且该结点取第k个离散值的概率。
若令：
\[\overline m_{ijk}=\sum_{s=1}^m\sum_z Q_z^{(s)}\mathbf {II}(\pi_i^{(s)}=j,x_i^{(s)}=k) \]
其中\(\mathbf {II}(·)\)是指示函数，与教材中含义相同，表示在\(\{Q_s^{(s)}\}\)分布下，结点xi取第k个离散值，对应的父节点取第j种组合的样本数的期望。当xi是隐结点或者其子结点时，πi或者xi的取值与z有关，分布Qz会影响到相应结果。对于其他结点，Qz分布与之无关，\(\overline m_{ijk}\)变为\(m_{ijk}\)。
那么，前面的式子可以继续表示为：
\[\begin{aligned} \max ELL(\theta|\theta^t)&=\sum_i\sum_{s=1}^m\sum_z Q^{(s)}(z) log\theta_{i,\pi_i^{(s)},x_i^{(s)}}\\ &=\sum_i\sum_j\sum_k\overline m_{ijk}log\theta_{ijk} \end{aligned}\]
与7.9题中求取网络参数的过程完全一样，我们可以得到网络参数为：
\[\begin{aligned} \theta_{ijk}^{t+1}&=argmax_\theta ELL(\theta|\theta^t)\\ &=\frac{\overline m_{ijk}}{\sum_k \overline m_{ijk}}\\ &=\frac{\overline m_{ijk}}{\overline m_{ij}} \end{aligned}\tag{6}\]
注意到，z的不同取值分布只会影响到隐结点z本身以及它的子节点处\(\overline m_{ijk}\)的数值，因此，仅隐结点及其子节点处的网络参数会进行更新，其余结点处的网络参数始终保持不变。
3.3 合并E-step和M-step过程
EM算法需要不断的重复E-step和M-step两个过程，在这个过程中，Qz和θ不断地被更新：
\[Q_z^0\rightarrow\theta^1\rightarrow Q_z^1\rightarrow\theta^2\rightarrow Q_z^2\rightarrow\theta^3\rightarrow\dots \]
Qz相较于θ参数方便于表示和存储，可以将Qz存储于m×|z|的矩阵中，|z|为z的状态数。我们可以将这两个过程合并，直接考察Qz随着迭代过程的变化规律：
Qz0→Qz1→Qz2→...
\[Q_z^0\rightarrow Q_z^1\rightarrow Q_z^2\rightarrow\dots \]
将(5)式和(6)式进行合并，比如对于第r个样本，z取值z1的分布更新为：
\[\begin{aligned} Q^{[t+1](r)}(z=z_1)&\propto P(X^{(r)},z=z_1;\theta)\\ &=\prod_i P(x_i^{(r)}|\pi_i^{(r)})\\ &=\prod_i\theta_{i,\pi_i^{(r)},x_i^{(r)}}\\ &=\prod_i\frac{\overline m_{i,\pi_i^{(r)},x_i^{(r)}}}{\overline m_{i,\pi_i^{(r)}}}\\ &\propto\overline m_{z,\pi_z^{(r)},z_1}\prod_{i\in son(z)}\frac{\overline m_{i,\pi_i^{(r)},x_i^{(r)}}}{\overline m_{i,\pi_i^{(r)}}} \end{aligned}\]
上式第二行开始对i的求和包括对Z结点，其中隐含了z=z1(暂未发现恰当的表示方法)。最后一行只取其中与z有关的项：z结点处的mijk与z有关，而mij与z无关，z子结点的mijk和mij都与z有关。比如具体来说，对于下图中的贝叶斯网络，Z结点的子结点只有x5，那么有：
\[\begin{aligned} Q^{[t+1](r)}(z=z_1)\propto&\overline m_{z,\pi_z^{(r)},z_1}\cdot\frac{\overline m_{5,\pi_5^{(r)},x_5^{(r)}}}{\overline m_{5,\pi_5^{(r)}}}\\ =&\sum_{s=1}^m Q^{[t](s)}(z=z_1)\mathbf {II}(x_1^{(s)}=x_1^{(r)},x_2^{(s)}=x_2^{(r)})\\ &\cdot\frac{\sum_{s=1}^mQ^{[t](s)}(z=z_1)\mathbf {II}(x_4^{(s)}=x_4^{(r)},x_5^{(s)}=x_5^{(r)})}{\sum_{s=1}^mQ^{[t](s)}(z=z_1)\mathbf {II}(x_4^{(s)}=x_4^{(r)})} \end{aligned}\]
在第t+1步更新Qz时的算法总结如下：
输入：上一步得到的Qz
数据集D
Z结点的子结点索引号son(z)
各个结点的父结点索引号πi
过程：
NewQz=ones(size(Qz)) #新的Qz初始化为全1矩阵
#考察Z结点
for j in |πz|： #遍历Z结点父结点的所有取值组合
Index=第j种取值的样本索引号
NewQz[index]*=sum(Qz[index])
#考察Z的子结点
for i in son(z): #遍历Z的子结点
for j in |πi+xi-z|: #遍历(父结点+xi结点-z结点)的所有取值组合
Index=第j种取值的样本索引号
NewQz[index]*=sum(Qz[index])
for j in |πi-z|: #遍历(父结点-z结点)的所有取值组合
Index=第j种取值的样本索引号
NewQz[index]/=sum(Qz[index])
NewQz/=sum(NewQz,axis=1)
Qz=NewQz
3.4 BIC评分函数
上题7.9中的BIC评分函数包含结构项和经验项，其中结构项保持不变，将经验项由“似然函数”变为“边际似然函数”：
\[\begin{aligned} BIC(B|D)&=\frac{logm}{2}|B|-LL(B|D)\\ &=\frac{logm}{2}\sum_{i=1}^n q_i(r_i-1)-\sum_{s=1}^m log\sum_z P(X^{(s)},z|\theta) \end{aligned}\]
其中结构项部分与上题完全一样，其中ri表示结点xi的可能取值数目，qi为父节点集取值可能组合数目。
下面着重分析经验项部分：
\[\begin{aligned} LL(B|D)&=\sum_{s=1}^m log\sum_z P(X^{(s)},z|\theta)\\ &=\sum_{s=1}^m log\sum_z Q_z^{(s)} \frac{P(X^{(s)},z|\theta)}{Q_z^{(s)}}\\ &\geq\sum_{s=1}^m\sum_z Q_z^{(s)} log\frac{P(X^{(s)},z|\theta)}{Q_z^{(s)}} \end{aligned}\]
当EM算法收敛时，对于最终的网络参数θ和隐变量分布\(\{Q_z^{(s)}\}\)，上式最后一行的不等式取等号。上式进一步化简为：
\[\begin{aligned} LL(B|D)&=\sum_{s=1}^m\sum_z Q_z^{(s)} log P(X^{(s)},z|\theta)-\sum_{s=1}^m\sum_z Q_z^{(s)} log Q_z^{(s)}\\ &=\sum_{s=1}^m\sum_z Q_z^{(s)} \sum_i log P(X_i^{(s)}|\pi_i^{(s)})-\sum_{s=1}^m\sum_z Q_z^{(s)} log Q_z^{(s)}\\ &=\sum_i\sum_j\sum_k\overline m_{ijk}log\theta_{ijk}-\sum_{s=1}^m\sum_z Q_z^{(s)} log Q_z^{(s)}\\ &=\sum_i\sum_j\sum_k\overline m_{ijk}log\frac{\overline m_{ijk}}{\overline m_{ij}}-\sum_{s=1}^m\sum_z Q_z^{(s)} log Q_z^{(s)}\\ &=\sum_i\sum_j\sum_k\overline m_{ijk}log\overline m_{ijk}-\sum_i\sum_j\overline m_{ij}log\overline m_{ij}-\sum_{s=1}^m\sum_z Q_z^{(s)} log Q_z^{(s)}\\ \end{aligned}\]
观察上式，前两项仍然是对各个结点求和的形式，第三项仅与隐结点z的取值分布函数Q(z)有关。
4.编程计算
代码附后，在上一题代码的基础上增加EM算法以及新的BIC评分函数计算方法。
说明几点：
① 设置不同的初始值，EM算法可能收敛到不同的结果，可以尝试不同的初始值，最后取边际似然最大的结果。
② 在更新Qz的过程中，将θ的更新隐含在其中，\(\theta_{ijk}=\frac{\overline m_{ijk}}{\overline m_{ij}}\)，当计算到的mij=0时，会出现0/0的错误。这可以通过laplace修正来解决，但在这里，我采取在mij上统一加一个较小的数，比如1E-100，来避免错误。
③ 在调整网络结构后，如果隐结点及其子结点处概率表不发生变化，可以无需运行EM算法来求取Qz。BIC评分只涉及到有限几个结点处的变化，可以利用这些特点来减少计算量。
5.结果讨论
从前面的分析可以看出，z的不同取值具有对称性，也就是我们让z1=“凹陷”还是“平坦”无关紧要，起影响作用的只有z的状态数(可能取值数目)，本题中状态数为3。计算到的Qz也可以在不同状态之间进行轮换，比如Qz→Qz[:,[2,0,1]]，不影响结果。
BIC评分包含两部分，结构项和经验项，与其他算法类似，可以通过系数来调整两者的相对比重：BIC_score=struct*alpha+emp。下面是设置不同alpha值的情况下搜索到的贝叶斯网络结构：
结果表明：alpha越大，结构项比重越大，搜索到的结构越简单。alpha越小，搜索到的结构越复杂；alpha=1时，结果表明“脐部”为一个独立属性，alpha取其他值时，“脐部”均作为多个属性的根属性。
附：编程代码
习题7.3(Python)
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 9 09:53:04 2020
@author: MS
"""
import numpy as np
def nb(x,X,Y,laplace=True):
#朴素贝叶斯分类器
# x 待预测样本
# X 训练样本特征值
# Y 训练样本标记
# laplace 是否采用“拉普拉斯修正”，默认为真
C=list(set(Y)) #所有可能标记
p=[] #存储概率值
print('===============朴素贝叶斯分类器===============')
for c in C:
Xc=[X[i] for i in range(len(Y)) if Y[i]==c] #c类样本子集
pc=(len(Xc)+laplace)/(len(X)+laplace*len(C)) #类先验概率
print('P(c=%s)=%.3f;\nP(xi|c=%s)='%(c,pc,c),end='')
logp=np.log(pc) #对数联合概率P(x,c)
for i in range(len(x)):
if type(x[i])!=type(X[0][i]):
print(
'样本数据第%d个特征的数据类型与训练样本数据类型不符,无法预测!'%(i+1))
return None
if type(x[i])==str:
# 若为字符型特征，按离散取值处理
Xci=[1 for xc in Xc if xc[i]==x[i]] #c类子集中第i个特征与待预测样本一样的子集
pxc=(len(Xci)+laplace)/(len(Xc) #pxc为类条件概率P(xi|c)
+laplace*len(set([x[i] for x in X])))
print('%.3f'%pxc,end=',')
elif type(x[i])==float:
# 若为浮点型特征，按高斯分布处理
Xci=[xc[i] for xc in Xc]
u=np.mean(Xci)
sigma=np.std(Xci,ddof=1)
pxc=1/np.sqrt(2*np.pi)/\
sigma*np.exp(-(x[i]-u)**2/2/sigma**2)
print('%.3f(~N(%.3f,%.3f^2))'%(pxc,u,sigma),end=',')
else:
print('目前只能处理字符型和浮点型数据，对于其他类型有待扩展相应功能。')
return None
logp+=np.log(pxc)
p.append(logp)
print(';\nlog(P(x,c=%s))=log(%.3E)=%.3f\n'%(c,np.exp(logp),logp))
predict=C[p.index(max(p))]
print('===============预测结果：%s类==============='%predict)
return predict
#====================================
# 主程序
#====================================
# 表4.3 西瓜数据集3.0
#FeatureName=['色泽','根蒂','敲声','纹理','脐部','触感','密度','含糖率']
#LabelName={1:'好瓜',0:'坏瓜'}
X=[['青绿','蜷缩','浊响','清晰','凹陷','硬滑',0.697,0.460],
['乌黑','蜷缩','沉闷','清晰','凹陷','硬滑',0.774,0.376],
['乌黑','蜷缩','浊响','清晰','凹陷','硬滑',0.634,0.264],
['青绿','蜷缩','沉闷','清晰','凹陷','硬滑',0.608,0.318],
['浅白','蜷缩','浊响','清晰','凹陷','硬滑',0.556,0.215],
['青绿','稍蜷','浊响','清晰','稍凹','软粘',0.403,0.237],
['乌黑','稍蜷','浊响','稍糊','稍凹','软粘',0.481,0.149],
['乌黑','稍蜷','浊响','清晰','稍凹','硬滑',0.437,0.211],
['乌黑','稍蜷','沉闷','稍糊','稍凹','硬滑',0.666,0.091],
['青绿','硬挺','清脆','清晰','平坦','软粘',0.243,0.267],
['浅白','硬挺','清脆','模糊','平坦','硬滑',0.245,0.057],
['浅白','蜷缩','浊响','模糊','平坦','软粘',0.343,0.099],
['青绿','稍蜷','浊响','稍糊','凹陷','硬滑',0.639,0.161],
['浅白','稍蜷','沉闷','稍糊','凹陷','硬滑',0.657,0.198],
['乌黑','稍蜷','浊响','清晰','稍凹','软粘',0.360,0.370],
['浅白','蜷缩','浊响','模糊','平坦','硬滑',0.593,0.042],
['青绿','蜷缩','沉闷','稍糊','稍凹','硬滑',0.719,0.103]]
Y=[1]*8+[0]*9
x=['青绿','蜷缩','浊响','清晰','凹陷','硬滑',0.697,0.460] #测试例"测1"
print("测1样本：")
nb(x,X,Y,False) #此时不用拉普拉斯修正，方便与教材对比计算结果
#若用拉普拉斯修正，可以去掉最后一个参数，或者设置为True
#任意设置一个新的"测例"
x=['浅白','蜷缩','沉闷','稍糊','平坦','硬滑',0.5,0.1]
print("\n任设的一个新测例，观察结果：")
nb(x,X,Y)
习题7.6(Python)
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 9 09:53:04 2020
@author: MS
"""
import numpy as np
def AODE(x,X,Y,laplace=True,mmin=3):
#平均独依赖贝叶斯分类器
# x 待预测样本
# X 训练样本特征值
# Y 训练样本标记
# laplace 是否采用“拉普拉斯修正”，默认为真
# mmin 作为父属性最少需要的样本数
C=list(set(Y)) #所有可能标记
N=[len(set([x[i] for x in X])) for i in range(len(x))] #各个属性的可能取值数
p=[] #存储概率值
print('===============平均独依赖贝叶斯分类器(AODE)===============')
for c in C:
#--------求取类先验概率P(c)--------
Xc=[X[i] for i in range(len(Y)) if Y[i]==c] #c类样本子集
pc=(len(Xc)+laplace)/(len(X)+laplace*len(C)) #类先验概率
print('P(c=%s)=%.3f;'%(c,pc))
#--------求取父属性概率P(c,xi)--------
p_cxi=np.zeros(len(x)) #将计算结果存储在一维向量p_cxi中
for i in range(len(x)):
if type(x[i])!=type(X[0][i]):
print(
'样本数据第%d个特征的数据类型与训练样本数据类型不符,无法预测!'%(i+1))
return None
if type(x[i])==str:
# 若为字符型特征，按离散取值处理
Xci=[1 for xc in Xc if xc[i]==x[i]] #c类子集中第i个特征与待预测样本一样的子集
p_cxi[i]=(len(Xci)+laplace)/(len(X)
+laplace*len(C)*N[i])
elif type(x[i])==float:
# 若为浮点型特征，按高斯分布处理
Xci=[xc[i] for xc in Xc]
u=np.mean(Xci)
sigma=np.std(Xci,ddof=1)
p_cxi[i]=pc/np.sqrt(2*np.pi)/\
sigma*np.exp(-(x[i]-u)**2/2/sigma**2)
else:
print('目前只能处理字符型和浮点型数据，对于其他类型有待扩展相应功能。')
return None
print(''.join(['p(c=%d,xi)='%c]+['%.3f'%p_cxi[i]+
(lambda i:';' if i==len(x)-1 else ',')(i) for i in range(len(x))]))
#--------求取父属性条件依赖概率P(xj|c,xi)--------
p_cxixj=np.eye(len(x)) #将计算结果存储在二维向量p_cxixj中
for i in range(len(x)):
for j in range(len(x)):
if i==j:
continue
#------根据xi和xj是离散还是连续属性分为多种情况-----
if type(x[i])==str and type(x[j])==str:
Xci=[xc for xc in Xc if xc[i]==x[i]]
Xcij=[1 for xci in Xci if xci[j]==x[j]]
p_cxixj[i,j]=(len(Xcij)+laplace)/(len(Xci)+laplace*N[j])
if type(x[i])==str and type(x[j])==float:
Xci=[xc for xc in Xc if xc[i]==x[i]]
Xcij=[xci[j] for xci in Xci]
#若子集Dc,xi数目少于2个，则无法用于估计正态分布参数，
#则将其设为标准正态分布
if len(Xci)==0:
u=0
else:
u=np.mean(Xcij)
if len(Xci)<2:
sigma=1
else:
sigma=np.std(Xcij,ddof=1)
p_cxixj[i,j]=1/np.sqrt(2*np.pi)/\
sigma*np.exp(-(x[j]-u)**2/2/sigma**2)
if type(x[i])==float and type(x[j])==str:
Xcj=[xc for xc in Xc if xc[j]==x[j]]
Xcji=[xcj[i] for xcj in Xcj]
if len(Xcj)==0:
u=0
else:
u=np.mean(Xcji)
if len(Xcj)<2:
sigma=1
else:
sigma=np.std(Xcji,ddof=1)
p_cxixj[i,j]=1/np.sqrt(2*np.pi)/sigma*np.exp(-(x[i]-u)**2/2/sigma**2)*p_cxi[j]/p_cxi[i]
if type(x[i])==float and type(x[j])==float:
Xcij=np.array([[xc[i],xc[j]] for xc in Xc])
u=Xcij.mean(axis=0).reshape(1,-1)
sigma2=(Xcij-u).T.dot(Xcij-u)/(Xcij.shape[0]-1)
p_cxixj[i,j]=1/2*np.pi/np.sqrt(np.linalg.det(sigma2))\
*np.exp(-0.5*([x[i],x[j]]-u).
dot(np.linalg.inv(sigma2)).
dot(([x[i],x[j]]-u).T))*pc/p_cxi[i]
print(''.join([(lambda j:'p(xj|c=%d,x%d)='%(c,i+1)if j==0 else '')(j)
+'%.3f'%p_cxixj[i][j]
+(lambda i:';\n' if i==len(x)-1 else ',')(j)
for i in range(len(x)) for j in range(len(x))]))
#--------求计算总的概率∑iP(c,xi)*∏jP(xj|c,xi)--------
sump=0
for i in range(len(x)):
if len([1 for xi in X if xi[1]==x[1]])>=mmin:
sump+=p_cxi[i]*p_cxixj[i,:].prod()
print('P(c=%d,x) ∝ %.3E\n'%(c,sump))
p.append(sump)
#--------做预测--------
predict=C[p.index(max(p))]
print('===============预测结果：%s类==============='%predict)
return predict
#====================================
# 主程序
#====================================
# 表4.3 西瓜数据集3.0
#FeatureName=['色泽','根蒂','敲声','纹理','脐部','触感','密度','含糖率']
#LabelName={1:'好瓜',0:'坏瓜'}
X=[['青绿','蜷缩','浊响','清晰','凹陷','硬滑',0.697,0.460],
['乌黑','蜷缩','沉闷','清晰','凹陷','硬滑',0.774,0.376],
['乌黑','蜷缩','浊响','清晰','凹陷','硬滑',0.634,0.264],
['青绿','蜷缩','沉闷','清晰','凹陷','硬滑',0.608,0.318],
['浅白','蜷缩','浊响','清晰','凹陷','硬滑',0.556,0.215],
['青绿','稍蜷','浊响','清晰','稍凹','软粘',0.403,0.237],
['乌黑','稍蜷','浊响','稍糊','稍凹','软粘',0.481,0.149],
['乌黑','稍蜷','浊响','清晰','稍凹','硬滑',0.437,0.211],
['乌黑','稍蜷','沉闷','稍糊','稍凹','硬滑',0.666,0.091],
['青绿','硬挺','清脆','清晰','平坦','软粘',0.243,0.267],
['浅白','硬挺','清脆','模糊','平坦','硬滑',0.245,0.057],
['浅白','蜷缩','浊响','模糊','平坦','软粘',0.343,0.099],
['青绿','稍蜷','浊响','稍糊','凹陷','硬滑',0.639,0.161],
['浅白','稍蜷','沉闷','稍糊','凹陷','硬滑',0.657,0.198],
['乌黑','稍蜷','浊响','清晰','稍凹','软粘',0.360,0.370],
['浅白','蜷缩','浊响','模糊','平坦','硬滑',0.593,0.042],
['青绿','蜷缩','沉闷','稍糊','稍凹','硬滑',0.719,0.103]]
Y=[1]*8+[0]*9
x=['青绿','蜷缩','浊响','清晰','凹陷','硬滑',0.697,0.460] #测试例"测1"
print("测1样本：")
AODE(x,X,Y) #此时不用拉普拉斯修正，方便与教材对比计算结果
#若用拉普拉斯修正，可以去掉最后一个参数，或者设置为True
#任意设置一个新的"测例"
x=['浅白','蜷缩','沉闷','稍糊','平坦','硬滑',0.5,0.1]
print("\n任设的一个新测例，观察结果：")
AODE(x,X,Y)
习题7.9(Python)
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 15 17:02:12 2020
@author: lsly
"""
import numpy as np
import matplotlib.pyplot as plt
#==============首先编写几个函数，主程序见后==============
def relationship(net):
# 计算网络中的每个结点的父母结点以及父母以上的祖辈结点
# 输入：
# net:array类型，网络结构，右上角元素ij表示各个连接边
# 取值0表示无边，取值1表示Xi->Xj,取值-1表示Xi<-Xj
# 输出：
# parents：list类型，存储各个结点的父节点编号，用取值1~N代表各个节点
# grands:list类型，存储各个结点更深的依赖节点，可以看成是“祖结点”
# circle:list类型，存储环节点编号，若图中存在环，那么这个结点将是它本身的“祖结点”
N=len(net) #节点数
#-----确定父结点-----
parents=[list(np.where(net[i,:]==-1)[0]+1)+
list(np.where(net[:,i]==1)[0]+1)
for i in range(N)]
grands=[]
#-----确定“祖结点”-----
for i in range(N):
grand=[]
#---爷爷辈---
for j in parents[i]:
for k in parents[j-1]:
if k not in grand:
grand.append(k)
#---曾祖及以上辈---
loop=True
while loop:
loop=False
for j in grand:
for k in parents[j-1]:
if k not in grand:
grand.append(k)
loop=True
grands.append(grand)
#-----判断环结点-----
circle=[i+1 for i in range(N) if i+1 in grands[i]]
return parents,grands,circle
def draw(net,name=None,title=''):
# 绘制贝叶斯网络的变量关系图
# net:array类型，网络结构，右上角元素ij表示各个连接边
# 取值0表示无边，取值1表示Xi->Xj,取值-1表示Xi<-Xj
# name:指定各个节点的名称，默认为None，用x1,x2...xN作为名称
N=net.shape[0]
Level=np.ones(N,dtype=int)
#-----确定层级-----
for i in range(N):
for j in range(i+1,N):
if net[i][j]==1 and Level[j]<=Level[i]:
Level[j]+=1
if net[i][j]==-1 and Level[i]<=Level[j]:
Level[i]+=1
#-----确定横向坐标-----
position=np.zeros(N)
for i in set(Level):
num=sum(Level==i)
position[Level==i]=np.linspace(-(num-1)/2,(num-1)/2,num)
#-----画图-----
plt.figure()
plt.title(title)
#设置出图显示中文
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
#--先画出各个结点---
for i in range(N):
if name==None:
text='x%d'%(i+1)
else:
text=name[i]
plt.annotate(text,xy=[position[i],Level[i]],bbox={'boxstyle':'circle','fc':'1'},ha='center')
#--再画连接线---
for i in range(N):
for j in range(i+1,N):
if net[i][j]==1:
xy=np.array([position[j],Level[j]])
xytext=np.array([position[i],Level[i]])
if net[i][j]==-1:
xy=np.array([position[i],Level[i]])
xytext=np.array([position[j],Level[j]])
if net[i][j]!=0:
L=np.sqrt(sum((xy-xytext)**2))
xy=xy-(xy-xytext)*0.2/L
xytext=xytext+(xy-xytext)*0.2/L
if (xy[0]==xytext[0] and abs(xy[1]-xytext[1])>1)or\
(xy[1]==xytext[1] and abs(xy[0]-xytext[0])>1):
arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=0.3')
#画弧线，避免遮挡(只考虑了横向和竖向边，暂未考虑斜向边遮挡的情况)
else:
arrowprops=dict(arrowstyle='->')
plt.annotate('',xy=xy,xytext=xytext,arrowprops=arrowprops,ha='center')
plt.axis([position.min(),position.max(),Level.min(),Level.max()+0.5])
plt.axis('off')
plt.show()
def BIC_score(net,D,consider=None):
# 计算评分函数
# 输入：
# net:贝叶斯网络
# D:数据集
# 输出：
# [struct,emp]:评分函数的结构项和经验项
#-----从数据集D中提取一些信息-----
m,N=D.shape #样本数和特征数
values=[np.unique(D[:,i])for i in range(len(D[0]))] #各个离散属性的可能取值
#-----父节点-----
parents=[list(np.where(net[i,:]==-1)[0]+1)+
list(np.where(net[:,i]==1)[0]+1)
for i in range(N)]
#-----计算BIC评分-----
struct,emp=0,0 #BIC评分的结构项和经验项
if consider==None:
consider=range(N)
for i in consider:
r=len(values[i]) #Xi结点的取值数
pa=parents[i] #Xi的父节点编号
nums=[len(values[p-1]) for p in pa] #父节点取值数
q=np.prod(nums) #父节点取值组合数
struct+=q*(r-1)/2*np.log(m) #对结构项的贡献
#-----如果父节点数目为零，亦即没有父节点
if len(pa)==0:
for value_k in values[i]:
Dk=D[D[:,i]==value_k] #D中Xi取值v_k的子集
mk=len(Dk) #Dk子集大小
if mk>0:
emp-=mk*np.log(mk/m) #对经验项的贡献
continue
#-----有父节点时，通过编码方式，遍历所有父节点取值组合
for code in range(q):
#解码：比如，父节点有2×3种组合，
#将0~5解码为[0,0]~[1,2]
code0=code
decode=[]
for step in range(len(pa)-1):
wight=np.prod(nums[step+1:])
decode.append(code0//wight)
code0=code0%wight
decode.append(code0)
# 父节点取某一组合时的子集
index=range(m) #子集索引号，初始为全集D
#起初误将m写为N，该错误极不容易发现，两天后发现并更正
for j in range(len(pa)):
indexj=np.where(D[:,pa[j]-1]==values[pa[j]-1][decode[j]])[0]
index=np.intersect1d(index,indexj)
Dij=D[index,:] #子集
mij=len(Dij) #子集大小
if mij>0: #仅当子集非空时才计算该种情况
for value_k in values[i]:
Dijk=Dij[Dij[:,i]==value_k] #Dij中Xi取值v_k的子集
mijk=len(Dijk) #Dijk子集大小
if mijk>0:
emp-=mijk*np.log(mijk/mij) #对经验项的贡献
return np.array([struct,emp])
#================================================
# 主程序
#================================================
#==============西瓜数据集2.0======================
# 将X和类标记Y放一起
D=[['青绿','蜷缩','浊响','清晰','凹陷','硬滑','是'],
['乌黑','蜷缩','沉闷','清晰','凹陷','硬滑','是'],
['乌黑','蜷缩','浊响','清晰','凹陷','硬滑','是'],
['青绿','蜷缩','沉闷','清晰','凹陷','硬滑','是'],
['浅白','蜷缩','浊响','清晰','凹陷','硬滑','是'],
['青绿','稍蜷','浊响','清晰','稍凹','软粘','是'],
['乌黑','稍蜷','浊响','稍糊','稍凹','软粘','是'],
['乌黑','稍蜷','浊响','清晰','稍凹','硬滑','是'],
['乌黑','稍蜷','沉闷','稍糊','稍凹','硬滑','否'],
['青绿','硬挺','清脆','清晰','平坦','软粘','否'],
['浅白','硬挺','清脆','模糊','平坦','硬滑','否'],
['浅白','蜷缩','浊响','模糊','平坦','软粘','否'],
['青绿','稍蜷','浊响','稍糊','凹陷','硬滑','否'],
['浅白','稍蜷','沉闷','稍糊','凹陷','硬滑','否'],
['乌黑','稍蜷','浊响','清晰','稍凹','软粘','否'],
['浅白','蜷缩','浊响','模糊','平坦','硬滑','否'],
['青绿','蜷缩','沉闷','稍糊','稍凹','硬滑','否']]
D=np.array(D)
FeatureName=['色泽','根蒂','敲声','纹理','脐部','触感','好瓜']
#=================初始化贝叶斯网结构=============
#构建贝叶斯网络，右上角元素ij表示各个连接边
#取值0表示无边，取值1表示Xi->Xj,取值-1表示Xi<-Xj
m,N=D.shape
net=np.zeros((N,N))
choose=4 #选择初始化类型，可选1,2,3,4
#分别代表全独立网络、朴素贝叶斯网络、全连接网络,随机网络
if choose==1: #全独立网络:图中没有任何连接边
pass
elif choose==2: #朴素贝叶斯网络:所有其他特征的父节点都是类标记"好瓜"
net[:-1,-1]=-1
elif choose==3: #全连接网络：任意两个节点都有连接边
again=True #若存在环，则重新生成
while again:
for i in range(N-1):
net[i,i+1:]=np.random.randint(0,2,N-i-1)*2-1
again=len(relationship(net)[2])!=0
elif choose==4: #随机网络：任意两个节点之间的连接边可有可无
again=True #若存在环，则重新生成
while again:
for i in range(N-1):
net[i,i+1:]=np.random.randint(-1,2,N-i-1)
again=len(relationship(net)[2])!=0
draw(net,FeatureName,'初始网络结构')
#==============下山法搜寻BIC下降的贝叶斯网==========
score0=BIC_score(net,D)
score=[score0]
print('===========训练贝叶斯网============')
print('初始BIC评分:%.3f(结构%.3f,经验%.3f)'%(sum(score0),score0[0],score0[1]))
eta,tao=0.1,50 #允许eta的概率调整到BIC评分增大的网络
#阈值随迭代次数指数下降
for loop in range(10000):
# 随机指定需要调整的连接边的两个节点：Xi和Xj
i,j=np.random.randint(0,N,2)
while i==j:
i,j=np.random.randint(0,N,2)
i,j=(i,j) if i<j else (j,i)
# 确定需要调整的结果
value0=net[i,j]
change=np.random.randint(2)*2-1 #结果为+1或-1
value1=(value0+1+change)%3 -1 #调整后的取值
net1=net.copy()
net1[i,j]=value1
if value1!=0 and len(relationship(net1)[2])!=0:
#如果value1取值非零，说明为转向或者增边
#若引入环，则放弃该调整
continue
delta_score=BIC_score(net1,D,[i,j])-BIC_score(net,D,[i,j])
if sum(delta_score)<0 or np.random.rand()<eta*np.exp(-loop/tao):
score0=score0+delta_score
score.append(score0)
print('调整后BIC评分:%.3f(结构%.3f,经验%.3f)'
%(sum(score0),score0[0],score0[1]))
net=net1
else:
continue
draw(net,FeatureName,'最终网络结构')
score=np.array(score)
x=np.arange(len(score))
plt.title('BIC贝叶斯网络结构搜索过程')
plt.xlabel('更新次数')
plt.ylabel('分值')
plt.plot(x,score[:,0],'.r-')
plt.plot(x,score[:,1],'.b-')
plt.plot(x,score.sum(axis=1),'.k-')
plt.legend(['struct','emp','BIC(struct+emp)'])
plt.show()
习题7.10(Python)
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 15 17:02:12 2020
@author: lsly
"""
import numpy as np
import matplotlib.pyplot as plt
#==============首先编写几个函数，主程序见后==============
def relationship(net):
# 计算网络中的每个结点的父母结点以及父母以上的祖辈结点
# 输入：
# net:array类型，网络结构，右上角元素ij表示各个连接边
# 取值0表示无边，取值1表示Xi->Xj,取值-1表示Xi<-Xj
# 输出：
# parents：list类型，存储各个结点的父节点编号，用取值1~N代表各个节点
# grands:list类型，存储各个结点更深的依赖节点，可以看成是“祖结点”
# circle:list类型，存储环节点编号，若图中存在环，那么这个结点将是它本身的“祖结点”
N=len(net) #节点数
#-----确定父结点-----
parents=[list(np.where(net[i,:]==-1)[0]+1)+
list(np.where(net[:,i]==1)[0]+1)
for i in range(N)]
grands=[]
#-----确定“祖结点”-----
for i in range(N):
grand=[]
#---爷爷辈---
for j in parents[i]:
for k in parents[j-1]:
if k not in grand:
grand.append(k)
#---曾祖及以上辈---
loop=True
while loop:
loop=False
for j in grand:
for k in parents[j-1]:
if k not in grand:
grand.append(k)
loop=True
grands.append(grand)
#-----判断环结点-----
circle=[i+1 for i in range(N) if i+1 in grands[i]]
return parents,grands,circle
def draw(net,name=None,title=''):
# 绘制贝叶斯网络的变量关系图
# net:array类型，网络结构，右上角元素ij表示各个连接边
# 取值0表示无边，取值1表示Xi->Xj,取值-1表示Xi<-Xj
# name:指定各个节点的名称，默认为None，用x1,x2...xN作为名称
N=net.shape[0]
Level=np.ones(N,dtype=int)
#-----确定层级-----
for i in range(N):
for j in range(i+1,N):
if net[i][j]==1 and Level[j]<=Level[i]:
Level[j]+=1
if net[i][j]==-1 and Level[i]<=Level[j]:
Level[i]+=1
#-----确定横向坐标-----
position=np.zeros(N)
for i in set(Level):
num=sum(Level==i)
position[Level==i]=np.linspace(-(num-1)/2,(num-1)/2,num)
#-----画图-----
plt.figure()
plt.title(title)
#设置出图显示中文
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
#--先画出各个结点---
for i in range(N):
if name==None:
text='x%d'%(i+1)
else:
text=name[i]
plt.annotate(text,xy=[position[i],Level[i]],bbox={'boxstyle':'circle','fc':'1'},ha='center')
#--再画连接线---
for i in range(N):
for j in range(i+1,N):
if net[i][j]==1:
xy=np.array([position[j],Level[j]])
xytext=np.array([position[i],Level[i]])
if net[i][j]==-1:
xy=np.array([position[i],Level[i]])
xytext=np.array([position[j],Level[j]])
if net[i][j]!=0:
L=np.sqrt(sum((xy-xytext)**2))
xy=xy-(xy-xytext)*0.2/L
xytext=xytext+(xy-xytext)*0.2/L
if (xy[0]==xytext[0] and abs(xy[1]-xytext[1])>1)or\
(xy[1]==xytext[1] and abs(xy[0]-xytext[0])>1):
arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=0.3')
#画弧线，避免遮挡(只考虑了横向和竖向边，暂未考虑斜向边遮挡的情况)
else:
arrowprops=dict(arrowstyle='->')
plt.annotate('',xy=xy,xytext=xytext,arrowprops=arrowprops,ha='center')
plt.axis([position.min(),position.max(),Level.min(),Level.max()+0.5])
plt.axis('off')
plt.show()
def coder(StateNums):
# 编码器,
# 设有结点x1,x2,...xN,各个结点的状态数为s1,s2,...sN,
# 那么结点取值的组合数目为s1*s2*...sN,
# 这些组合状态可以编码表示为[0,0,...0]~[s1-1,s2-1,...sN-1]
# 输出：
# StateNums：各个结点状态数,
# 比如[2,3,2]意为x1,x2,x3分别有2,3,2种状态，
# 组合起来便有12种状态。
# 输出：
# codes:用于遍历所有状态的索引编号，
# 比如，对于[2,3],总共6种组合状态，遍历这6种组合状态的编码为：
# [0,0],[0,1],[0,2],[1,0],[1,1],[1,2]
Nodes=len(StateNums) #结点数
states=np.prod(StateNums) #组合状态数
codes=[]
for s in range(states):
s0=s
code=[]
for step in range(Nodes-1):
wight=np.prod(StateNums[step+1:])
code.append(s0//wight)
s0=s0%wight
code.append(s0)
codes.append(code)
return codes
def EM(net,D,ZStateNum,Try=1):
# EM算法计算隐变量概率分布Q(z),这里仅考虑单个隐变量的简单情况
# 输入：
# net:贝叶斯网络结构，以矩阵右上角元素表示连接关系，
# 约定将隐变量排在最后一个。
# D:可观测变量数据集
# ZStateNum:隐变量状态数(离散取值数目)结果
# Try:尝试次数，由于EM算法收敛到的结果受初始值影响较大，
# 因此，尝试不同初始值，最终选择边际似然最大的。
# 输出：
# Qz：隐变量概率分布
#=====网络性质=====
parents=[list(np.where(net[i,:]==-1)[0])+
list(np.where(net[:,i]==1)[0])
for i in range(len(net))] #计算各个结点的父节点
#=====可观测变量参数=====
m,Nx=D.shape #样本数和可观测变量数
values=[np.unique(D[:,i])for i in range(Nx)] #可观测变量的离散取值
#=====隐变量子节点=====
Zsonindex=list(np.where((net[:Nx,Nx:]==-1).any(axis=1))[0]) #隐结点子节点索引号
#=====运行多次EM，每次随机初始化Qz，最终选择边际似然最大的结果=====
for t in range(Try):
#=====隐变量分布初始化=====
Qz=np.random.rand(m,ZStateNum) #初始化隐变量概率分布
Qz=Qz/Qz.sum(axis=1).reshape(-1,1) #概率归一化
#Qz=np.c_[np.ones([m,1]),np.zeros([m,2])]
#=====迭代更新Qz=====
dif=1 #两次Qz的差别
while dif>1E-8:
NewQz=np.ones(Qz.shape)
#-----对于隐结点-----
pa=parents[-1] #隐结点的父结点
if len(parents[-1])==0: #如果隐结点没有父节点
NewQz*=Qz.sum(axis=0)
else:
ValueNums=[len(values[p]) for p in pa] #各个父结点的状态数
codes=coder(ValueNums) #用于遍历各种取值的编码
for code in codes:
#父结点取值组合
CombValue=[values[pa[p]][code[p]] for p in range(len(pa))]
index=np.where((D[:,pa]==CombValue).all(axis=1))[0]
NewQz[index]*=Qz[index].sum(axis=0) if len(index)!=0 else 1
#-----对于隐结点的子结点-----
for son in Zsonindex:
#分子部分
pa=parents[son] #父结点,
Nodes=pa+[son] #加上该结点本身,
Nodes.remove(Nx) #然后，移去隐结点作为考察结点
ValueNums=[len(values[N]) for N in Nodes] #各个结点的状态数
codes=coder(ValueNums) #用于遍历各种取值的编码
for code in codes:
CombValue=[values[Nodes[N]][code[N]] for N in range(len(Nodes))]
index=np.where((D[:,Nodes]==CombValue).all(axis=1))[0]
NewQz[index]*=Qz[index].sum(axis=0) if len(index)!=0 else 1
#分母部分
pa=parents[son]+[] #仅考察父结点
pa.remove(Nx) #移去隐结点
if len(pa)==0: #如果父结点只有隐结点一个
NewQz/=Qz.sum(axis=0)
else:
ValueNums=[len(values[p]) for p in pa] #各个父结点的状态数
codes=coder(ValueNums) #用于遍历各种取值的编码
for code in codes:
#父结点取值组合
CombValue=[values[pa[p]][code[p]] for p in range(len(pa))]
index=np.where((D[:,pa]==CombValue).all(axis=1))[0]
NewQz[index]/=Qz[index].sum(axis=0)+1E-100 if len(index)!=0 else 1
NewQz=NewQz/NewQz.sum(axis=1).reshape(-1,1) #归一化
dif=np.sum((Qz-NewQz)**2,axis=1).mean() #新旧Qz的差别
Qz=NewQz #更新Qz
if t==0:
BestQz=Qz
maxLL=LL(net,D,Qz,consider=(Zsonindex+[Nx]))
else:
newLL=LL(net,D,Qz,consider=(Zsonindex+[Nx]))
if newLL>maxLL:
maxLL=newLL
BestQz=Qz
return BestQz
def LL(net,D,Qz,consider=None):
# 含有单个隐变量的情况下，计算边际似然
# 输入：
# net:贝叶斯网络结构，以矩阵右上角元素表示连接关系，
# 约定将隐变量排在最后一个。
# D:可观测变量数据集
# Qz：隐变量概率分布
# consider:所考察的结点。根据分析，
# 边际似然中部分项可以表示为各个结点求和的形式，
# 因此可以指定求和所包含的结点
# 输出：
# LL：边际似然
#=====网络性质=====
parents=[list(np.where(net[i,:]==-1)[0])+
list(np.where(net[:,i]==1)[0])
for i in range(len(net))] #计算各个结点的父节点
#=====可观测变量参数=====
m,Nx=D.shape #样本数和可观测变量数
values=[np.unique(D[:,i])for i in range(Nx)] #可观测变量的离散取值
#=====待考察结点=====
if consider is None:
consider=range(Nx+1)
#=====计算边际似然的求和项=====
LL=0
#print(consider)
for i in consider:
#print(i)
pa=parents[i] #父结点
sign=1
for nodes in [pa+[i],pa]: #nodes为当前所考察的结点
if len(nodes)==0: #考虑当前xi没有父结点的情况
LL+=sign*m*np.log(m)
continue
zin=nodes.count(Nx)!=0 #是否含有隐结点
if zin:
nodes.remove(Nx)
if len(nodes)==0: #除了隐结点外没有其他结点
mz=Qz.sum(axis=0)
LL+=sign*sum(np.log(mz**mz))
else:
StateNums=[len(values[nd]) for nd in nodes]
for code in coder(StateNums):
CombValue=[values[nodes[N]][code[N]] for N in range(len(nodes))]
index=np.where((D[:,nodes]==CombValue).all(axis=1))[0]
if zin:
mz=Qz[index].sum(axis=0)
LL+=sign*sum(np.log(mz**mz))
else:
mz=len(index)
LL+=sign*(np.log(mz**mz))
sign*=-1
#=====计算隐变量概率分布项=====
LL-=np.sum(np.log(Qz**Qz))
return LL
def BIC(net,D,Qz,alpha=1,consider=None):
# 计算BIC评分
# 输入：
# net:贝叶斯网络结构，以矩阵右上角元素表示连接关系，
# 约定将隐变量排在最后一个。
# D:可观测变量数据集
# Qz：隐变量分布
# alpha:结构项的比重系数
# consider:所考察的结点。根据分析，
# BIC评分中前两项可以表示为各个结点求和的形式，
# 因此可以指定求和所包含的结点
# 输出：
# np.array([struct,emp])：BIC评分结构项和经验项两部分的分值
#-----从数据集D中提取一些信息-----
m,Nx=D.shape #样本数和可观测变量数特征数
values=[list(np.unique(D[:,i])) for i in range(len(D[0]))] #各个离散属性的可能取值
values.append(range(ZStateNum)) #再加上隐变量的取值数
#-----父节点-----
parents=[list(np.where(net[i,:]==-1)[0])+
list(np.where(net[:,i]==1)[0])
for i in range(len(net))]
#-----计算BIC评分-----
emp=-LL(net,D,Qz,consider) #经验项部分调用LL函数来计算
struct=0 #下面计算结构项
if consider is None:
consider=range(Nx+1)
for i in consider:
r=len(values[i]) #Xi结点的取值数
pa=parents[i] #Xi的父节点编号
nums=[len(values[p]) for p in pa] #父节点取值数
q=np.prod(nums) #父节点取值组合数
struct+=q*(r-1)/2*np.log(m) #对结构项的贡献
return np.array([struct*alpha,emp])
def BIC_change(net0,D,Qz,change,alpha=1):
# 计算贝叶斯网络结构发生变化后BIC评分的变化量
# 输入：net0:变化前的网络结构
# D:数据集
# Qz:隐结点分布
# change:网络结构的变化，内容为[i,j,value],
# 意为xi到xj之间的连接边变为value值
# alpha：计算BIC评分时，结构项的比重系数
# 输出：dscore:BIC评分的改变,内容为[struct,emp],
# 分别表示结构项和经验项的变化
# NewQz:新的隐结点分布
#=====网络结构的改变
i,j,value=change
consider=[i,j]
net1=net0.copy()
net1[i,j]=value
#=====隐变量子节点
son0=list(np.where(net0[:-1,-1]==-1)[0])
son1=list(np.where(net1[:-1,-1]==-1)[0])
#=====判断是否需要重新运行EM
if j==len(net0)-1 or (i in son0) or (j in son0):
Qz1=EM(net1,D,Qz.shape[1],12)
consider=consider+son0+son1+[len(net0)-1]
consider=np.unique(consider)
else:
Qz1=Qz
dscore=BIC(net1,D,Qz1,alpha,consider)-BIC(net0,D,Qz,alpha,consider)
return dscore,Qz1
#================================================
# 主程序
#================================================
#==============西瓜数据集2.0======================
# 将X和类标记Y放一起
D=[['青绿','蜷缩','浊响','清晰','凹陷','硬滑','是'],
['乌黑','蜷缩','沉闷','清晰','凹陷','硬滑','是'],
['乌黑','蜷缩','浊响','清晰','凹陷','硬滑','是'],
['青绿','蜷缩','沉闷','清晰','凹陷','硬滑','是'],
['浅白','蜷缩','浊响','清晰','凹陷','硬滑','是'],
['青绿','稍蜷','浊响','清晰','稍凹','软粘','是'],
['乌黑','稍蜷','浊响','稍糊','稍凹','软粘','是'],
['乌黑','稍蜷','浊响','清晰','稍凹','硬滑','是'],
['乌黑','稍蜷','沉闷','稍糊','稍凹','硬滑','否'],
['青绿','硬挺','清脆','清晰','平坦','软粘','否'],
['浅白','硬挺','清脆','模糊','平坦','硬滑','否'],
['浅白','蜷缩','浊响','模糊','平坦','软粘','否'],
['青绿','稍蜷','浊响','稍糊','凹陷','硬滑','否'],
['浅白','稍蜷','沉闷','稍糊','凹陷','硬滑','否'],
['乌黑','稍蜷','浊响','清晰','稍凹','软粘','否'],
['浅白','蜷缩','浊响','模糊','平坦','硬滑','否'],
['青绿','蜷缩','沉闷','稍糊','稍凹','硬滑','否']]
D=np.array(D)
FeatureName=['色泽','根蒂','敲声','纹理','脐部','触感','好瓜']
#======将“脐部”视为隐变量，对数据进行相应的修改=====
D=D[:,[0,1,2,3,5,6]] #可观测数据集
XName=['色泽','根蒂','敲声','纹理','触感','好瓜'] #x变量名称
ZName=['脐部'] #隐变量名称
ZStateNum=3 #隐变量的状态数(离散取值数目)
FeatureName=XName+ZName #包括可观测变量和隐变量的所有变量的名称
#=================初始化为朴素贝叶斯网=============
#构建贝叶斯网络，右上角元素ij表示各个连接边
#取值0表示无边，取值1表示Xi->Xj,取值-1表示Xi<-Xj
m=D.shape[0] #样本数
N=len(XName)+1 #结点数
net=np.zeros((N,N))
choose=4 #选择初始化类型，可选1,2,3,4
#分别代表全独立网络、朴素贝叶斯网络、全连接网络,随机网络
if choose==1: #全独立网络:图中没有任何连接边
pass
elif choose==2: #朴素贝叶斯网络:所有其他特征的父节点都是类标记"好瓜"
net[:-1,-1]=-1
elif choose==3: #全连接网络：任意两个节点都有连接边
again=True #若存在环，则重新生成
while again:
for i in range(N-1):
net[i,i+1:]=np.random.randint(0,2,N-i-1)*2-1
again=len(relationship(net)[2])!=0
elif choose==4: #随机网络：任意两个节点之间的连接边可有可无
again=True #若存在环，则重新生成
while again:
for i in range(N-1):
net[i,i+1:]=np.random.randint(-1,2,N-i-1)
again=len(relationship(net)[2])!=0
draw(net,FeatureName,'初始网络结构')
#==============下山法搜寻BIC下降的贝叶斯网==========
alpha=0.1 #BIC评分的结构项系数
Qz=EM(net,D,ZStateNum,12)
score0=BIC(net,D,Qz,alpha)
score=[score0]
print('===========训练贝叶斯网============')
print('初始BIC评分:%.3f(结构%.3f,经验%.3f)'%(sum(score0),score0[0],score0[1]))
eta,tao=0.1,50 #允许eta的概率调整到BIC评分增大的网络
#阈值随迭代次数指数下降
for loop in range(500):
# 随机指定需要调整的连接边的两个节点：Xi和Xj
i,j=np.random.randint(0,N,2)
while i==j:
i,j=np.random.randint(0,N,2)
i,j=(i,j) if i<j else (j,i)
# 确定需要调整的结果
value0=net[i,j] #可能为0,1,-1
change=np.random.randint(2)*2-1 #结果+1或-1
value1=(value0+1+change)%3 -1 #调整后的取值
net1=net.copy()
net1[i,j]=value1
if value1!=0 and len(relationship(net1)[2])!=0:
#如果value1取值非零，说明为转向或者增边
#若引入环，则放弃该调整
continue
chage,NewQz=BIC_change(net,D,Qz,[i,j,value1],alpha)
if sum(chage)<0 or np.random.rand()<eta*np.exp(-loop/tao):
score0=score0+chage
score.append(score0)
net=net1
Qz=NewQz
print('调整后BIC评分:%.3f(结构%.3f,经验%.3f)'
%(sum(score0),score0[0],score0[1]))
else:
continue
draw(net,FeatureName,'最终网络结构,alpha=%.3f'%(alpha))
Qz=EM(net,D,ZStateNum)
score=np.array(score)
x=np.arange(len(score))
plt.title('BIC贝叶斯网络结构搜索过程,alpha=%.3f'%(alpha))
plt.xlabel('更新次数')
plt.ylabel('分值')
plt.plot(x,score[:,0],'.r-')
plt.plot(x,score[:,1],'.b-')
plt.plot(x,score.sum(axis=1),'.k-')
plt.legend(['struct','emp','BIC(struct+emp)'])
plt.show()