1 ## Copyright (C) 2006, 2007 Arno Onken <asnelt@asnelt.org>
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
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
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/>.
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.
24 ## @subheading Arguments
28 ## @var{len} is the number of steps to generate. @var{sequence} and
29 ## @var{states} will have @var{len} entries each.
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}.
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}.
42 ## @subheading Return values
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)}.
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)}.
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.
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
65 ## @subheading Examples
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)
75 ## symbols = @{'A', 'B', 'C'@};
76 ## statenames = @{'One', 'Two'@};
77 ## [sequence, states] = hmmgenerate (25, transprob, outprob,
78 ## 'symbols', symbols, 'statenames', statenames)
82 ## @subheading References
86 ## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics
87 ## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC,
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.
97 ## Author: Arno Onken <asnelt@asnelt.org>
98 ## Description: Output sequence and hidden states of a hidden Markov model
100 function [sequence, states] = hmmgenerate (len, transprob, outprob, varargin)
103 if (nargin < 3 || mod (length (varargin), 2) != 0)
107 if (! isscalar (len) || len < 0 || round (len) != len)
108 error ("hmmgenerate: len must be a non-negative scalar integer")
111 if (! ismatrix (transprob))
112 error ("hmmgenerate: transprob must be a non-empty numeric matrix");
114 if (! ismatrix (outprob))
115 error ("hmmgenerate: outprob must be a non-empty numeric matrix");
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
122 noutput = columns (outprob);
124 # Check whether transprob and outprob are feasible for a hidden Markov
126 if (columns (transprob) != nstate)
127 error ("hmmgenerate: transprob must be a square matrix");
129 if (rows (outprob) != nstate)
130 error ("hmmgenerate: outprob must have the same number of rows as transprob");
135 # Flag for statenames
139 for i = 1:2:length (varargin)
140 # There must be an identifier: 'symbols' or 'statenames'
141 if (! ischar (varargin{i}))
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");
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");
159 # Use the following argument as statenames
160 statenames = varargin{i + 1};
162 error ("hmmgenerate: expected 'symbols' or 'statenames' but found '%s'", varargin{i});
166 # Each row in transprob and outprob should contain probabilities
167 # => scale so that the sum is 1
168 # A zero row remains zero
170 s = sum (transprob, 2);
172 transprob = transprob ./ repmat (s, 1, nstate);
174 s = sum (outprob, 2);
176 outprob = outprob ./ repmat (s, 1, noutput);
178 # Generate sequences of uniformly distributed random numbers between 0 and
180 # - for the state transitions
181 transdraw = rand (1, len);
183 outdraw = rand (1, len);
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);
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
197 transprob(:, end:-1:1) = cumsum (transprob(:, end:-1:1), 2);
200 outprob(:, end:-1:1) = cumsum (outprob(:, end:-1:1), 2);
203 # cstate is the current state
204 # Start in state 1 but do not include it in the states vector
207 # Compare the randon number i of transdraw to the cumulated
208 # probability of the state transition and set the transition
210 states(i) = sum (transdraw(i) <= transprob(cstate, :));
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);
218 # Transform default matrices into symbols/statenames if requested
220 sequence = reshape (symbols(sequence), 1, len);
223 states = reshape (statenames(states), 1, len);
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));
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));