前言：
本次實驗主要任務是學習CRF模型的參數，實驗例子和PGM練習3中的一樣，用CRF模型來預測多張圖片所組成的單詞，我們知道在graph model的推理中，使用較多的是factor，而在graph model參數的學習中，則使用較多的是指數線性模型，本實驗的CRF使用的是log-linear模型，實驗內容請參考 coursera課程：Probabilistic Graphical Models 中的assignmnet 7. 實驗code可參考網友的:code實驗對應的模型示意圖如下：
CRF參數求解過程：
本實驗中CRF模型所表示的條件概率計算公式為:
其中的分母為划分函數，表達式為：
采用優化方法訓練CRF模型的參數時，主要任務是計算模型的cost和grad表達式。其中cost表達式為：
grad表達式為：
公式中的2個期望值表示模型對特征的期望以及數據對特征的期望，其表達式如下：
在計算cost和grad時需要分別計算下面6個中間量：
關於這幾個中間量的計算方法，可以參考實驗教程中的介紹，或者直接看博文后面貼出的代碼，這里簡單介紹下其計算方法：
log partition function: 需先建立CRF模型對應的clique tree，並對其校正，校正過程中時需要message passing，而最后passing的消息和最后一個clique factor相乘后的val之和取對數就logZ了。
weighted feature counts: 當訓練樣本中某個樣本及標簽的值符合CRF模型的某一個特征時，就將該特征對應的參數值累加，最后求和即可。
regularization cost: 直接計算。
model expected feature counts: 計算模型對特征的期望，同樣需要用到前面校正好了的clique tree. 當某個特征的變量全部屬於clique tree中某個clique變量時，求出該clique對應的factor中符合這些特征變量值的和，注意歸一化。
data feature counts: 在計算weighted feature counts的同時，如果某個特征在樣本中出現，則對相應特征計數。
regularization gradient term: 直接計算。
matlab知識：
ndx = sub2ind(siz,varargin):
siz為一個矩陣的維度向量，varargin輸入的向量表示在矩陣size的位置，返回的是linear index的值。比如sub2ind([3 4],2,4)返回11,表示在3×4大小的矩陣中，第2行第4列為矩陣的第11個元素。
C = horzcat(A1,...,AN):
將參數表示的矩陣在水平方向合成一個大矩陣。
實驗中一些函數簡單說明：
[cost, grad] = LRCostSGD(X, y, theta, lambda, i):
計算帶L2懲罰項的LR cost. 其中X是輸入矩陣，每一行代表一個樣本，y為對應的標簽向量。theta為LR模型的權值參數，lambda為權值懲罰系數。i表示選擇X矩陣中第i個樣本來計算(循環取,mod實現)。cost和grad分別為這個樣本的誤差值和輸出對權值的導數值。
thetaOpt = StochasticGradientDescent (gradFunc, theta0, maxIter):
實驗1的內容。gradFunc是函數句柄，[cost, grad] = gradFunc(theta, i),計算logistic regression在第i個樣本theta處的cost和grad. theta0為權值初始值,maxIter為最大迭代次數。這里的每次迭代只使用了一個樣本，采用隨機梯度下降法(SGD)更新權值。
thetaOpt = LRTrainSGD(X, y, lambda):
該函數完成的是用訓練樣本X和標簽y對LR進行參數優化，迭代次數和初始學習率等超參數在函數內部給定，實現該函數時需調用StochasticGradientDescent().
pred = LRPredict (X, theta):
計算樣本矩陣X在參數theta下的預測標簽pred.
acc = LRAccuracy(GroundTruth, Predictions):
計算關於真實標簽GroundTruth和預測標簽Predictions的准確率。
allAcc = LRSearchLambdaSGD(Xtrain, Ytrain, Xvalidation, Yvalidation, lambdas):
實驗2的內容。該函數是計算lambdas向量中每個lambda在驗證集Xvalidation,Yvalidation上的錯誤率，這些錯誤率保存在輸出變量allAcc中。
[P, logZ] = CliqueTreeCalibrate(P, isMax):
實驗3的內容。對clique tree P進行校正，在校正過程中同時求得划分函數Z的log值:logZ. 求logZ時只能用sum-product不能使用max-product,所以此時的isMax=0. 其方法是將最后一次傳送的message和最后一個clique相乘得到的factor，然后將factor中的val求和即可。
featureSet = GenerateAllFeatures(X, modelParams):
X是訓練樣本矩陣，因為在CRF中需要同時輸入多張圖片(本實驗中多張圖片構成一個單詞)，所以這里X中的每一行代表一張圖片。結構體modelParams有3個成員：numHiddenStates，表示CRF中隱含節點的個數，這里為26(26個字母); numObservedStates,表示CRF中觀察節點的個數，這里為2(每個像素要么為0,要么為1); lambda,權值懲罰系數。返回值featureSet包括2個成員：numParams, CRF中參數的個數，需考慮權值共享情況，即可能有多個特征共用一個權值; features, 裝有多個feature的向量，且每個feature又是一個結構體。該feature結構體中有3個成員，如下：var,特征所包含的變量;assignment,因為特征一般為指示函數，所以表示特征只在assignment處的值為1,其它處為0; paramIdx,特征所對應的參數在theta中的索引。
features = ComputeConditionedSingletonFeatures (X, modelParams):
計算輸入圖像中單個像素的特征，如果輸入X為3*32大小，因為有26個字母，所以總共的特征數為3*32*26=2496. 又假設每個像素可取0或1兩個值，所以總共的參數個數為2*32*26=1664. 很明顯有些特征是共用相同參數的。得到的features.var為圖片序列的編號，features.assignment為對應字母的編號，features.paramIdx由像素值決定參數的位置。
features = ComputeUnconditionedSingletonFeatures (len, modelParams):
len為圖片序列中圖片的個數。如果len=3,則該函數實現后的特征向量長度為3*26=78，參數個數為26. features.var為圖片編號，features.assignment為字母編號的，features.paramIdx位置和字母編號一樣。
features = ComputeUnconditionedPairFeatures (len, modelParams):
len為圖片序列中圖片的個數。如果len=3, 則該函數實現后的特征向量長度為2*26*26=1352，參數個數為26×26=676. 其中的features.var為圖片序列編號的組合，features.assignment為字母序號編號的組合，features.paramIdx為所在字母序號中對應的位置。
由上面學到的3種特征可知，特征的var都是與輸入圖片序列的標號有關，特征的assignment都是與字母的序號有關，paramIdx可能與字母序列以及圖片序列編號都 有關。實驗教程中給出的3種特征如下：
上面第2種特征為應該為conditionedPairFeatures，但和實驗code沒有對應起來，實驗code中該特征被替換成了conditionedSigleFeatures.
VisualizeCharacters (X):
可視化樣本X，因為X中一個樣本序列，可能包含多個字母，該函數將X中所含的字母顯示在一張圖上。
[nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams):
實驗4和5的內容。其中參數X,y,modelParams和前面介紹的一樣，注意X矩陣對CRF來說只算一個樣本。參數theta為列向量，大小numParams x 1，是整個CRF模型中共享的參數。這2個實驗的實現主要按照博文前面介紹的算法來計算，代碼如下：
% function [nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams) % returns the negative log-likelihood and its gradient, given a CRF with parameters theta, % on data (X, y). % % Inputs: % X Data. (numCharacters x numImageFeatures matrix) % X(:,1) is all ones, i.e., it encodes the intercept/bias term. % y Data labels. (numCharacters x 1 vector) % theta CRF weights/parameters. (numParams x 1 vector) % These are shared among the various singleton / pairwise features. % modelParams Struct with three fields: % .numHiddenStates in our case, set to 26 (26 possible characters) % .numObservedStates in our case, set to 2 (each pixel is either on or off) % .lambda the regularization parameter lambda % % Outputs: % nll Negative log-likelihood of the data. (scalar) % grad Gradient of nll with respect to theta (numParams x 1 vector) % % Copyright (C) Daphne Koller, Stanford Univerity, 2012 function [nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams) % featureSet is a struct with two fields: % .numParams - the number of parameters in the CRF (this is not numImageFeatures % nor numFeatures, because of parameter sharing) % .features - an array comprising the features in the CRF. % % Each feature is a binary indicator variable, represented by a struct % with three fields: % .var - a vector containing the variables in the scope of this feature % .assignment - the assignment that this indicator variable corresponds to % .paramIdx - the index in theta that this feature corresponds to % % For example, if we have: % % feature = struct('var', [2 3], 'assignment', [5 6], 'paramIdx', 8); % % then feature is an indicator function over X_2 and X_3, which takes on a value of 1 % if X_2 = 5 and X_3 = 6 (which would be 'e' and 'f'), and 0 otherwise. % Its contribution to the log-likelihood would be theta(8) if it's 1, and 0 otherwise. % % If you're interested in the implementation details of CRFs, % feel free to read through GenerateAllFeatures.m and the functions it calls! % For the purposes of this assignment, though, you don't % have to understand how this code works. (It's complicated.) featureSet = GenerateAllFeatures(X, modelParams); %因為擬合的是條件概率，所以需要使用X % Use the featureSet to calculate nll and grad. % This is the main part of the assignment, and it is very tricky - be careful! % You might want to code up your own numerical gradient checker to make sure % your answers are correct. % % Hint: you can use CliqueTreeCalibrate to calculate logZ effectively. % We have halfway-modified CliqueTreeCalibrate; complete our implementation % if you want to use it to compute logZ. nll = 0; grad = zeros(size(theta)); %%% % Your code here: % 計算cost ctree = CliqueTreeFromFeatrue(featureSet.features, theta, modelParams); %對整個展開的CRF對應的graph而言 [ctree,logZ] = CliqueTreeCalibrate(ctree, 0); %對tree進行校正，並求出划分函數的對數 [featureCnt,weightCnt] = WeightFeatureCnt(y, featureSet.features, theta); weightedFeatureCnt = sum(weightCnt); regCost = (modelParams.lambda/2)*(theta * theta'); nll = logZ-weightedFeatureCnt+regCost; % 計算grad mFeatureCnt = ModelFeatureCount (ctree, featureSet.features, theta);%求模型期望時，不能使用y信息 regGrad = modelParams.lambda* theta; grad = mFeatureCnt-featureCnt+regGrad; end %% 該函數實現的功能是對每個特征都建立一個factor function ctree = CliqueTreeFromFeatrue(features, theta, modelParams) n = length(features); factors = repmat(EmptyFactorStruct(),n,1); for i=1:n factors(i).var = features(i).var; factors(i).card = ones(1,length(features(i).var))*modelParams.numHiddenStates; %難道都是y變量？ factors(i).val = ones(1, prod(factors(i).card)); % 給該factor賦特征值 factors(i) = SetValueOfAssignment(factors(i), features(i).assignment, exp(theta(features(i).paramIdx))); end ctree = CreateCliqueTree(factors); end %% 該函數是求輸入樣本是否滿足各個特征，如果滿足特征i，則counts(i)=1,且weighted(i)填入相應的權值。 function [counts, weighted] = WeightFeatureCnt(y, features, theta) %這里要使用y值，因為要用y來計算指示特征 %注意特征向量的長度和參數向量的長度並不相同,因為多個特征可以共用一個參數，所以一般參數向量要短些 counts = zeros(1,length(theta)); weighted = zeros(length(theta), 1); for i = 1:length(features) feature = features(i); if all(y(feature.var)==feature.assignment) %判斷所給的y是否滿足特征所描述的 counts(feature.paramIdx) = 1; weighted(i) = theta(feature.paramIdx); end end end %% 該函數是計算模型的特征期望值，利用模型對應校正好的clique tree來計算，每個特征由其對應的clique中歸一化的val構成 function mFeatureCnt = ModelFeatureCount (ctree, features, theta) mFeatureCnt = zeros(1,length(theta)); %提前開辟空間有利於matlab運算速度 for i = 1:length(features) mIdx = features(i).paramIdx; cliqueIdx = 0; for j = 1:length(ctree.cliqueList) if all(ismember(features(i).var,ctree.cliqueList(j).var)) cliqueIdx = j; %在clique tree上找到包含第i個特征所有元素的clique break; end end eval = setdiff(ctree.cliqueList(cliqueIdx).var, features(i).var); featureFactor = FactorMarginalization(ctree.cliqueList(cliqueIdx),eval); %得到只包含該特征變量的factor idx = AssignmentToIndex(features(i).assignment,featureFactor.card); mFeatureCnt(mIdx) = mFeatureCnt(mIdx) + featureFactor.val(idx) / sum(featureFactor.val);%歸一化 end end
相關理論知識點：
learning按照模型結構是否已知，數據是否完全可以分為4類。比如HMM屬於結構已知但數據不完全那一類(因為模型中的狀態變量不能觀測)。
PGM learning的任務有：probabilistic queries(for new instance);Specific prediction task(for new instance);Knowledge discovery(for distribution).
overfitting分為參數的overfitting和結構的overfitting.
PGM中比較容易整合先驗知識。
MLE(最大似然估計)與充分統計量密切相關。
在BN中，如果有disjoint set的參數，則可將似然函數分解成局部的似然函數乘積。如果是table CPD的話，則局部似然函數又能進一步按照變量中每一維分解。
數據越少需使用越簡單的模型，這樣泛化性能才好，否則很容易過擬合。
MLE的缺點是不能很好的判斷其參數估計的可信度。比如在下列兩種情況下用拋硬幣估計硬幣朝上的概率時，使用MLE有結果：(a). 10次有7次朝上，這時估計硬幣朝上的概率為0.7;(b). 10000次有7000次朝上,硬幣朝上的概率也被判為0.7. 雖然估計的結果都為0.7,但很明顯第二種情形的估計結果更可信，第一種情形過擬合。
為了克服MLE的缺點，可采用Bayesian估計(也叫最大后驗估計)，Bayesian估計的抗噪能力更強，類似於權值懲罰。貝葉斯估計是把模型中的參數也看成是一個隨機變量，這樣在估計該參數時會引入該參數的先驗。為了達到共軛分布的目的，該先驗分布一般取Dirichlet分布。
2個變量的Dirichlet分布曲線可以在2維平面上畫出，是因為Dirichlet變量之間有和為1的約束，相當於減少了一個自由量。
在貝葉斯網絡中，如果參數在先驗中是相互獨立的，則這些參數在后驗中也是相互獨立的。
划分函數的對數對參數求導，直接把划分函數按照定義代入公式，可得求導結果為模型對特征的期望。在采用MLE估計時，可以推出loss對參數的導數為數據對特征的期望減去模型對特征的期望, MRF下的證明如下：
所以最終MRF的梯度公式為：
相對應CRF是用的log條件似然函數，其梯度計算公式為：
MRF和CRF的loss函數的優化都屬於凸優化。並且兩者計算梯度時都需要圖模型的inference,其中MRF每次迭代需inference一次，而CRF每次迭代的每個樣本處需inference一次。從這點看貌似CRF的計算量要大些。不過由於MRF需要擬合聯合概率，而CRF不需要，所以總的來說，CRF計算量要小些。
在MRF和CRF中也可以采用MAP估計，這里的先驗函數一般為高斯先驗和拉普拉斯先驗，兩者的公式分別為：
其中高斯先驗類似L2懲罰，使得很多參數在0附近(因為接近0時導數變小，懲罰力度變小)，但不一定為0,所以對應的模型很dense，而拉普拉斯先驗類似L1懲罰，使得很多參數都為0(因為其導數不變，即懲罰力度沒變)，對應模型比較sparse.
參考資料：
coursera課程：Probabilistic Graphical Models