西瓜书习题试答-第07章-贝叶斯分类器

时间 2020-09-26

标签西瓜习题贝叶分类器繁體版

原文原文链接

试答系列：“西瓜书”-周志华《机器学习》习题试答
系列目录
[第01章：绪论]
[第02章：模型评估与选择]
[第03章：线性模型]
[第04章：决策树]
[第05章：神经网络]
[第06章：支持向量机]
第07章：贝叶斯分类器
 第08章：集成学习
 第09章：聚类
 第10章：降维与度量学习
 第11章：特征选择与稀疏学习
第12章：计算学习理论(暂缺)
第13章：半监督学习
 第14章：几率图模型
(后续章节更新中...)html

7.1 试使用极大似然法估算西瓜数据集3.0中前三个属性的类条件几率。

7.2 试证实：条件独立性假设不成立时，朴素贝叶斯分类器仍有可能产生最优贝叶斯分类器。

答：能够参考“阅读材料”部分的解释。(待完善)python

一、对分类任务来讲，只需各种别的条件几率排序正确、无须精准几率值便可致使正确分类结果。
二、若属性间依赖对全部类别影响相同，或依赖关系的影响能相互抵消，则属性条件独立性假设在下降计算开销的同时不会对性能产生负面影响。面试

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} \]

经过最大化对数似然能够求得参数估计：app

\[\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\)。
那么，该贝叶斯分类器的决策函数为：dom

\[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准则构建一个贝叶斯网。

答：

一、贝叶斯网络结构在计算机中的表示。

首先，咱们约定一下如何在计算机中表达一个贝叶斯网络：设有\(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\)的链接边。好比，将下图中的贝叶斯网络用矩阵表达为：

二、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} \]

三、搜索算法。

咱们的目标是搜寻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分值进行比较。

四、计算结果。

下面是在西瓜数据集2.0上某一次运行的结果：

初始BIC评分：277.572(结构203.991,经验73.580)
最终BIC评分：119.360(结构 43.915,经验75.445)
屡次运行，在不一样初始网络的状况下，结果均表现出一个规律：“色泽”是一个独立的属性，是不是“好瓜”主要由“纹理”和“触感”两个属性决定。

7.10 以西瓜数据集2.0中属性“脐部”为隐变量，试基于EM算法构建一个贝叶斯网。

答：

一、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()

王书海等, BIC评分贝叶斯网络模型及其应用, 计算机工程, 2008.08 ↩︎ ↩︎