前言:html
此次練習完成的是圖模型的近似推理,參考的內容是coursera課程:Probabilistic Graphical Models . 上次實驗PGM練習四:圖模型的精確推理 中介紹的是圖模型的精確推理,但在大多數graph上,其精確推理是NP-hard的,因此有必要採用計算上可行的近似推理。本實驗中的近似推理分爲2個部分,LBP(loop belief propagation算法)和MCMC採樣。實驗code可參考:實驗code可參考網友的:code.git
算法流程:github
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的.dom
(d): 按照某種策略選擇一條message來發射,並更新message passing先後的MESSAGE矩陣。ide
(e): 繼續步驟(d)直到收斂,收斂是指先後MESSAGE矩陣中的message的val值都很是接近.函數
(f): 步驟(e)完成後,說明bethe cluster graph已是calibration的. 求某個變量的邊緣機率時,找到該graph中含有該變量的某個cluster,並對其求和積分掉其它變量,最後歸一化便可。tornado
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狀態轉移的示意圖:
由圖能夠看出,其轉移過程主要有下面幾個步驟:
套用公式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章