[project @ 1996-07-25 21:02:03 by partain]
[nofib.git] / real / HMMS / HmmDensities.lhs
1
2         In Chapter~\ref{ch:HmmDigraphs} we were concerned only with
3 the HMM state transition structure.  In this chapter we are concerned
4 with developing and using statistical models for the observation
5 vectors.
6         \begin{haskell}{HmmDensities}
7
8 > module HmmDensities(
9 >       module Native, module MathTypes, module Phones,
10 >       GaussianComponent(..), TiedMixture(..), TmTable(..),
11 >       LogDensityTable(..),
12 >       eval_log_densities, readMixtures, readMixture, extern_to_intern
13 >       ) where
14
15 > import Native         -- hbc library modules
16 > --partain: import Maybe
17 > import Printf
18
19 > import Lists          -- general library modules
20 > import MathTypes
21 > import MaybeStateT
22
23 > import Phones         -- application-specific modules
24 > import HmmConstants
25 > import Array--1.3
26
27 \end{haskell}
28
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 \section {Multivariate Gaussian Densities}
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32
33         Let $x$ be a vector of dimension $N$ with real-valued
34 components, i.e., $x \in \reals^N$.  The formula for a multivariate
35 Gaussian probability density function having mean vector $\mu \in
36 \reals^N$ and covariance matrix $\Sigma \in \reals^{N \times N}$ is
37         \begin{equation}
38         \scriptN(x \verticalbar \mu, \Sigma)
39         = \frac{1}{ (2\pi)^{N/2}
40           (\det \Sigma)^{1/2} } \exp \{ -\frac{1}{2} (x - \mu)^{T}
41           \Sigma^{-1} (x - \mu) \}.
42         \label{eq:gauss-gen}
43         \end{equation}
44         The superscript $T$ denotes the matrix/vector transpose;
45 vectors without the superscript $T$ are taken to be column vectors.
46
47
48         If the covariance matrix is diagonal, i.e., $\Sigma =
49 \mbox{diag}\{\sigma_1^2,\sigma_2^2,\ldots,\sigma_N^2\}$, and if we
50 take $x = (x_1, x_2, \ldots, x_N)$ and $\mu = (\mu_1, \mu_2, \ldots,
51 \mu_N)$, then (\ref{eq:gauss-gen}) can be simplified to
52         \begin{equation}
53         \scriptN(x \verticalbar \mu, \Sigma) = 
54           \frac{1}{ (2\pi)^{N/2}
55                     \prod_{j=1}^N \sigma_j}
56           \exp \{ -\frac{1}{2}
57                    \sum_{j=1}^N (x_j - \mu_j)^2 / \sigma_j^2 \} .
58         \label{eq:gauss-diag}
59         \end{equation}
60
61
62         We will assume diagonal covariance matrices for the rest of
63 this chapter.
64
65
66 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67 \section {Mixture Densities}
68 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69
70         Given $M$ density functions $f_1, f_2, \ldots, f_M$ defined
71 over $\reals^N$ and $M$
72 nonnegative real numbers $c_1, c_2, \ldots, c_M$ for which
73         \[
74         \sum_{i=1}^M c_i \ = \ 1,
75         \]
76         we can define a new density function $g$ over $\reals^N$ by
77 the equation
78         \[
79         g = \sum_{i=1}^M c_i f_i.
80         \]
81         The function $g$ is called a {\em mixture density function}.
82 The $c_i$'s are the {\em mixture coefficients\/} or the {\em component
83 probabilities}, and the $f_i$'s are the {\em component densities}.
84
85
86         An excellent introduction to the theory of mixture
87 distributions is provided in~\cite{EverHand81}.
88
89
90 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 \section {Gaussian Mixture Densities}
92 \label{sc:gms}
93 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94
95         A {\em Gaussian mixture density\/} is a mixture density for
96 which each component density is Gaussian.  Thus, a Gaussian mixture
97 density has the form
98         \begin{equation}
99         g(x) = \sum^M_{i=1} c_i \scriptN(x \verticalbar \mu_i,\Sigma_i)
100         \label{eq:mixG-gen}
101         \end{equation} where the $i$th component density has mean
102 $\mu_i$ and covariance matrix $\Sigma_i$.
103
104
105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 \subsection {Representation}
107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108
109
110         Recall that we are assuming that the covariance matrix is
111 diagonal; hence, it can be represented as a vector.  The parameters of
112 a Gaussian mixture can be stored in a list of triples.  The first
113 component of each triple is the mixture coefficient, the second
114 component is the mean vector, and the third component is the variance
115 vector.
116         \begin{haskell}{GaussianComponent}
117
118 > type GaussianComponent = (Float, Vector, Vector)
119
120 \end{haskell}
121         \fixhaskellspacing\begin{haskell}{GaussianMixture}
122
123 > type GaussianMixture   = [GaussianComponent]
124
125 \end{haskell}
126         For example, a 3-component mixture density would be
127 represented by the list
128         \[
129         \startlist (c_1, \mu_1, \Sigma_1),\ 
130                    (c_2, \mu_2, \Sigma_2),\ 
131                    (c_3, \mu_3, \Sigma_3) \stoplist
132         \]
133
134
135         This is how we will store the mixture parameters in external
136 files.  However, to make the evaluation of a mixture density more
137 efficient, it pays to modify the {\em internal\/} representation
138 somewhat.
139
140
141         Using a subscript $i$ to index the mixture components and
142 another subscript $j$ to index the vector components, we can replace
143 each $\scriptN(x \verticalbar \mu_i, \Sigma_i)$ in (\ref{eq:mixG-gen})
144 by its definition (\ref{eq:gauss-diag}) to get
145         \begin{equation}
146            g(x) = \sum^M_{i=1}
147              \frac{c_i}{  (2\pi)^{N/2} \prod_{j=1}^N \sigma_{i,j} }
148             \exp \{ -\frac{1}{2} \sum_{j=1}^N (x_j - \mu_{i,j})^2 /
149             \sigma_{i,j}^2 \}.
150         \end{equation}
151         Denote the ``normalized'' mixture coefficient by $c'_i$:
152         \begin{equation}
153           c'_i = \frac{c_i}{  (2\pi)^{N/2} \prod_{j=1}^N \sigma_{i,j} }.
154         \label{eq:normalized-coefficient}
155         \end{equation}
156         Then
157         \begin{eqnarray}
158           g(x) & = & \sum^M_{i=1} c'_i \exp \{ -\frac{1}{2}
159                      \sum^N_{j=1} (x_j - \mu_{i,j})^2 / \sigma_{i,j}^2 \}\\
160                & = & \sum^M_{i=1} c'_i \exp \{ 
161                     - \sum^N_{j=1} (x_j - \mu_{i,j})^2 / (2\sigma_{i,j}^2) \}.
162         \label{eq:mixture-final-form}
163         \end{eqnarray}
164
165
166         The internal representation has the same structure as the
167 external representation, i.e., a list of triples of type
168 \verb~GaussianComponent~, except the mixture coefficient is replaced
169 by the normalized mixture coefficient and each of the variances is
170 multiplied by two.  For example, a 3-component mixture density would
171 be represented internally by the list
172         \[
173         \startlist (c'_1, \mu_1, \Sigma'_1),\ 
174                    (c'_2, \mu_2, \Sigma'_2),\ 
175                    (c'_3, \mu_3, \Sigma'_3) \stoplist
176         \]
177         where
178         \[
179         \Sigma_i = \mbox{diag}\{\sigma_{i,1}^2, \sigma_{i,2}^2,
180                                 \ldots, \sigma_{i,N}^2 \}
181         \ \Rightarrow\ 
182         \Sigma'_i = \mbox{diag}\{2\sigma_{i,1}^2, 2\sigma_{i,2}^2,
183                                 \ldots, 2\sigma_{i,N}^2 \}
184         \]
185
186
187         We need to know the dimension of the observation vectors to
188 calculate the Gaussian normalization coefficient efficiently.  We
189 could do the calculation without explicit knowledge of the observation
190 dimension by putting a $\sqrt{2\pi}$ term inside the multiplication
191 over the $\sigma_{i,j}$'s (see (\ref{eq:normalized-coefficient})), but
192 this increases the amount of computation.  Since we need the
193 observation dimension for other functions, we assume that we have
194 access to it via the module \verb~HmmConstants~.
195         \begin{haskell}{gaussian_factor}
196
197 > gaussian_factor = sqrt (2.0*pi) ^ observation_dimen :: Double
198
199 \end{haskell}
200         The normalizing calculation is performed using double
201 precision arithmetic.  The expression \verb~fromRational.toRational~
202 is a Haskell idiom for conversion of floating point types.
203         \begin{haskell}{extern_to_intern}
204
205 > extern_to_intern :: GaussianMixture -> GaussianMixture
206 > extern_to_intern = map extern_to_intern'
207
208 > extern_to_intern' :: GaussianComponent -> GaussianComponent
209 > extern_to_intern' (c, mu, vars) =
210 >       let
211 >         vars_d   = map (fromRational.toRational) vars :: [Double]
212 >         t        = gaussian_factor * sqrt (product vars_d)
213 >         denom    = (fromRational.toRational) t :: Float
214 >       in
215 >         (c/denom, mu, map (2*) vars)
216
217 \end{haskell}
218
219
220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 \subsection {Evaluation}
222 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
223
224         In the Viterbi alignment algorithm~(Chapter~\ref{ch:Viterbi}),
225 we are actually going to do our calculations on a logarithmic scale.
226 Thus, the function we really want to compute is
227         \begin{equation}
228           \ln g(x) = \ln\ \sum^M_{i=1} c'_i \exp \{ 
229                 - \sum^N_{j=1} (x_j - \mu_{i,j})^2 / (2\sigma_{i,j}^2) \}.
230         \label{eq:log-mixture-final-form}
231         \end{equation}
232
233
234         The function for evaluating (\ref{eq:log-mixture-final-form})
235 with parameters \verb~igm~ (for ``internal Gaussian mixture'') at the
236 point \verb~x~ can be defined as follows:
237         \begin{haskell}{eval_log_mixture}
238
239 > eval_log_mixture :: Vector -> GaussianMixture -> Float
240
241 > eval_log_mixture  x  igm  =
242 >       density_log (sum [c' * exp( eval_exponent x mu vars' )
243 >                            | (c', mu, vars') <- igm])
244
245 \end{haskell}
246         \fixhaskellspacing\begin{haskell}{eval_exponent}
247
248 > eval_exponent (x:rxs) (u:rus) (v:rvs) =
249 >       let  t = x - u
250 >       in   (eval_exponent rxs rus rvs) - (t*t)/v
251
252 > eval_exponent [] [] [] = 0.0
253
254 \end{haskell}
255
256
257         The use of the function \verb~density_log~ is an attempt to
258 avoid an unnecessary logarithm calculation if the density value is too
259 small.
260         \begin{haskell}{density_log}
261
262 > density_log x = if x <= density_threshold
263 >                 then smallest_log_density
264 >                 else log x
265
266 \end{haskell}
267         \fixhaskellspacing\begin{haskell}{density_threshold}
268
269 > density_threshold = 0.1023 ^ observation_dimen :: Float
270
271 \end{haskell}
272         \fixhaskellspacing\begin{haskell}{smallest_log_density}
273
274 > smallest_log_density = log density_threshold
275
276 \end{haskell}
277         We now explain how we picked the formula for the density
278 threshold.  Equation (\ref{eq:gauss-diag}) is the probability density
279 function for $N$ mutually independent Gaussian random variables.  The
280 probability that an individual zero-mean unit variance Gaussian random
281 variable takes on an absolute value exceeding 1.65 is about 0.10.  The
282 value of the density function at this point is about 0.1023.  Since
283 the Gaussian density function is unimodal, if all $N$ components had
284 taken values exceeding 1.65 in absolute value, the value of the
285 density function at such a point would have to be smaller than
286 $0.1023^N$.  But the probability of this happing is not greater than
287 $0.10^N$.  So, if our density is this small, it is quite unlikely that
288 the observed vector came from this particular mixture, so we don't
289 bother to compute the logarithm precisely.\footnote{Of course, there
290 are other ways that the density function can be that small; this
291 calculation is just an engineering approximation.}  Note that, like
292 the Gaussian normalizing factor, these constants depend on the
293 dimension of the observation vector, which is defined in the module
294 \verb~HmmConstants~; see Chapter~\ref{ch:Constants}.
295
296
297 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
298 \section {The Log-Density Table}
299 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
300
301         The states of each HMM are indexed by integers.
302         \begin{verbatim}
303
304 > type HmmState = Int
305
306 \end{verbatim}
307
308
309         Each state of a phonetic HMM either has its own mixture
310 density or it is {\em tied\/} to the density of another state,
311 possibly belonging to another HMM.
312         \begin{haskell}{TiedMixture}
313
314 > data TiedMixture = Gm GaussianMixture | Tie Phone HmmState
315
316 \end{haskell}
317
318
319         The tied mixtures are stored in a two-dimensional table that
320 is implemented as an array of arrays in order to save storage space
321 when building tables for HMMs with different numbers of states.
322         \begin{haskell}{TmTable}
323
324 > type TmTable = Array Phone (Array HmmState TiedMixture)
325
326 \end{haskell}
327
328
329         The log-density values are stored in a table that has the same
330 structure as the mixture parameter table.
331         \begin{haskell}{LogDensityTable}
332
333 > type LogDensityTable = Array Phone (Array HmmState Float)
334
335 \end{haskell}
336
337
338         The function \verb~eval_log_densities~ computes a table of
339 log-density values from an observation vector and a tied mixture
340 table.  In other words, this function evaluates the log-densities for
341 all of the states of all of the HMMs, placing the values in an array
342 for efficient retrieval.
343         \begin{haskell}{eval_log_densities}
344
345 > eval_log_densities :: TmTable -> Vector -> LogDensityTable
346
347 > eval_log_densities tmt x = ldt
348 >       where ldt = amap (amap eval_tied_mixture) tmt
349 >             eval_tied_mixture (Gm gm)   = eval_log_mixture x gm
350 >             eval_tied_mixture (Tie p k) = ldt!p!k
351
352 \end{haskell}
353
354
355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
356 \section {Reading Gaussian Mixture Densities from Binary Files}
357 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
358
359
360         The function readMixtures reads all of the mixtures in a
361 binary file that is assumed to contain only mixture data.
362         \begin{haskell}{readMixtures}
363
364 > readMixtures          :: Bytes -> [GaussianMixture]
365
366 > readMixtures bs       =
367 >       case readMixture bs of
368 >       Nothing         -> if null bs
369 >                          then []
370 >                          else error "readMixtures: left-over bytes"
371 >       Just (m, bs')   -> m : readMixtures bs'
372
373 \end{haskell}
374
375
376         The function \verb~readMixture~ reads the parameters for one
377 mixture density.  First it reads the number of component densities in
378 the mixture, then it reads each of them in.  This function and the
379 next two that follow are implemented using the Maybe State Transformer
380 (Chapter~\ref{ch:MaybeStateT}).
381         \begin{haskell}{readMixture}
382
383 > readMixture :: MST Bytes GaussianMixture
384
385 > readMixture =
386 >       readBytes                     `bindMST` \ num_components ->
387 >       readGaussians num_components  `bindMST` \ gs ->
388 >       returnMST gs
389
390 \end{haskell}
391
392
393         The first argument to the function \verb~readGaussians~ is the
394 number of component Gaussian densities in the mixture.
395         \begin{haskell}{readGaussians}
396
397 > readGaussians         :: Int -> MST Bytes GaussianMixture
398
399 > readGaussians  0      = returnMST []
400
401 > readGaussians m       = readGaussian        `bindMST` \ g ->
402 >                         readGaussians (m-1) `bindMST` \ gs ->
403 >                         returnMST (g:gs)
404
405 \end{haskell}
406
407
408         The function \verb~readGaussian~ reads the mixture
409 coefficient, the mean, and the variance vector for a multivariate
410 Gaussian density.
411         \begin{haskell}{readGaussian}
412
413 > readGaussian      :: MST Bytes GaussianComponent
414
415 > readGaussian =
416 >       readBytes                         `bindMST` \ c ->
417 >       listReadBytes observation_dimen   `bindMST` \ m ->
418 >       listReadBytes observation_dimen   `bindMST` \ v ->
419 >       returnMST (c, m, v)
420
421 \end{haskell}
422
423
424 %%%%%%%%%%  End of HmmDensities.lhs  %%%%%%%%%%