西瓜書習題試答-第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()