hmm前向算法的matlab,隐马尔可夫模型(HMM)的MATLAB实现——前向后向算法

由于研究需要,在网上找了不少关于隐马尔可夫模型的MATLAB程序,可能是没有耐下心去看,总之感觉看懂别人写的程序很费劲,所以就自己动手写了一下。

主要的参考书目是李航的《统计学习方法》,并实现了书中部分例题和习题的解答,关于隐马尔可夫模型的理论知识不再赘述,这方面的经典论文和网上的优秀博客很多,这里介绍两篇经典的论文:

1,Bilmes L,et al. A Gentle Tutorial of EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models(1997)

2,Rabiner L. A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition(1989)

下面按照HMM的三个问题:概率计算,参数学习和解码顺序介绍我的代码。

1,概率计算问题--前向后向算法

function [P1,P2,p1] = ForwardBackwardAlgo(A,B,Pi,O)

% 函数功能:前向后向算法,计算alpha值矩阵P1和beta值矩阵P2

%

% 参考文献:李航《统计学习方法》

%

% 思路:

% 给定模型lambda和观测序列O,在时刻t处于状态qi的概率,记做

% gamma_t(i) = P(it = qi|O,lambda) = P(it = qi,O|lambda) / P(O|lambda)

% P(it = qi,O|lambda) = alpha_t(i) * beta_t(i)

% 在这里,alpha矩阵表示前向概率矩阵,放在矩阵P1中,且为i*t阶

% beta矩阵表示后向概率矩阵,放在矩阵P2中,且为i*t阶

%

% 输入:隐马尔可夫模型lambda(A,B,Pi),观测序列O;

% 其中A为转移概率,B为观测概率(行数表示状态集合中的元的个数,列数表示观测集中的元的个数);

% Pi为初始概率,P为1*2阶矩阵,存放所求概率P(it = qi,O|lambda)中it和qi的角标t和i,即P=[t,i]

% 【!!注意:初始概率Pi、观测序列O以及矩阵P均须为列矩阵】

% 输出:P1:alpha值构成的矩阵,为i*t阶

% P2:beta值构成的矩阵,为i*t阶

% p1:观测序列概率P(O|lambda)

%

% 示例:(李航《统计学虚方法》P188习题10.1)

% 输入:

% A1 =

% 0.5000 0.1000 0.4000

% 0.3000 0.5000 0.2000

% 0.2000 0.2000 0.6000

% B1 =

% 0.5000 0.5000

% 0.4000 0.6000

% 0.7000 0.3000

% Pi =

% 0.2000

% 0.3000

% 0.5000

% O =

% 1

% 2

% 1

% 2

% 输出

% P1 =

% 0.1000 0.0770 0.0419 0.0211

% 0.1600 0.1104 0.0355 0.0252

% 0.2800 0.0606 0.0528 0.0138

% P2 =

% 0.1125 0.2461 0.4600 1.0000

% 0.1217 0.2312 0.5100 1.0000

% 0.1049 0.2577 0.4300 1.0000

% p1 =

% 0.0601

alpha_1 = zeros();

sum_alpha_j = zeros();

sum_beta_i = zeros();

% sum_p = zeros();

P1 = zeros();

A_size = size(A);

O_size = size(O);

N = A_size(1,1);

M = A_size(1,2);

K = O_size(1,1);

P2 = zeros(N,K);

% ---------------前向算法---------------

%求解初值

for i = 1:N

alpha_1(i,1) = Pi(i,1) * B(i,O(1,1));

end

P1= alpha_1;

%递推求解前向概率,存放在矩阵P1中

for k = 1:K-1

for j = 1:M

for i = 1:N

sum_alpha_j(i,1) = P1(i,k) .* A(i,j);

end

P1(j,k+1) = sum(sum_alpha_j) .* B(j,O(k+1,1));

end

end

%计算观测概率P(O|lambda)

p1 = sum(P1(:,K));

% ----------------后向算法---------------

% 对最终时刻所有状态有beta_T(i) = 1

P2(:,K) = 1;

% 递推求解后向概率,存放在矩阵P2中

for k = K-1:-1:1

for i = 1:M

for j = 1:N

sum_beta_i(j,1) = A(i,j) .* B(j,O(k+1,1)) .* P2(j,k+1);

end

P2(i,k) = sum(sum_beta_i);

end

end上述程序中的例子用的是后向概率的算法,程序中前向算法和后向算法可以单独拿出来进行计算,另外,李航《统计学习方法》中P189的习题10.2也可以计算得到,其代码中调用上述程序函数,代码如下:function [p1,p2] = ForwardBackwardAppli(A,B,Pi,O,P)

% 函数功能:前向后向算法应用,计算P(it = qi|O,lambda)

%

% 参考文献:李航《统计学习方法》

%

% 思路:

% 给定模型lambda和观测序列O,在时刻t处于状态qi的概率,记做

% gamma_t(i) = P(it = qi|O,lambda) = P(it = qi,O|lambda) / P(O|lambda)

% P(it = qi,O|lambda) = alpha_t(i) * beta_t(i)

% 在这里,alpha矩阵表示前向概率矩阵,放在矩阵P1中,且为i*t阶

% beta矩阵表示后向概率矩阵,放在矩阵P2中,且为i*t阶

%

% 输入:隐马尔可夫模型lambda(A,B,Pi),观测序列O;

% 其中A为转移概率,B为观测概率(行数表示状态集合中的元的个数,列数表示观测集中的元的个数),Pi为初始概率;

% P为1*2阶矩阵,存放所求概率P(it = qi,O|lambda)中it和qi的角标t和i,即P=[t,i]

% 【!!注意:初始概率Pi、观测序列O以及矩阵P均须为列矩阵】

% 输出:p1:观测序列概率P(O|lambda)

% p2:所求概率P(it = qi|O,lambda)

%

% 示例:(李航《统计学习方法》P189习题10.2)

% 输入:

% A1 =

% 0.5000 0.1000 0.4000

% 0.3000 0.5000 0.2000

% 0.2000 0.2000 0.6000

% B1 =

% 0.5000 0.5000

% 0.4000 0.6000

% 0.7000 0.3000

% Pi =

% 0.2000

% 0.3000

% 0.5000

% O3 =

% 1

% 2

% 1

% 1

% 2

% 1

% 2

% 2

% P =

% 4

% 3

% 输出

% [a7,b7] = ForwardBackwardAlgo(A1,B1,Pi,O3,P)

% a7 =

% 0.0035

% b7 =

% 0.5370

P1 = zeros();

A_size = size(A);

O_size = size(O);

N = A_size(1,1);

M = A_size(1,2);

K = O_size(1,1);

P2 = zeros(N,K);

% 调用函数前向后向算法ForwardBackwardAlgo_1,计算得到alpha矩阵P1,beta矩阵P2观测序列O的概率p1

[P1,P2,p1] = ForwardBackwardAlgo(A,B,Pi,O);

% 计算所求值

p2 = P1(P(2,1),P(1,1)) * P2(P(2,1),P(1,1)) / p1;

下一篇介绍参数学习的Baum-Welch算法程序。

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
hmm算法matlab实现和实例 hmm_em.m function [LL, prior, transmat, obsmat, nrIterations] = ... dhmm_em(data, prior, transmat, obsmat, varargin) % LEARN_DHMM Find the ML/MAP parameters of an HMM with discrete outputs using EM. % [ll_trace, prior, transmat, obsmat, iterNr] = learn_dhmm(data, prior0, transmat0, obsmat0, ...) % % Notation: Q(t) = hidden state, Y(t) = observation % % INPUTS: % data{ex} or data(ex,:) if all sequences have the same length % prior(i) % transmat(i,j) % obsmat(i,o) % % Optional parameters may be passed as 'param_name', param_value pairs. % Parameter names are shown below; default values in [] - if none, argument is mandatory. % % 'max_iter' - max number of EM iterations [10] % 'thresh' - convergence threshold [1e-4] % 'verbose' - if 1, print out loglik at every iteration [1] % 'obs_prior_weight' - weight to apply to uniform dirichlet prior on observation matrix [0] % % To clamp some of the parameters, so learning does not change them: % 'adj_prior' - if 0, do not change prior [1] % 'adj_trans' - if 0, do not change transmat [1] % 'adj_obs' - if 0, do not change obsmat [1] % % Modified by Herbert Jaeger so xi are not computed individually % but only their sum (over time) as xi_summed; this is the only way how they are used % and it saves a lot of memory. [max_iter, thresh, verbose, obs_prior_weight, adj_prior, adj_trans, adj_obs] = ... process_options(varargin, 'max_iter', 10, 'thresh', 1e-4, 'verbose', 1, ... 'obs_prior_weight', 0, 'adj_prior', 1, 'adj_trans', 1, 'adj_obs', 1); previous_loglik = -inf; loglik = 0; converged = 0; num_iter = 1; LL = []; if ~iscell(data) data = num2cell(data, 2); % each row gets its own cell end while (num_iter <= max_iter) & ~converged % E step [loglik, exp_num_trans, exp_num_visits1, exp_num_emit] = ... compute_ess_dhmm(prior, transmat, obsmat, data, obs_prior_weight); % M step if adj_prior prior = normalise(exp_num_visits1); end if adj_trans & ~isempty(exp_num_trans) transmat = mk_stochastic(exp_num_trans); end if adj_obs obsmat = mk_stochastic(exp_num_emit); end
HMM是一种常见的统计模型,用于描述随机生成的序列。Viterbi算法是一种在HMM中进行解码的动态规划算法,用于寻找最可能的隐藏状态序列。下面是HMMMATLAB实现,包括Viterbi算法。 首先,我们需要定义一个HMM模型。我们假设这个模型有3个隐藏状态和2个可见状态。我们使用矩阵A表示隐藏状态之间的转移概率,矩阵B表示每个隐藏状态生成每个可见状态的概率,向量pi表示初始隐藏状态的概率分布。 ```matlab % 定义HMM模型参数 A = [0.5 0.2 0.3; 0.3 0.5 0.2; 0.2 0.3 0.5]; B = [0.5 0.5; 0.4 0.6; 0.7 0.3]; pi = [0.2 0.4 0.4]; ``` 接下来,我们需要生成一个可见序列。我们使用HMM模型中的随机过程生成一个长度为10的可见序列。 ```matlab % 生成可见序列 T = 10; q = zeros(1, T); o = zeros(1, T); q(1) = randsrc(1, 1, [1:3; pi]); o(1) = randsrc(1, 1, [1:2; B(q(1), :)]); for t = 2:T q(t) = randsrc(1, 1, [1:3; A(q(t-1), :)]); o(t) = randsrc(1, 1, [1:2; B(q(t), :)]); end ``` 接下来,我们使用Viterbi算法解码这个可见序列,得到最可能的隐藏状态序列。我们定义一个矩阵V表示每个时间步的最大概率,以及一个矩阵path表示每个时间步的最大概率对应的一个状态。 ```matlab % Viterbi算法解码 V = zeros(3, T); path = zeros(3, T); V(:, 1) = pi' .* B(:, o(1)); for t = 2:T for j = 1:3 [V(j, t), path(j, t)] = max(V(:, t-1) .* A(:, j)); V(j, t) = V(j, t) * B(j, o(t)); end end ``` 最后,我们找到最可能的隐藏状态序列。我们首先找到最后一个时间步的最大概率对应的隐藏状态,然后从后往依次寻找每个时间步的最大概率对应的隐藏状态,最终得到整个隐藏状态序列。 ```matlab % 找到最可能的隐藏状态序列 [~, q(T)] = max(V(:, T)); for t = T-1:-1:1 q(t) = path(q(t+1), t+1); end ``` 完整代码如下: ```matlab % 定义HMM模型参数 A = [0.5 0.2 0.3; 0.3 0.5 0.2; 0.2 0.3 0.5]; B = [0.5 0.5; 0.4 0.6; 0.7 0.3]; pi = [0.2 0.4 0.4]; % 生成可见序列 T = 10; q = zeros(1, T); o = zeros(1, T); q(1) = randsrc(1, 1, [1:3; pi]); o(1) = randsrc(1, 1, [1:2; B(q(1), :)]); for t = 2:T q(t) = randsrc(1, 1, [1:3; A(q(t-1), :)]); o(t) = randsrc(1, 1, [1:2; B(q(t), :)]); end % Viterbi算法解码 V = zeros(3, T); path = zeros(3, T); V(:, 1) = pi' .* B(:, o(1)); for t = 2:T for j = 1:3 [V(j, t), path(j, t)] = max(V(:, t-1) .* A(:, j)); V(j, t) = V(j, t) * B(j, o(t)); end end % 找到最可能的隐藏状态序列 [~, q(T)] = max(V(:, T)); for t = T-1:-1:1 q(t) = path(q(t+1), t+1); end ```
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值