--- /dev/null
+## Copyright (C) 2006, 2007 Arno Onken <asnelt@asnelt.org>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{transprobest}, @var{outprobest}] =} hmmestimate (@var{sequence}, @var{states})
+## @deftypefnx {Function File} {} hmmestimate (@dots{}, 'statenames', @var{statenames})
+## @deftypefnx {Function File} {} hmmestimate (@dots{}, 'symbols', @var{symbols})
+## @deftypefnx {Function File} {} hmmestimate (@dots{}, 'pseudotransitions', @var{pseudotransitions})
+## @deftypefnx {Function File} {} hmmestimate (@dots{}, 'pseudoemissions', @var{pseudoemissions})
+## Estimate the matrix of transition probabilities and the matrix of output
+## probabilities of a given sequence of outputs and states generated by a
+## hidden Markov model. The model assumes that the generation starts in
+## state @code{1} at step @code{0} but does not include step @code{0} in the
+## generated states and sequence.
+##
+## @subheading Arguments
+##
+## @itemize @bullet
+## @item
+## @var{sequence} is a vector of a sequence of given outputs. The outputs
+## must be integers ranging from @code{1} to the number of outputs of the
+## hidden Markov model.
+##
+## @item
+## @var{states} is a vector of the same length as @var{sequence} of given
+## states. The states must be integers ranging from @code{1} to the number
+## of states of the hidden Markov model.
+## @end itemize
+##
+## @subheading Return values
+##
+## @itemize @bullet
+## @item
+## @var{transprobest} is the matrix of the estimated transition
+## probabilities of the states. @code{transprobest(i, j)} is the estimated
+## probability of a transition to state @code{j} given state @code{i}.
+##
+## @item
+## @var{outprobest} is the matrix of the estimated output probabilities.
+## @code{outprobest(i, j)} is the estimated probability of generating
+## output @code{j} given state @code{i}.
+## @end itemize
+##
+## If @code{'symbols'} is specified, then @var{sequence} is expected to be a
+## sequence of the elements of @var{symbols} instead of integers.
+## @var{symbols} can be a cell array.
+##
+## If @code{'statenames'} is specified, then @var{states} is expected to be
+## a sequence of the elements of @var{statenames} instead of integers.
+## @var{statenames} can be a cell array.
+##
+## If @code{'pseudotransitions'} is specified then the integer matrix
+## @var{pseudotransitions} is used as an initial number of counted
+## transitions. @code{pseudotransitions(i, j)} is the initial number of
+## counted transitions from state @code{i} to state @code{j}.
+## @var{transprobest} will have the same size as @var{pseudotransitions}.
+## Use this if you have transitions that are very unlikely to occur.
+##
+## If @code{'pseudoemissions'} is specified then the integer matrix
+## @var{pseudoemissions} is used as an initial number of counted outputs.
+## @code{pseudoemissions(i, j)} is the initial number of counted outputs
+## @code{j} given state @code{i}. If @code{'pseudoemissions'} is also
+## specified then the number of rows of @var{pseudoemissions} must be the
+## same as the number of rows of @var{pseudotransitions}. @var{outprobest}
+## will have the same size as @var{pseudoemissions}. Use this if you have
+## outputs or states that are very unlikely to occur.
+##
+## @subheading Examples
+##
+## @example
+## @group
+## transprob = [0.8, 0.2; 0.4, 0.6];
+## outprob = [0.2, 0.4, 0.4; 0.7, 0.2, 0.1];
+## [sequence, states] = hmmgenerate (25, transprob, outprob);
+## [transprobest, outprobest] = hmmestimate (sequence, states)
+## @end group
+##
+## @group
+## symbols = @{'A', 'B', 'C'@};
+## statenames = @{'One', 'Two'@};
+## [sequence, states] = hmmgenerate (25, transprob, outprob,
+## 'symbols', symbols, 'statenames', statenames);
+## [transprobest, outprobest] = hmmestimate (sequence, states,
+## 'symbols', symbols,
+## 'statenames', statenames)
+## @end group
+##
+## @group
+## pseudotransitions = [8, 2; 4, 6];
+## pseudoemissions = [2, 4, 4; 7, 2, 1];
+## [sequence, states] = hmmgenerate (25, transprob, outprob);
+## [transprobest, outprobest] = hmmestimate (sequence, states, 'pseudotransitions', pseudotransitions, 'pseudoemissions', pseudoemissions)
+## @end group
+## @end example
+##
+## @subheading References
+##
+## @enumerate
+## @item
+## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics
+## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC,
+## 2001.
+##
+## @item
+## Lawrence R. Rabiner. A Tutorial on Hidden Markov Models and Selected
+## Applications in Speech Recognition. @cite{Proceedings of the IEEE},
+## 77(2), pages 257-286, February 1989.
+## @end enumerate
+## @end deftypefn
+
+## Author: Arno Onken <asnelt@asnelt.org>
+## Description: Estimation of a hidden Markov model for a given sequence
+
+function [transprobest, outprobest] = hmmestimate (sequence, states, varargin)
+
+ # Check arguments
+ if (nargin < 2 || mod (length (varargin), 2) != 0)
+ print_usage ();
+ endif
+
+ len = length (sequence);
+ if (length (states) != len)
+ error ("hmmestimate: sequence and states must have equal length");
+ endif
+
+ # Flag for symbols
+ usesym = false;
+ # Flag for statenames
+ usesn = false;
+
+ # Variables for return values
+ transprobest = [];
+ outprobest = [];
+
+ # Process varargin
+ for i = 1:2:length (varargin)
+ # There must be an identifier: 'symbols', 'statenames',
+ # 'pseudotransitions' or 'pseudoemissions'
+ if (! ischar (varargin{i}))
+ print_usage ();
+ endif
+ # Upper case is also fine
+ lowerarg = lower (varargin{i});
+ if (strcmp (lowerarg, 'symbols'))
+ usesym = true;
+ # Use the following argument as symbols
+ symbols = varargin{i + 1};
+ # The same for statenames
+ elseif (strcmp (lowerarg, 'statenames'))
+ usesn = true;
+ # Use the following argument as statenames
+ statenames = varargin{i + 1};
+ elseif (strcmp (lowerarg, 'pseudotransitions'))
+ # Use the following argument as an initial count for transitions
+ transprobest = varargin{i + 1};
+ if (! ismatrix (transprobest))
+ error ("hmmestimate: pseudotransitions must be a non-empty numeric matrix");
+ endif
+ if (rows (transprobest) != columns (transprobest))
+ error ("hmmestimate: pseudotransitions must be a square matrix");
+ endif
+ elseif (strcmp (lowerarg, 'pseudoemissions'))
+ # Use the following argument as an initial count for outputs
+ outprobest = varargin{i + 1};
+ if (! ismatrix (outprobest))
+ error ("hmmestimate: pseudoemissions must be a non-empty numeric matrix");
+ endif
+ else
+ error ("hmmestimate: expected 'symbols', 'statenames', 'pseudotransitions' or 'pseudoemissions' but found '%s'", varargin{i});
+ endif
+ endfor
+
+ # Transform sequence from symbols to integers if necessary
+ if (usesym)
+ # sequenceint is used to build the transformed sequence
+ sequenceint = zeros (1, len);
+ for i = 1:length (symbols)
+ # Search for symbols(i) in the sequence, isequal will have 1 at
+ # corresponding indices; i is the right integer for that symbol
+ isequal = ismember (sequence, symbols(i));
+ # We do not want to change sequenceint if the symbol appears a second
+ # time in symbols
+ if (any ((sequenceint == 0) & (isequal == 1)))
+ isequal *= i;
+ sequenceint += isequal;
+ endif
+ endfor
+ if (! all (sequenceint))
+ index = max ((sequenceint == 0) .* (1:len));
+ error (["hmmestimate: sequence(" int2str (index) ") not in symbols"]);
+ endif
+ sequence = sequenceint;
+ else
+ if (! isvector (sequence))
+ error ("hmmestimate: sequence must be a non-empty vector");
+ endif
+ if (! all (ismember (sequence, 1:max (sequence))))
+ index = max ((ismember (sequence, 1:max (sequence)) == 0) .* (1:len));
+ error (["hmmestimate: sequence(" int2str (index) ") not feasible"]);
+ endif
+ endif
+
+ # Transform states from statenames to integers if necessary
+ if (usesn)
+ # statesint is used to build the transformed states
+ statesint = zeros (1, len);
+ for i = 1:length (statenames)
+ # Search for statenames(i) in states, isequal will have 1 at
+ # corresponding indices; i is the right integer for that statename
+ isequal = ismember (states, statenames(i));
+ # We do not want to change statesint if the statename appears a second
+ # time in statenames
+ if (any ((statesint == 0) & (isequal == 1)))
+ isequal *= i;
+ statesint += isequal;
+ endif
+ endfor
+ if (! all (statesint))
+ index = max ((statesint == 0) .* (1:len));
+ error (["hmmestimate: states(" int2str (index) ") not in statenames"]);
+ endif
+ states = statesint;
+ else
+ if (! isvector (states))
+ error ("hmmestimate: states must be a non-empty vector");
+ endif
+ if (! all (ismember (states, 1:max (states))))
+ index = max ((ismember (states, 1:max (states)) == 0) .* (1:len));
+ error (["hmmestimate: states(" int2str (index) ") not feasible"]);
+ endif
+ endif
+
+ # Estimate the number of different states as the max of states
+ nstate = max (states);
+ # Estimate the number of different outputs as the max of sequence
+ noutput = max (sequence);
+
+ # transprobest is empty if pseudotransitions is not specified
+ if (isempty (transprobest))
+ # outprobest is not empty if pseudoemissions is specified
+ if (! isempty (outprobest))
+ if (nstate > rows (outprobest))
+ error ("hmmestimate: not enough rows in pseudoemissions");
+ endif
+ # The number of states is specified by pseudoemissions
+ nstate = rows (outprobest);
+ endif
+ transprobest = zeros (nstate, nstate);
+ else
+ if (nstate > rows (transprobest))
+ error ("hmmestimate: not enough rows in pseudotransitions");
+ endif
+ # The number of states is given by pseudotransitions
+ nstate = rows (transprobest);
+ endif
+
+ # outprobest is empty if pseudoemissions is not specified
+ if (isempty (outprobest))
+ outprobest = zeros (nstate, noutput);
+ else
+ if (noutput > columns (outprobest))
+ error ("hmmestimate: not enough columns in pseudoemissions");
+ endif
+ # Number of outputs is specified by pseudoemissions
+ noutput = columns (outprobest);
+ if (rows (outprobest) != nstate)
+ error ("hmmestimate: pseudoemissions must have the same number of rows as pseudotransitions");
+ endif
+ endif
+
+ # Assume that the model started in state 1
+ cstate = 1;
+ for i = 1:len
+ # Count the number of transitions for each state pair
+ transprobest(cstate, states(i)) ++;
+ cstate = states (i);
+ # Count the number of outputs for each state output pair
+ outprobest(cstate, sequence(i)) ++;
+ endfor
+
+ # transprobest and outprobest contain counted numbers
+ # Each row in transprobest and outprobest should contain estimated
+ # probabilities
+ # => scale so that the sum is 1
+ # A zero row remains zero
+ # - for transprobest
+ s = sum (transprobest, 2);
+ s(s == 0) = 1;
+ transprobest = transprobest ./ (s * ones (1, nstate));
+ # - for outprobest
+ s = sum (outprobest, 2);
+ s(s == 0) = 1;
+ outprobest = outprobest ./ (s * ones (1, noutput));
+
+endfunction
+
+%!test
+%! sequence = [1, 2, 1, 1, 1, 2, 2, 1, 2, 3, 3, 3, 3, 2, 3, 1, 1, 1, 1, 3, 3, 2, 3, 1, 3];
+%! states = [1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1];
+%! [transprobest, outprobest] = hmmestimate (sequence, states);
+%! expectedtransprob = [0.88889, 0.11111; 0.28571, 0.71429];
+%! expectedoutprob = [0.16667, 0.33333, 0.50000; 1.00000, 0.00000, 0.00000];
+%! assert (transprobest, expectedtransprob, 0.001);
+%! assert (outprobest, expectedoutprob, 0.001);
+
+%!test
+%! sequence = {'A', 'B', 'A', 'A', 'A', 'B', 'B', 'A', 'B', 'C', 'C', 'C', 'C', 'B', 'C', 'A', 'A', 'A', 'A', 'C', 'C', 'B', 'C', 'A', 'C'};
+%! states = {'One', 'One', 'Two', 'Two', 'Two', 'One', 'One', 'One', 'One', 'One', 'One', 'One', 'One', 'One', 'One', 'Two', 'Two', 'Two', 'Two', 'One', 'One', 'One', 'One', 'One', 'One'};
+%! symbols = {'A', 'B', 'C'};
+%! statenames = {'One', 'Two'};
+%! [transprobest, outprobest] = hmmestimate (sequence, states, 'symbols', symbols, 'statenames', statenames);
+%! expectedtransprob = [0.88889, 0.11111; 0.28571, 0.71429];
+%! expectedoutprob = [0.16667, 0.33333, 0.50000; 1.00000, 0.00000, 0.00000];
+%! assert (transprobest, expectedtransprob, 0.001);
+%! assert (outprobest, expectedoutprob, 0.001);
+
+%!test
+%! sequence = [1, 2, 1, 1, 1, 2, 2, 1, 2, 3, 3, 3, 3, 2, 3, 1, 1, 1, 1, 3, 3, 2, 3, 1, 3];
+%! states = [1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1];
+%! pseudotransitions = [8, 2; 4, 6];
+%! pseudoemissions = [2, 4, 4; 7, 2, 1];
+%! [transprobest, outprobest] = hmmestimate (sequence, states, 'pseudotransitions', pseudotransitions, 'pseudoemissions', pseudoemissions);
+%! expectedtransprob = [0.85714, 0.14286; 0.35294, 0.64706];
+%! expectedoutprob = [0.178571, 0.357143, 0.464286; 0.823529, 0.117647, 0.058824];
+%! assert (transprobest, expectedtransprob, 0.001);
+%! assert (outprobest, expectedoutprob, 0.001);
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff