试答系列：“西瓜书”-周志华《机器学习》习题试答
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时，可以将计算结果与教材计算结果进行对比，对比发现，教材中对于连续型属性-----密度和含糖率的正态分布参数估计中，对于方差的估计采用了无偏估计，亦即。
7.4 实践中使用式(7.15)决定分类类别时，若数据的维度非常高，则概率连乘)
的结果通常会非常接近于0从而导致下溢。试述防止下溢的可能方案。
答：通常采用取对数的方法将连乘变为连加：
7.5 试证明：二分类任务中两类数据满足高斯分布且方差相同时，线性判别分析产生贝叶斯最优分类器。
答：假设数据满足高斯分布：，模型中需要确定的参数有均值
和
，以及共同的方差
，
为对称正定矩阵。数据集的对数似然为：
通过最大化对数似然可以求得参数估计：
上式中求取梯度时应用了关系
。
那么，该贝叶斯分类器的决策函数为：
对于二分类任务，这等价于：
上式中第二行采取了同先验假设，亦即P(0)=P(1)=1/2.
在3.4节线性判别分析(LDA)中，关于的定义与上面求得的
完全相同，而根据(3.33)式可知，
。在3.4节中求得最优投影直线方向为
， LDA对于新数据的分类是根据投影点距离两个投影中心的距离远近决定的，可以将其表达为：
注意到上式第二行右边项，而
为对称正定矩阵，因此该项恒为负，因此，上式可以进一步化简为：
对比可知，决策函数和
完全相同，因此可以说LDA产生了最优Bayes分类。
7.6 试编程实现AODE分类器，并以西瓜数据集3.0为训练集，对p.151的“测1”样本进行判别。
答：编程代码附后。
在实现AODE分类器过程中，需要应用(7.23)~(7.25)式计算各个概率值。由于西瓜数据集3.0中含有连续属性，需要对(7.24)和(7.25)计算式进行相应的拓展。
- 对于(7.24)式，若
为连续属性，则按下式计算：
其中为正态分布的概率密度函数，参数为
。
- 对于(7.25)式，存在多种情况：
(a). 若为连续属性，此时：
(b). 若
(c). 若
需要注意的是，对于上面(a),(b)两种情况下，求取正态分布参数时，可能因为子集中样本数太少而无法估计。在编程中，若样本数为零，则假定为正态分布，亦即
；若仅一个样本，则将平均值设为该唯一取值，方差设为1.
7.7 给定d个二值属性的二分类任务，假设对于任何先验概率项的估算至少需30个样例，则在朴素贝叶斯分类器式(7.15)中估算先验概率项)
需30×2=60个样例。试估计在AODE式(7.23)中估算先验概率)
所需的样例数(分别考虑最好和最坏情形)。
答：需要计算的先验概率的数目共有2×2×d个，因为c有2个取值，i有d个取值，
有2个取值。
首先考虑属性，对于参数
，分别需要30个样例来估算，总计需要120个样例。
在最好情形下，对于所有其他属性，假设刚好能在c=0的60个样例中找到30个
=0和30个
=1的样例，同样地，对于c=1的情况也如此，则无需更多样例，此时对于所有属性总计需要120个样例；
在最坏情形下，比如对于属性，假设60个c=0的样例中
取值全部为0，则需要另外补充30个
的样例，同样的，对于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不成立。
证明：1、在同父结构中。假设x3⊥x4成立，则有：
另一方面有：
两式相减有：
然而上式未必成立，因此x3⊥x4成立未必成立。
2、在顺序结构中。首先有联合概率为：
可以证明：
上式即表明 z⊥y | x.
至于z⊥y是否成立，先假设其成立，则有：
又因为：
两式相减有：
然而上式未必成立，因此z⊥y未必成立。
7.9 以西瓜数据集2.0为训练集，试基于BIC准则构建一个贝叶斯网。
答：首先，我们约定一下如何在计算机中表达一个贝叶斯网络：设有个节点，
，利用
矩阵的右上角元素来表达连接边，右上角元素
可能取值为0,1,-1，取值0意味着
与
之间没有连接边，等于1意味着
的连接边，等于-1意味着
的连接边。比如，将下图中的贝叶斯网络用矩阵表达为：
贝叶斯网络: 图形表达和矩阵表达
本题的解答主要参考了文献：王书海等, BIC评分贝叶斯网络模型及其应用, 计算机工程, 2008.08
在贝叶斯网络中，对于某个节点，设它的可能取值数目为
，它的父节点集为
,父节点取值的可能组合数目为
。在数据集D中，父节点集
取值组合为第
种组合的样本数目为
，与此同时，节点
取第
个值的样本数为
。那么，对于某个网络结构B，它的BIC评分函数计算公式可以表示为：
上式的等式右边第一项可以视为结构风险项，第二项为经验风险项。第一项中是由于
个概率之和为1，因此独立参数个数为
。
我们的目标是搜寻BIC评分最小化的贝叶斯网络结构，这里采用下山法进行搜索，每一步随机选取一个边进行调整，调整包括三种操作：加边、减边、转向。需要注意的是，加边和转向后可能会引入环，要予以避免。
三种基本调整操作
如果新网络的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算法构建一个贝叶斯网。
答：
附：编程代码
习题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])
def show_parameter(net,D,name=None,laplace=True):
# 计算并输出贝叶斯网络参数，也就是各个概率值，方法为简单计数
# 输入：
# net:贝叶斯网络结构
# D:数据集
# name:指定各个节点的名称，默认为None，用x1,x2...xN作为名称
# laplace:计算概率时是否考虑laplace修正
m,N=D.shape #样本数和特征数
values=[np.unique(D[:,i])for i in range(len(D[0]))] #各个离散属性的可能取值
if name==None:
name=['x%d'%(i+1) for i in range(N)] #特征名称
parents=[list(np.where(net[i,:]==-1)[0]+1)+ #父节点
list(np.where(net[:,i]==1)[0]+1)
for i in range(N)]
for i in range(N):
print('{:^20}'.format(name[i])) #首先显示所考察特征名称
r=len(values[i]) #Xi结点的取值数
pa=parents[i] #Xi的父节点编号
nums=[len(values[p-1]) for p in pa] #父节点取值数
q=np.prod(nums) #父节点取值组合数
#-----如果父节点数目为零，亦即没有父节点
if len(pa)==0:
line_side='-' #边框行
line_value='|' #取值行
line_probility='|' #概率行
for value_k in values[i]:
Dk=D[D[:,i]==value_k] #D中Xi取值v_k的子集
mk=len(Dk) #Dk子集大小
line_side+="{:-^7}".format('')+'-'
line_value+="{:^5}".format(value_k)+'|'
line_probility+="{:^7.2f}".format((mk+laplace)/(m+laplace*r))+'|'
print('\n'.join([line_side,line_value,line_side,
line_probility,line_side]))
continue
#-----有父节点时，通过编码方式，遍历所有父节点取值组合
line_side='-' #边框行
line_title='|' #标题行
for p in pa:
line_side+="{:-^7}".format('')+'-'
line_title+="{:^5}".format(name[p-1])+'|'
for value_k in values[i]:
line_side+="{:-^7}".format('')+'-'
line_title+="{:^5}".format(value_k)+'|'
print('\n'.join([line_side,line_title,line_side]))
line_probility=['|']*q #取值概率行
for code in range(q):
#解码：比如，父节点有2×3种组合，
#将code=0,1,2,3,4,5解码为
#[0,0],[0,1],[0,2],[1,0],[1,1],[1,2]
code0=code
decode=[]
for step in range(len(nums)-1):
wight=np.prod(nums[step+1:])
decode.append(code0//wight)
code0=code0%wight
decode.append(code0)
#父节点取某一组合时的子集
index=range(m) #子集索引号，初始为全集D
condition=''
for j in range(len(pa)):
value=values[pa[j]-1][decode[j]]
indexj=np.where(D[:,pa[j]-1]==value)[0]
index=np.intersect1d(index,indexj)
line_probility[code]+="{:^5}".format(value)+'|'
Dij=D[index,:] #子集
mij=len(Dij) #子集大小
for value_k in values[i]:
Dijk=Dij[Dij[:,i]==value_k] #Dij中Xi取值v_k的子集
mijk=len(Dijk) #Dijk子集大小
line_probility[code]+="{:^7.2f}".format(
(mijk+laplace)/(mij+laplace*r))+'|'
show=[]
for qq in range(q):
show.append(line_probility[qq])
show.append(line_side)
print('\n'.join(show))
#================================================
# 主程序
#================================================
#==============西瓜数据集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)
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)
if 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
print('调整后BIC评分:%.3f(结构%.3f,经验%.3f)'
%(sum(score0),score0[0],score0[1]))
net=net1
if sum(delta_score)>0:
eta=np.sqrt(eta)
else:
continue
draw(net,FeatureName,'最终网络结构')
show_parameter(net,D,FeatureName)
