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.
8 > module HmmDensities(
9 >       module Native, module MathTypes, module Phones,
10 >       GaussianComponent(..), TiedMixture(..), TmTable(..),
11 >       LogDensityTable(..),
13 >       ) where
15 > import Native         -- hbc library modules
16 > --partain: import Maybe
17 > import Printf
19 > import Lists          -- general library modules
20 > import MathTypes
21 > import MaybeStateT
23 > import Phones         -- application-specific modules
24 > import HmmConstants
25 > import Array--1.3
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 \section {Multivariate Gaussian Densities}
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
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}
62         We will assume diagonal covariance matrices for the rest of
63 this chapter.
66 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67 \section {Mixture Densities}
68 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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}.
86         An excellent introduction to the theory of mixture
87 distributions is provided in~\cite{EverHand81}.
90 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 \section {Gaussian Mixture Densities}
92 \label{sc:gms}
93 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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$.
105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 \subsection {Representation}
107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
118 > type GaussianComponent = (Float, Vector, Vector)
123 > type GaussianMixture   = [GaussianComponent]
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$
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.
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}
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$
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~.
197 > gaussian_factor = sqrt (2.0*pi) ^ observation_dimen :: Double
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.
205 > extern_to_intern :: GaussianMixture -> GaussianMixture
206 > extern_to_intern = map extern_to_intern'
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)
220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 \subsection {Evaluation}
222 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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}
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:
239 > eval_log_mixture :: Vector -> GaussianMixture -> Float
241 > eval_log_mixture  x  igm  =
242 >       density_log (sum [c' * exp( eval_exponent x mu vars' )
243 >                            | (c', mu, vars') <- igm])
248 > eval_exponent (x:rxs) (u:rus) (v:rvs) =
249 >       let  t = x - u
250 >       in   (eval_exponent rxs rus rvs) - (t*t)/v
252 > eval_exponent [] [] [] = 0.0
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.
262 > density_log x = if x <= density_threshold
263 >                 then smallest_log_density
264 >                 else log x
269 > density_threshold = 0.1023 ^ observation_dimen :: Float
274 > smallest_log_density = log density_threshold
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}.
297 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
298 \section {The Log-Density Table}
299 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
301         The states of each HMM are indexed by integers.
302         \begin{verbatim}
304 > type HmmState = Int
306 \end{verbatim}
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.
314 > data TiedMixture = Gm GaussianMixture | Tie Phone HmmState
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.
324 > type TmTable = Array Phone (Array HmmState TiedMixture)
329         The log-density values are stored in a table that has the same
330 structure as the mixture parameter table.
333 > type LogDensityTable = Array Phone (Array HmmState Float)
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.
345 > eval_log_densities :: TmTable -> Vector -> LogDensityTable
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
355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
356 \section {Reading Gaussian Mixture Densities from Binary Files}
357 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
361 binary file that is assumed to contain only mixture data.
364 > readMixtures          :: Bytes -> [GaussianMixture]
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'
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}).
383 > readMixture :: MST Bytes GaussianMixture
386 >       readBytes                     bindMST \ num_components ->
387 >       readGaussians num_components  bindMST \ gs ->
388 >       returnMST gs
393         The first argument to the function \verb~readGaussians~ is the
394 number of component Gaussian densities in the mixture.
397 > readGaussians         :: Int -> MST Bytes GaussianMixture
399 > readGaussians  0      = returnMST []
401 > readGaussians m       = readGaussian        bindMST \ g ->
402 >                         readGaussians (m-1) bindMST \ gs ->
403 >                         returnMST (g:gs)
416 >       readBytes                         bindMST \ c ->
417 >       listReadBytes observation_dimen   bindMST \ m ->
418 >       listReadBytes observation_dimen   bindMST \ v ->