前言：
這次練習完成的是圖模型的近似推理，參考的內容是coursera課程：Probabilistic Graphical Models . 上次實驗PGM練習四：圖模型的精確推理 中介紹的是圖模型的精確推理，但在大多數graph上，其精確推理是NP-hard的，所以有必要采用計算上可行的近似推理。本實驗中的近似推理分為2個部分，LBP(loop belief propagation算法)和MCMC采樣。實驗code可參考：實驗code可參考網友的：code.
算法流程：
LBP(loop belief propagation)算法近似推理：
(a): 輸入factor list F和Evidence E.
(b): 用evidence E掉F中的每個factor, 得到新的F.
(c): 利用F中的factor構造bethe cluster graph，該graph中的factor要么是single的要么是pairwise的.
(d): 按照某種策略選擇一條message來發射，並更新message passing前后的MESSAGE矩陣。
(e): 繼續步驟(d)直到收斂，收斂是指前后MESSAGE矩陣中的message的val值都非常接近.
(f): 步驟(e)完成后，說明bethe cluster graph已經是calibration的. 求某個變量的邊緣概率時，找到該graph中含有該變量的某個cluster，並對其求和積分掉其它變量，最后歸一化即可。
Gibbs算法圖模型的近似推理過程如下：
由此可知，圖摸型的狀態指的是所有節點的一個assignment,並不是指單獨分析某一個節點。
Metropolis-Hastings(簡稱MH算法)算法流程如下：
由上面的流程可知，MH算法與Gibbs不同之處在於，它每次更新state時是將圖中所有節點都更新，而不是一個個來。其最關鍵之處是建議分布的設計。
matlab知識：
在matlab中盡量少使用for-loop,因為速度非常慢。
[I,J] = ind2sub(siz,IND):
其中siz是一個線性遞增矩陣(從1開始按列遞增)的長和寬，IND為需要進行index的向量值。返回的I和J就是這些值所在的行和列。
ind = find(X, k):
找到矩陣X中元素不為0的最多前k個元素，ind為這k個元素的下標(linear indices ).記住find函數找到的是下標index即可。
S = sparse(i,j,s,m,n):
生成一個稀疏矩陣，矩陣中的下標分別存在向量I,j中，對應的值存在向量s中。生成的稀疏矩陣大小為m*n.
n = nnz(X):
算出矩陣X中非0元素的個數。
實驗中一些函數簡單說明：
V2F = VariableToFactorCorrespondence(V, F):
該函數的作用是找出變量所對應的factor集。其中V就是變量構成的向量，F是所有的factor構成的向量。
[toy_network, toy_factors] = ConstructToyNetwork(on_diag_weight, off_diag_weight):
該函數的作用是構造一個pair-wise的小網絡，網絡大小4×4,每個節點都只能取0和1. 網格中每條edge上的2個節點如果取值相同的話，則邊的權值大小為on_diag_weight, 反之為off_diag_weight. 返回值toy_network就是所構成的網絡結構體，里面包含的成員有names、card、edges、var2factors. 下面是一個toy network的例子：
[i, j] = NaiveGetNextClusters(P, m):
實驗1的內容。參數P為cluster graph, 里面包含的變量有clusterList, edges. 參數m表示已經傳遞了m條消息。I→j表示需要傳遞的第m+1條消息。按照j下標的升序來選edge.
[i j] = SmartGetNextClusters(P,Messages,oldMessages,m):
另一種選擇需要傳送message的方法。不過在函數內部參數Messages和oldMessages主要用於smoothing messages和informed message scheduling兩種方法，這兩種方法都是為了提高BP收斂速度的。這里沒有去實現它。課件中有對應的理論。
[i j] = GetNextClusters(P,Messages,oldMessages,m,useSmart):
用參數useSmart來選擇上面的兩種message passing方法。
P = CreateClusterGraph(F, Evidence):
實驗2的內容。用factorlist F來構造cluster graph，這比PGM練習4中構造clique tree要簡單很多，因為不要求是tree，且只需edge是兩相鄰節點交集的子集。Evidence為觀察到的變量。返回的P包含clusterList和edges兩部分。該函數內部構造的是bethe cluster tree.
converged = CheckConvergence(mNew, mOld):
實驗3的內容。該函數是用來判斷message passing過程是否收斂，如果message passing前后的message矩陣中所有元素對應的value差值都很小(e-6)時，則代表已差不多收斂。
[P MESSAGES] = ClusterGraphCalibrate(P,useSmartMP):
實驗4的內容。和PGM練習4里面的一樣，也是找到一條message，然后發送，不同的是，PGM練習4中message passing的次數固定了，而這里其次數需達到cluster graph收斂的次數才行。
M = ComputeApproxMarginalsBP(F,E):
實驗5的內容。利用LBP算法進行對變量的條件概率進行近似推理。具體過程見前面的LBP算法流程。
LogBS = BlockLogDistribution(V, G, F, A):
實驗6的內容。V為需要計算條件概率的變量集(如果在Gibbs情況下，則一般每次只取一個變量，此時V的長度就為1; 如果長度不為1,則需它們之間的值是相同的)，每個變量的取值范圍必須一樣，G為cluster graph，F為factorlist，A為當前graph的一個assignment. 該函數的作用是計算V中每個變量的log條件概率值. 計算依據如下：
由上面的公式可知，采樣Xi時只將與包含Xi的那些factor相乘(不考慮歸一化時).同時，程序中為了防止數值下溢(numerical underflow)問題,上面值的計算將在log空間進行。
[v] = randsample(vals,numSamp,replace,weightIncrements):
該函數實際上完成的是一個多項式分布采樣。其中參數vals一般表示采樣的最大值，numSamp表示需要采樣的樣本數。weightIncrements表示多項式分布的區間(因為是用0～1之間的隨機數來投票).但GibbsTrans()一直通過不了測試，課程老師說matlab 2012后的版本都不行，是個bug。但從函數內部的代碼也沒看出有和matlab版本相關的代碼，所以只要用到這個函數的，提交答案時都通過不了。下面是PGM課程團隊在新的實驗教程中對bug的提示:
A = GibbsTrans(A, G, F):
實驗7的內容。轉移矩陣采用Gibbs方法。A為圖G的一個assignment, F為factorlist.該函數的作用是從當前的A得到下一個assignment,也保持在A中。內部調用函數BlockLogDistribution()得到了其分布，然后使用randi()函數實現多項式的隨機采樣。
M = ExtractMarginalsFromSamples(G, samples, collection_indx):
參數G仍為graph, sample為采樣到的所有樣本(包括mix-time前面的樣本), collection_indx為samples中用於inference樣本的下標。輸出M為每個樣本的margin值。內部實現很簡單，其實就是求估計變量的概率，用符合要求的樣本個數除以樣本總數。
E2F = EdgeToFactorCorrespondence(V, F):
算出圖中每條edge所對應的factor編號。
[M, all_samples] = MCMCInference(G, F, E, TransName, mix_time, num_samples, sampling_interval, A0):
實驗8和實驗12的內容。MCMC的接口函數。其中的G、F、E與前面的一樣;TransName代表MCMC轉移矩陣的類型; mix_time表示分布達到穩定前的時間。num_samples為所需sampling的樣本。sampling_interval為采樣間隔，一般都設置為1. A為Markov Chain的初始狀態，即初始的assignment. 該函數輸出值M求的就是每一個變量的條件概率，是調用函數ExtractMarginalsFromSamples()來實現的。
ogp = LogProbOfJointAssignment(F, A):
求出assignment A與F中所有相關factor對應var處的乘積，就是聯合概率。這個乘積與MH算法公式:
里那個分子分母有關，如果A為當前時刻的，則logp值為公式中的分母;如果A為狀態轉移得到的下一時刻，則lopg為公式中的分子。為什么這里2個會與實驗教程公式1中相對應呢?主要是因為公式中的轉移矩陣T此時為1,直接可以約掉。另外因為公式中是相除，所以概率值不用歸一化。
A = MHUniformTrans(A, G, F):
實驗9的內容。轉移矩陣采用MH方法。同樣是從當前的assignment得到下一個assignment，只不過這里采用的是均勻分布取下一狀態. 內部可采用類似的LogProbOfJointAssignment()思想實現。
[sci sizes] = scomponents(A):
生成的sci,sizes是2個列向量。其中的sci是返回的每個頂點所對應連通區域的標號, sizes是對應連通區域中元素的個數。
A = MHSWTrans(A, G, F, variant):
實驗10和實驗11的內容，實現MH轉移概率。采用Swendsen-Wang建議分布的MH算法。該算法有2種建議分布類型, 兩者的主要不同之處體現在狀態轉移概率那里。下面是Swendsen-Wang狀態轉移的示意圖：
由圖可以看出，其轉移過程主要有下面幾個步驟：
- 對graph隨機取一個初始狀態;
- 將graph中不同值兩節點之間的邊斷開，這樣相同值的節點就構成了一個個小連通圖;
- 將上面得到的每個小連通圖中的每條邊以概率(1-qij)隨機斷開，這樣就得到了更多小小連通圖;
- 以均等概率選取圖中的一個小連通圖，並將該小圖記為Y;
- 以R概率分布對Y中節點值進行改變，不管是改變前還是改變后，Y中所有節點的值都是一樣的;
套用公式1中的MH理論，此時建議分布為：
兩個的比為：
其中記：
而練習中的2個Swendsen-Wang算法不同主要體現在qij選取不同以及R分布的不同。算法2中的qij計算公式為：
其R采用的是BlockLogDistribution.m中的方法。
相關理論知識點：
LBP算法不僅限於clique trees，它對任意的cluster graph都適應。只不過是不同的cluster graph的收斂速度和精度不同。由於cluster graph中有feedback環，所以它只能用於近似推理。在LBP算法中，message passing order並沒有嚴格要求，不用規定root節點。
如果給定樣本(IID樣本)后，則可以利用這些樣本來估計一些統計量，且由Hoeffding Bound和Chernoff Bound理論保證，當樣本數越多時，其估計越准確。Hoeffding Bound也稱為additive bound, Chernoff Bound被稱為Chernoff Bound.
如果給定統計量估計的錯誤率，則由這2個bound都可以計算出小於該錯誤率所需的最小樣本數。
兩個bound雖然存在，但是作用不大，因為有些統計量需要的樣本數呈指數增加。
對BNs進行forward sampling比較簡單(網絡參數給定的前提下)，先sampling父節點，然后sampling子節點。如果已知某些變量的觀察值，則直接將那些不符合觀測值的樣本點去掉即可。
Forward sampling不適用於Mns.
Regular Markov Chains是指任意兩個狀態之間都有連線，並且每一個狀態都有自己的一個環(self-transition).
不管初始化狀態如何，Regular Markov Chains都收斂於一個唯一的穩定分布(離散的話則是多項式分布)，對這個分布進行采樣比較簡單。
如果某個分布P的樣本不好直接采樣，則我們可以構造一個Regular Markov Chains T，使得T的穩定分布就是分布P，這樣從T上采樣就可得到P的樣本。
mixing: Regular Markov Chains達到穩定狀態時則稱為mixed. 但是很難證明Regular Markov Chains是否已mixed，頂多能驗證它有沒有mixed. 采用MC sampling方法的難點就在於對mixed狀態的判斷。一般可從多個不同初始狀態值出發，對chain采用窗口計算一些統計量，通過可視化這些統計量來判斷。
即使Regular Markov Chains已經mixed，穩定分布下連續采樣得到的樣本也是相關的，並不屬於IID，所以這些樣本不能直接用於統計推斷。按照經驗來說，如果Regular Markov Chains收斂速度越快，則這些相鄰樣本的相關性越弱，也就越趨於IID，越有用。
前面的sampling都是假設我們已經構造好了P所對應的Regular Markov Chain. 那么到底該怎樣來構造呢？對應有很多的理論。不過對於PGM中的兩種模型BNs和MNs而言，已有理論證明：對Gibbs chian的采樣過程是可行的。
如果graph中所有的factor都是positive的，則Gibbs chain是regular的。並且Gibbs chain采樣某個元素值時只與該元素所在的factor有關，與其它factor無關。當然，前提是我們知道怎樣從這些factor中來采樣。
精確推理在general networks上是NP-hard的。
Reversible Chains: 滿足細致平衡的chain. 細致平衡是指當前狀態為x，且下一狀態為x'的概率與當前狀態為x', 且下一狀態為x的概率相等。
構造分布P對應的MCMC鏈的一個通用的理論框架是MH(Metropolis-Hasngs)算法(Gibbs也屬於MH算法的一種)。
MH框架下需要一個建議分布Q，以及一個接受率A。該建議分布必須是reversible的。為了縮短mix time，Q應該越寬越好，但是如果太寬則接受率又會過低。
一般來說，當graph越稀疏就越適合用message passing的方法來inference，如果越dense，則越適合用sampling的方法。
為了使inference更准確，message passing方法中應該使cluster graph 中的cluster更大;在sampling方法中應該使建議分布每次轉移覆蓋更多變量。
參考資料：
Daphne Koller，Probabilistic Graphical Models Principles and Techniques書籍第12章
coursera課程：Probabilistic Graphical Models