隐马尔可夫模型（二）--- viterbi算法

2. 根据可观察序列找到最可能的隐藏状态序列

 

需要理解的是这样的一个过程跟混淆矩阵所代表的一个过程的区别在于哪里？

混淆矩阵是我知道当前的海藻状态是dry，要得到当前可能的天气状况需要乘以混淆矩阵。

但是根据可观察序列找到最可能的隐藏状态序列不是由观察状态得到隐藏状态，而是由一个观察状态序列得到隐藏状态序列。目标是序列。

 

2.1.穷举搜索

            

可观察的序列已经给出（dry->damp->soggy）,寻找最可能的隐藏序列就是找出上面说有可能路径中概率最大的。所以需要计算每一个路径的概率。

 

需要用到HMM中的状态转移矩阵。这种方法是可行，但是计算的代价很高。

 

2.2.使用递归降低复杂度

与之前的前向算法一样，先计算部分概率，然后在计算剩下的概率。与前向算法不同的是在t时刻最可能到达的某个状态的一条路径的概率，而不是所有的概率之和。

 

每一个状态都存在一个概率最大的路径，这个路径是最优路径。称之为部分最优路径，概率称之为部分概率。对于每一个时刻的每一个状态都存在这样一个路径和这样一个概率。

最后需要的得到的是t=T时刻的每一个状态的部分最优路径和这个路劲的部分概率。选择概率最大的状态的部分最优路径便是全局的最优路径。

 

1.初始状态t=1时刻部分概率

  这个时候不存路径，所以需要初始状态概率π={πi}乘以混淆矩阵。

                                                                                

 

2.计算t>1时刻的部分概率

采用的是递归的计算方法，所以在计算t时刻的部分概率，只需要知道t-1时刻的部分概率。

                                                              

 

计算所有到状态X的概率，找到其中最可能的路径，只可能是下列三种路径之一。

1.(状态序列)，. . .，A，X

2.(状态序列)，. . .，B，X

3.(状态序列)，. . .，C，X

 

计算公式：Pr (到达A的最优路径) . Pr (X | A) . Pr (观察状态 | X)

                          

 

Pr（it-1）：t-1时刻某一状态的部分概率

Pr（X|i）：t时刻在t-1时刻状态为i的条件下，状态为x的概率。也就是HMM中的状态转移矩阵。

Pr（observationt |X）：t时刻的观察状态在隐藏状态为x条件下的概率，也就是HMM中的混淆矩阵。

 

于是得到下列公式：

                                                                     

 

2.3后向指针

通过之前的过程，得出的是最优路径的概率，但是最优的路径是怎样的，还是不清楚，我们需要一个方法得到最优路径的每一个节点，这便需要用到后向指针。

                                                              

 

δt-1 (j)：t-1时刻状态为j的最优概率

aij ：由状态t-1时刻的状态j转移到t时刻状态i的状态转移概率

Argmaxj（δt-1 (j)aij ）：找出t-1时刻的所有的状态转移到t时刻状态i的最大的值，便对应着最优路径节点。

 

 

 

2.4在matlab中使用viterbi算法

TRANS = [0.50 0.375 0.125;0.25 0.125 0.625;0.25 0.375 0.375];
EMIS = [0.60 0.20 0.15 0.05;0.25 0.25 0.25 0.25;0.05 0.10 0.35 0.50];
%seq是输出序列  
%States是状态序列
%[seq,states]=hmmgenerate(10,TRANS,EMIS);
seq = [1 1 1 1 2 3 1 1 1 1 ];
%函数hmmviterbi使用Viterbi算法计算模型给定输出序列seq最有可能通过的状态序列

 

viterbi算法解析

 

[currentState,logP] = hmmviterbi(seq,TRANS,EMIS);

numStates = size(tr,1);
checkTr = size(tr,2);
%判断是否是方阵
if checkTr ~= numStates
    error(message('stats:hmmviterbi:BadTransitions'));
end
 
% number of rows of e must be same as number of states
 
checkE = size(e,1);
%判断是否是方阵
if checkE ~= numStates
    error(message('stats:hmmviterbi:InputSizeMismatch'));
end
 
numEmissions = size(e,2);
customStatenames = false;
 
% deal with options
if nargin > 3
    okargs = {'symbols','statenames'};
    [symbols,statenames] = ...
        internal.stats.parseArgs(okargs, {[] []}, varargin{:});
    
    if ~isempty(symbols)
        numSymbolNames = numel(symbols);
        if ~isvector(symbols) || numSymbolNames ~= numEmissions
            error(message('stats:hmmviterbi:BadSymbols'));
        end
        [~, seq]  = ismember(seq,symbols);
        if any(seq(:)==0)
            error(message('stats:hmmviterbi:MissingSymbol'));
        end
    end
    if ~isempty(statenames)
        numStateNames = length(statenames);
        if numStateNames ~= numStates
            error(message('stats:hmmviterbi:BadStateNames'));
        end
        customStatenames = true;
    end
end
 
L = length(seq);
%判断输出序列是否有误，是否小于最小值，是否大于最大值（状态空间的个数），是否存在小数
if any(seq(:)<1) || any(seq(:)~=round(seq(:))) || any(seq(:)>numEmissions)
     error(message('stats:hmmviterbi:BadSequence', numEmissions));
end
%创建一个新的序列
currentState = zeros(1,L);
if L == 0
    return
end
%取对数
%为什么这里取对数，因为对数中的加法运算，相当于计算公式中的乘法运算
logTR = log(tr);
logE = log(e);
 
 
pTR = zeros(numStates,L);
%假设模型出于第0步的状态1
%创建一个矩阵，所有元素都为-Inf
v = -Inf(numStates,1);
 
v(1,1) = 0;
vOld = v;
 
%前一步状态为inner
%后一步状态法为states
for count = 1:L
    for state = 1:numStates
        % for each state we calculate
        % v(state) = e(state,seq(count))* max_k(vOld(:)*tr(k,state));
        bestVal = -inf;
        bestPTR = 0;
        % use a loop to avoid lots of calls to max
 
        for inner = 1:numStates 
            val = vOld(inner) + logTR(inner,state);
            if val > bestVal
                bestVal = val;
                %保存前一步状态，在后来的后向指针算法中要用到
                bestPTR = inner;
            end
        end
        %保存最大的状态转移概率
        % save the best transition information for later backtracking
        pTR(state,count) = bestPTR;
        % update v
        %取出混淆矩阵中状态为state，输出为seq（count）的元素
        v(state) = logE(state,seq(count)) + bestVal;
    end
    vOld = v;
end
 
%得出最优路径的概率和最终的状态
[logP, finalState] = max(v);
 
%下面要使用后向指针
currentState(L) = finalState;
for count = L-1:-1:1
    currentState(count) = pTR(currentState(count+1),count+1);
    if currentState(count) == 0
        error(message('stats:hmmviterbi:ZeroTransitionProbability', currentState( count + 1 )));
    end
end
 
if customStatenames
    currentState = reshape(statenames(currentState),1,L);
end