octave_packages/statistics-1.1.3/hmmgenerate.m

   1 ## Copyright (C) 2006, 2007 Arno Onken <asnelt@asnelt.org>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{sequence}, @var{states}] =} hmmgenerate (@var{len}, @var{transprob}, @var{outprob})
  18 ## @deftypefnx {Function File} {} hmmgenerate (@dots{}, 'symbols', @var{symbols})
  19 ## @deftypefnx {Function File} {} hmmgenerate (@dots{}, 'statenames', @var{statenames})
  20 ## Generate an output sequence and hidden states of a hidden Markov model.
  21 ## The model starts in state @code{1} at step @code{0} but will not include
  22 ## step @code{0} in the generated states and sequence.
  23 ##
  24 ## @subheading Arguments
  25 ##
  26 ## @itemize @bullet
  27 ## @item
  28 ## @var{len} is the number of steps to generate. @var{sequence} and
  29 ## @var{states} will have @var{len} entries each.
  30 ##
  31 ## @item
  32 ## @var{transprob} is the matrix of transition probabilities of the states.
  33 ## @code{transprob(i, j)} is the probability of a transition to state
  34 ## @code{j} given state @code{i}.
  35 ##
  36 ## @item
  37 ## @var{outprob} is the matrix of output probabilities.
  38 ## @code{outprob(i, j)} is the probability of generating output @code{j}
  39 ## given state @code{i}.
  40 ## @end itemize
  41 ##
  42 ## @subheading Return values
  43 ##
  44 ## @itemize @bullet
  45 ## @item
  46 ## @var{sequence} is a vector of length @var{len} of the generated
  47 ## outputs. The outputs are integers ranging from @code{1} to
  48 ## @code{columns (outprob)}.
  49 ##
  50 ## @item
  51 ## @var{states} is a vector of length @var{len} of the generated hidden
  52 ## states. The states are integers ranging from @code{1} to
  53 ## @code{columns (transprob)}.
  54 ## @end itemize
  55 ##
  56 ## If @code{'symbols'} is specified, then the elements of @var{symbols} are
  57 ## used for the output sequence instead of integers ranging from @code{1} to
  58 ## @code{columns (outprob)}. @var{symbols} can be a cell array.
  59 ##
  60 ## If @code{'statenames'} is specified, then the elements of
  61 ## @var{statenames} are used for the states instead of integers ranging from
  62 ## @code{1} to @code{columns (transprob)}. @var{statenames} can be a cell
  63 ## array.
  64 ##
  65 ## @subheading Examples
  66 ##
  67 ## @example
  68 ## @group
  69 ## transprob = [0.8, 0.2; 0.4, 0.6];
  70 ## outprob = [0.2, 0.4, 0.4; 0.7, 0.2, 0.1];
  71 ## [sequence, states] = hmmgenerate (25, transprob, outprob)
  72 ## @end group
  73 ##
  74 ## @group
  75 ## symbols = @{'A', 'B', 'C'@};
  76 ## statenames = @{'One', 'Two'@};
  77 ## [sequence, states] = hmmgenerate (25, transprob, outprob,
  78 ##                      'symbols', symbols, 'statenames', statenames)
  79 ## @end group
  80 ## @end example
  81 ##
  82 ## @subheading References
  83 ##
  84 ## @enumerate
  85 ## @item
  86 ## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics
  87 ## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC,
  88 ## 2001.
  89 ##
  90 ## @item
  91 ## Lawrence R. Rabiner. A Tutorial on Hidden Markov Models and Selected
  92 ## Applications in Speech Recognition. @cite{Proceedings of the IEEE},
  93 ## 77(2), pages 257-286, February 1989.
  94 ## @end enumerate
  95 ## @end deftypefn
  96
  97 ## Author: Arno Onken <asnelt@asnelt.org>
  98 ## Description: Output sequence and hidden states of a hidden Markov model
  99
 100 function [sequence, states] = hmmgenerate (len, transprob, outprob, varargin)
 101
 102   # Check arguments
 103   if (nargin < 3 || mod (length (varargin), 2) != 0)
 104     print_usage ();
 105   endif
 106
 107   if (! isscalar (len) || len < 0 || round (len) != len)
 108     error ("hmmgenerate: len must be a non-negative scalar integer")
 109   endif
 110
 111   if (! ismatrix (transprob))
 112     error ("hmmgenerate: transprob must be a non-empty numeric matrix");
 113   endif
 114   if (! ismatrix (outprob))
 115     error ("hmmgenerate: outprob must be a non-empty numeric matrix");
 116   endif
 117
 118   # nstate is the number of states of the hidden Markov model
 119   nstate = rows (transprob);
 120   # noutput is the number of different outputs that the hidden Markov model
 121   # can generate
 122   noutput = columns (outprob);
 123
 124   # Check whether transprob and outprob are feasible for a hidden Markov
 125   # model
 126   if (columns (transprob) != nstate)
 127     error ("hmmgenerate: transprob must be a square matrix");
 128   endif
 129   if (rows (outprob) != nstate)
 130     error ("hmmgenerate: outprob must have the same number of rows as transprob");
 131   endif
 132
 133   # Flag for symbols
 134   usesym = false;
 135   # Flag for statenames
 136   usesn = false;
 137
 138   # Process varargin
 139   for i = 1:2:length (varargin)
 140     # There must be an identifier: 'symbols' or 'statenames'
 141     if (! ischar (varargin{i}))
 142       print_usage ();
 143     endif
 144     # Upper case is also fine
 145     lowerarg = lower (varargin{i});
 146     if (strcmp (lowerarg, 'symbols'))
 147       if (length (varargin{i + 1}) != noutput)
 148         error ("hmmgenerate: number of symbols does not match number of possible outputs");
 149       endif
 150       usesym = true;
 151       # Use the following argument as symbols
 152       symbols = varargin{i + 1};
 153     # The same for statenames
 154     elseif (strcmp (lowerarg, 'statenames'))
 155       if (length (varargin{i + 1}) != nstate)
 156         error ("hmmgenerate: number of statenames does not match number of states");
 157       endif
 158       usesn = true;
 159       # Use the following argument as statenames
 160       statenames = varargin{i + 1};
 161     else
 162       error ("hmmgenerate: expected 'symbols' or 'statenames' but found '%s'", varargin{i});
 163     endif
 164   endfor
 165
 166   # Each row in transprob and outprob should contain probabilities
 167   # => scale so that the sum is 1
 168   # A zero row remains zero
 169   # - for transprob
 170   s = sum (transprob, 2);
 171   s(s == 0) = 1;
 172   transprob = transprob ./ repmat (s, 1, nstate);
 173   # - for outprob
 174   s = sum (outprob, 2);
 175   s(s == 0) = 1;
 176   outprob = outprob ./ repmat (s, 1, noutput);
 177
 178   # Generate sequences of uniformly distributed random numbers between 0 and
 179   # 1
 180   # - for the state transitions
 181   transdraw = rand (1, len);
 182   # - for the outputs
 183   outdraw = rand (1, len);
 184
 185   # Generate the return vectors
 186   # They remain unchanged if the according probability row of transprob
 187   # and outprob contain, respectively, only zeros
 188   sequence = ones (1, len);
 189   states = ones (1, len);
 190
 191   if (len > 0)
 192     # Calculate cumulated probabilities backwards for easy comparison with
 193     # the generated random numbers
 194     # Cumulated probability in first column must always be 1
 195     # We might have a zero row
 196     # - for transprob
 197     transprob(:, end:-1:1) = cumsum (transprob(:, end:-1:1), 2);
 198     transprob(:, 1) = 1;
 199     # - for outprob
 200     outprob(:, end:-1:1) = cumsum (outprob(:, end:-1:1), 2);
 201     outprob(:, 1) = 1;
 202
 203     # cstate is the current state
 204     # Start in state 1 but do not include it in the states vector
 205     cstate = 1;
 206     for i = 1:len
 207       # Compare the randon number i of transdraw to the cumulated
 208       # probability of the state transition and set the transition
 209       # accordingly
 210       states(i) = sum (transdraw(i) <= transprob(cstate, :));
 211       cstate = states(i);
 212     endfor
 213
 214     # Compare the random numbers of outdraw to the cumulated probabilities
 215     # of the outputs and set the sequence vector accordingly
 216     sequence = sum (repmat (outdraw, noutput, 1) <= outprob(states, :)', 1);
 217
 218     # Transform default matrices into symbols/statenames if requested
 219     if (usesym)
 220       sequence = reshape (symbols(sequence), 1, len);
 221     endif
 222     if (usesn)
 223       states = reshape (statenames(states), 1, len);
 224     endif
 225   endif
 226
 227 endfunction
 228
 229 %!test
 230 %! len = 25;
 231 %! transprob = [0.8, 0.2; 0.4, 0.6];
 232 %! outprob = [0.2, 0.4, 0.4; 0.7, 0.2, 0.1];
 233 %! [sequence, states] = hmmgenerate (len, transprob, outprob);
 234 %! assert (length (sequence), len);
 235 %! assert (length (states), len);
 236 %! assert (min (sequence) >= 1);
 237 %! assert (max (sequence) <= columns (outprob));
 238 %! assert (min (states) >= 1);
 239 %! assert (max (states) <= rows (transprob));
 240
 241 %!test
 242 %! len = 25;
 243 %! transprob = [0.8, 0.2; 0.4, 0.6];
 244 %! outprob = [0.2, 0.4, 0.4; 0.7, 0.2, 0.1];
 245 %! symbols = {'A', 'B', 'C'};
 246 %! statenames = {'One', 'Two'};
 247 %! [sequence, states] = hmmgenerate (len, transprob, outprob, 'symbols', symbols, 'statenames', statenames);
 248 %! assert (length (sequence), len);
 249 %! assert (length (states), len);
 250 %! assert (strcmp (sequence, 'A') + strcmp (sequence, 'B') + strcmp (sequence, 'C') == ones (1, len));
 251 %! assert (strcmp (states, 'One') + strcmp (states, 'Two') == ones (1, len));