机器学习&数据挖掘笔记_PGM练习七:CRF中参数的学习



前言:

  本次实验主要任务是学习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.

 

  参考资料:

       实验code可参考网友的:code

       PGM练习三:马尔科夫网络在OCR上的简单应用

       coursera课程:Probabilistic Graphical Models 

 

作者:tornadomeet
出处:http://www.cnblogs.com/tornadomeet 欢迎转载或分享,但请务必声明文章出处。
新浪微博:tornadomeet,欢迎交流!


0

下一篇:机器学习&数据挖掘笔记_PGM练习六:制定决策
[按键盘← → 左右方向键也可以翻页哦~]