马氏链的MATLAB程序

关于马尔可夫链的MATLAB程序峄做信道分析及离散事件仿真度的同事有用。

% File: c15_bwa.m
% Software given here is to accompany the textbook: W.H. Tranter 
% K.S. Shanmugan T.S. Rappaport and K.S. Kosbar Principles of 
% Communication Systems Simulation with Wireless Applications 
% Prentice Hall PTR 2004.
%
function [p pye b] = c15_bwa（iterationsstatesout）;
len = length（out）;
p = input（‘Enter the initial state transition matrix P > ‘）;
pye = input（‘Enter the initial state probability vector pye > ‘）;
b = input（‘Enter the initial output symbol probability matrix B > ‘） ;
alpha=zeros（lenstates）; beta=zeros（lenstates）;
eta=zeros（statesstates）; gamma=zeros（1states）; scale=zeros（len1）;
loglikelihood = zeros（1iterations）;
iplot = 1;							% likelihood plot if 1
%
p									% display initial p
pye									% display initial pye
%
pye_rec = zeros（states1）;
pye_rec（:1） = pye‘;
sum_gamma = 0;
sum_eta = 0;
%
for cycle = 1:iterations
   cycle		                    % display iteration index
   %
   %     alpha generation 
   %
   alpha（1:） = pye.*b（1:）;
   scale（1） = sum（alpha（1:））;
   alpha（1:） = alpha（1:）/scale（1）;
   for t = 2:len
       alpha（t:） = （alpha（t-1:）*p）.*b（out（t）+1:）;
       scale（t） = sum（alpha（t:））;
       alpha（t:） = alpha（t:）/scale（t）;
   end
   %
   %     beta generation 
   %
   beta（len:） = 1/scale（len）;
   for t = len-1:-1:1
       beta（t:） = （beta（t+1:）.*b（out（t+1）+1:））*（p‘）/scale（t）;
   end
   %
   %     eta generation 
   %
   sum_eta = zeros（states）;
   for t = 1:len-1
      for i = 1:states
         eta（i:） = （（alpha（ti）*（p（i:）.*（b（out（t+1）+1:）））.*beta（t+1:）））...
            /sum（alpha（t:）.*beta（t:））;
       end
       sum_eta  = sum_eta + eta;
   end
   % 
   %     gamma generation 
   %
   gamma_sum = zeros（1states）;
   for t = 1:len
      gamma_sum = gamma_sum + alpha（t:）.*beta（t:）;
   end
   % 
   %		calculate and display the log_likelihood function
   %
   loglikelihood = sum（log10（scale））;
   log_likelihood（cycle） = loglikelihood;
   loglikelihood									% display result
   %
   %     Re-estimation of the intial state probability vector pye        
   %  
   pye（1:） = alpha（1:）.*beta（1:）/sum（alpha（1:）.*beta（1:））;
   %
   pye							% display pye
   %
   pye_rec（:cycle+1） = pye‘;	% Save for plot illustrating convergence	
   %
   %     Re-estimation of the state transition matrix P
   %
   for i = 1:states
      for j = 1:states
         p_estimate（ij） = sum_eta（ij）/（gamma_sum（i）-alpha（leni）.*beta（leni）...
            /（sum（alpha（len:）.*beta（len:））））;
      end
      p_estimate（i:） = p_estimate（i:）/sum（p_estimate（i:））;
   end
   %
   p = p_estimate				% display p
   %
   %     Re-estimation of output symbol probability matrix B
   %
   out_0 = find（out == 0）;
   out_1 = find（out == 1）;
   sum_0 = zeros（1states）;
   sum_1 = zeros（1states）;
   for i = 1:length（out_0）
      sum_0 = sum_0 + alpha（out_0（i）:）.*beta（out_0（i）:）...
         /sum（alpha（

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2004-08-04 15:22  Chapter_15\
     文件        3667  2004-08-02 13:29  Chapter_15\c15_bwa.m
     文件        1773  2004-08-02 13:29  Chapter_15\c15_errvector.m
     文件        1612  2004-08-02 13:29  Chapter_15\c15_hmm2.m
     文件        1011  2004-08-02 13:29  Chapter_15\c15_hmmtest.m
     文件         804  2004-08-02 13:29  Chapter_15\c15_intervals1.m
     文件        2215  2004-08-02 13:29  Chapter_15\c15_intervals2.m
     文件         702  2004-08-02 13:29  Chapter_15\c15_MMtransient.m
     文件        1170  2004-08-02 13:29  Chapter_15\c15_seglength.m
     文件        5243  2004-08-02 13:29  Chapter_15\c15_semiMarkov.m