前言:
本次實驗包含了2部分:貝葉斯模型參數的學習以及貝葉斯模型結構的學習,在前面的博文PGM練習七:CRF中參數的學習 中我們已經知道怎樣學習馬爾科夫模型(CRF)的參數,那個實驗采用的是優化方法,而這里貝葉斯模型參數的學習是先假定樣本符合某種分布,然后使用統計的方法去學習這些分布的參數,來達到學習模型參數的目的。實驗內容請參考 coursera課程:Probabilistic Graphical Models 中的assignmnet 8,實驗code可參考網友的:code
實驗中所用到的body pose表現形式如下:
共有10個節點,每個節點由2個坐標變量和1個方向變量構成。關節名稱和坐標系可直接參考上圖。
matlab知識:
cov(x): 求的是X無偏差的協方差,即分母除的是n-1.
cov(X,1):求的是X的最大似然估計協方差,即分母除的是n.
實驗中一些函數簡單說明:
[mu sigma] = FitGaussianParameters(X):
實驗1內容。輸入X是一個向量,該函數是計算X的均值和標准差,直接按照均值方程公式計算就OK了。
[Beta sigma] = FitLinearGaussianParameters(X, U):
實驗2內容。X: 大小為(M x 1),表示輸入的樣本向量X含有M個樣本,是樹結構中的子節點。U: 大小為(M x N), 表示有N個父節點,每個父節點也對應M個樣本。該函數主要是計算N+1個beta系數和一個sigma系數,其計算方法可參考前面的公式。給出本節點的數據和其所有父節點的數據,就可以計算出在父節點下的條件概率參數。本實驗中,如果某個節點只含有一個父節點,在使用公式計算時也相當於有3個父節點,這是因為每個節點由3個維度構成,需4個參數。
關於CLG(conditional linear Gaussian)模型有:
如果節點X有n個父節點U1..Un,則條件概率為:
求出樣本的log似然表達式,然后分別這些表達式中的n+1個參數求導,令導數為0,化簡后可得到n+1個方程,其中與X直接相關的為:
其它X分別與n個父節點相關的方程為:
利用上面n+1個方程就可以求出Beta的n+1個參數了,最后將log似然函數對sigda求導,並化簡有下面式子(怎么我推導出來不同呢):
這樣的話,CLG模型中所有的參數都可以從數據樣本中學到了。
VisualizeDataset(Dataset):
DataSet大小為n*10*3,即有n個樣本,每個樣本包含了骨架中的10個節點信息。該函數是將這n個樣本對應的骨架圖一張張顯示出來,有動畫感,大家可以跑一下, 有點意思。
val = lognormpdf(x, mu, sigma):
求高斯分布在x處的log值,高斯分布的參數為mu和sigma,可以是多維的高斯分布。
loglikelihood = ComputeLogLikelihood(P, G, dataset):
實驗3的內容。P為模型參數的結構體,P.c為長度為2的向量,其值分別代表human和alien的先驗概率。P.clg為CLG模型中的參數,包含有mu_x, mu_y, mu_angle, sigma_x, sigma_y, sigma_angle, 以及theta(theta長度為12,因為每個節點有3個坐標信息,而每個坐標在其父節點情況下需4個參數表示). 當然P.clg中應該同時包含human和alien的這些參數。參數G為父節點的描述,具體表示細節可參考實驗教材。dataset為訓練樣本,大小為N x 10 x 3. 輸出參數loglikelihood為整個dataset在模型下的log似然值。似然函數值表示模型對數據的擬合程度,可作為后續樣本分類的依據。loglikeihood的值按下列公式計算:
其中的聯合概率可以按類別的全概率公式分解計算,如下:
[P loglikelihood] = LearnCPDsGivenGraph(dataset, G, labels):
實驗4的內容。在實驗3中,我們是已知了每個節點的CLG參數,然后求在這些參數下數據的似然值的。那么這些參數到底怎樣得到呢?本函數就是完成該功能的。對每個節點的3個坐標而言,每次計算一個坐標對應的CLG參數,使用的樣本為:該數據對應的樣本以及其父節點對應樣本(包含全部3個坐標),然后調用FitLinearGaussianParameters()來計算模型參數,結果都保存在P中,另外輸出的loglikelihood是直接調用ComputeLogLikelihood()來函數實現。
accuracy = ClassifyDataset(dataset, labels, P, G):
實驗5的內容。給定樣本dataset和標簽labels,給定模型結構G,給定模型參數P,求對應樣本分類的正確率。實現時需記住:有多少個類別,相當於計算了多少個模型,計算哪個模型下產生該樣本的概率大,就將該樣本分類為對應的類別。
I = GaussianMutualInformation(X, Y):
輸入的矩陣X大小為N*D1,表示N個樣本,每個樣本D1維;矩陣Y大小為N*D2,同樣表示N個樣本,每個樣本D2維。這兩個矩陣X和Y的每一維樣本都服從高斯分布。該函數是計算兩個多維高斯分布樣本的互信息(互信息是個標量,如果2個矩陣相同,則互信息為0),計算公式如下:
adj = MaxSpanningTree (weights):
weights為N*N大小的對稱矩陣,weights(I,j)的值表示節點i和節點j之間的權重。返回的adj也是一個N*N大小的矩陣,但不是對稱矩陣。adj(I,j)=1表示有從節點i指向節點j的邊。該函數的功能就是從權值矩陣中學出tree的結構,並輸出。
[A W] = LearnGraphStructure(dataset):
實驗6的內容。該函數是從樣本矩陣dataset(沒有使用標簽信息)中學習模型的權值矩陣W和模型的樹形結構矩陣A.其中每條邊的權值計算公式如下:
該公式第一個等號的推導過程可參考Daphne Koller,Probabilistic Graphical Models Principles and Techniques書籍第18章的18.4小節。它的物理意義是:在i和j之間加入邊前后得分score的變化量,如果為正,則表示加入該邊后score增加,且說明該結構有利,score增加越大,結構就越有利。第二個等號的推導可參考coursera課程:Probabilistic Graphical Models 在structure learning課件中的第7頁,等號成立是在likelihood score前提下。
G = ConvertAtoG(A):
該函數是將max spanning tree A轉換為圖的結構矩陣G,G的含義與前面的一樣。
[P G loglikelihood] = LearnGraphAndCPDs(dataset, labels):
實驗7的內容。該函數實現的是使用dataset學習模型的結構,然后結合labels來學習模型的參數,最后利用該參數來計算樣本的log似然。
相關理論知識點:
這部分可參考Corsera中的課件以及網友demonstrate的blog:PGM 讀書筆記節選(十五)
Structure learning的好處有:可發現數據中的相關信息;當專業領域的知識不完善而建立不起完全正確的graph model時可采用。structure learning常見的有2種:constraint-based structure learning和score-based structure learning,BNs結構學習中一般采用后面那種比較多,即在預設的幾個結構中找出得分最高的那個(通過 搜索)。
Likelihood Score:指的是給定模型結構G和數據D,求出最優的參數,然后求出在G和該參數下對數據D的log似然值。其表達式如下:
式中第一個求和表示子節點與父節點之間的互信息,第2個求和為節點的熵,關於互信息計算公式如下:
互信息的值>=0,且當變量x和y相互獨立時互信息為0,但通常情況下即使x和y是從相互獨立中的分布采用得來的,其互信息也大於0。
Likelihood score容易overfitting,為了減少overfitting可采用BIC score(BIC:Bayesian Information Criterion),公式為:
該score對模型中相互獨立參數的個數進行了限制,能夠平衡模型對數據的擬合度和模型的復雜度。BIC score具有漸進一致性(consistency),因為隨着樣本的增多,它能夠收斂到與true graph的I-equivalent網絡上,並且score也越來越大。
Bayesian score是利用所有參數的加權平均,其表達式如下:
第一項叫Marginal Likelihood,包含了參數的先驗知識在里面,要求是能夠分解。第二項叫結構先驗,它對邊的數目,參數的個數進行懲罰,即對模型的復雜度進行了懲罰。在Bayesian score 中,有唯一一個導致 consistency 的參數先驗BDe prior,所以一般情況下選擇BDe prior. BIC score是Bayeisan score的一種近似。
在擁有score函數以及一些候選網絡結構后,需要用搜索算法來尋找最優網絡結構。一般可采用候選網絡結構為tree(由tree構成forest), 因為tree能夠很好的被優化,且擁有稀疏參數,泛化能力較強。另外tree中每個節點頂多有1個父節點,所以我們可以對它進行分解來計算score值:先將有父節點和無父節點的節點分開,然后根據邊的權值公式在網絡加入邊直到score不再增加,或者出現環時停止。最后采用max-weight spanning trees(MST)搜索算法找出權值大的那些邊。
對於Likelihood score來說,邊的權值非負,對於BIC score和Bayesian score來說,邊的權值可以為負數(因為score函數對邊數有懲罰,相當於邊權值為負)。
實際中很多網絡並不是tree結構,也就是說可以有多個父節點,更復雜。有理論證明,如果父節點個數>1,則搜索得分最高的網絡的過程是NP-hard的。此時一般采用啟發式搜索,以greedy爬山法為例:給定一個初始的網絡結構(可以是空的,隨機的,或是學習出來的tree,亦或是先驗設計的),在該網絡基礎上采用添加、刪除、反向等操作逐步選擇得分最高的那個操作,依次迭代直到score不再變。雖然邊的反向操作可以由刪除和添加2個步驟構成,但是由於在局部最優點附近采用greedy的方法時,刪除邊會導致score降低(意味着不會采用該操作),這樣就實現不了邊的反向了,因此需要引入Edge Reversal操作。
Naive Computational Analysis: 每一步edge的添加、刪除、反向的時間復雜度為O(n^2);如果網絡中有n個component,且計算每個component的充分統計量需O(M)時間;另外環檢測(Acyclicity check)時間和邊數m成正比,為O(m); 則結構搜索中的每一步代價為O(n^2*(Mn+m)).
錯誤提示:
matlab2013b+ubuntu13.10下矩陣運算出現如下錯誤:
BLAS loading error: dlopen: cannot load any more object matlab.
網上查了下有不少人碰見過,我這里是重啟下matlab解決的。
參考資料:
coursera課程:Probabilistic Graphical Models
網友demonstrate的blog:PGM 讀書筆記節選(十五)
Daphne Koller,Probabilistic Graphical Models Principles and Techniques書籍第18章