labman2.tex

一本关于dsp的原著教材

likelihood of an observation sequence $X=\{x_1,x_2,\cdots,x_T\}$ with
respect to a Hidden Markov Model with parameters $\Theta$ expands as
follows\,:
\[
	p(X|\Theta) = \sum_{\mbox{\scriptsize every possible }Q} p(X,Q|\Theta)
\]
i.e. it is the sum of the joint likelihoods of the sequence over all
possible state sequence allowed by the model.
%
\item[-] the {\em forward recursion}\,: in practice, the enumeration of
every possible state sequence is infeasible. Nevertheless, $p(X|\Theta)$
can be computed in a recursive way by the {\em forward recursion}. This
algorithm defines a forward variable $\alpha_t(i)$ corresponding to\,:
\[
	\alpha_t(i) = p(x_1,x_2,\cdots x_t,q^t=q_i|\Theta)
\]
i.e. $\alpha_t(i)$ is the probability of having observed the partial
sequence $\{x_1,x_2,\cdots,x_t\}$ {\em and} being in the state $i$ at time
$t$ (event denoted $q_i^t$ in the course), given the parameters
$\Theta$. For a HMM with 5 states (where states 1 and N are the
non-emitting initial and final states, and states $2 \cdots N-1$ are
emitting), $\alpha_t(i)$ can be computed recursively as follows\,:

\pagebreak
\noindent \fbox{\bf The Forward Recursion}
\begin{enumerate}
\item {\bf Initialization}
%
	\[
		\alpha_1(i) = a_{1i} \cdot b_i(x_1), \;\;\;\; 2 \leq i \leq N-1
	\]
%
where $a_{1i}$ are the transitions from the initial non-emitting state to
the emitting states with pdfs $b_{i,\,i = 2 \cdots N-1}(x)$. Note that
$b_1(x)$ and $b_N{x}$ do not exist since they correspond to the
non-emitting initial and final states.
\item {\bf Recursion}
	\[
		\alpha_{t+1}(j) = \left[ \sum_{i=2}^{N-1} \alpha_{t}(i) \cdot a_{ij} \right] b_j(x_{t+1}),
		\;\;\;\; \begin{array}{l} 1 \leq t \leq T \\ 2 \leq j \leq N-1 \end{array}
	\]
\item {\bf Termination}
	\[
		p(X|\Theta) = \left[ \sum_{i=2}^{N-1} \alpha_{T}(i) \cdot a_{iN} \right]
	\]
%
i.e. at the end of the observation sequence, sum the probabilities of the
paths converging to the final state (state number $N$).
\end{enumerate}
(For more detail about the forward procedure, refer to~\cite{RAB93},
chap.6.4.1).

This procedure raises a very important implementation issue. As a matter of
fact, the computation of the $\alpha_t$ vector consists in products of a
large number of values that are less than 1 (in general, {\em
significantly} less than 1). Hence, after a few observations ($t \approx$
10), the values of $\alpha_t$ head exponentially to 0, and the floating
point arithmetic precision is exceeded (even in the case of double
precision arithmetics). Two solutions exist for that problem. One consists
in scaling the values and undo the scaling at the end of the procedure\,:
see \cite{RAB93} for more explanations. The other solution consists in
using log-likelihoods and log-probabilities, and to compute $\log
p(X|\Theta)$ instead of $p(X|\Theta)$.
\end{itemize}

\subsubsection*{Questions\,:}
\begin{enumerate}
\item The following formula can be used to compute the log of a sum given
the logs of the sum's arguments\,:
\[
	\log(a+b) = f(\log a,\log b) = \log a + \log \left( 1 + e^{(\log b - \log a)} \right)
\]
%
Demonstrate its validity.

Naturally, one has the choice between using $\log(a+b) = \log a + \log
\left( 1 + e^{(\log b - \log a)} \right)$ or $\log(a+b) = \log b + \log
\left( 1 + e^{(\log a - \log b)} \right)$, which are equivalent in theory.
If $\log a > \log b$, which version leads to the most precise
implementation ?
\item Express the log version of the forward recursion. (Don't fully
develop the log of the sum in the recursion step, just call it
``logsum''\,: $\sum_{i=1}^{N} x_i \stackrel{\log}{\longmapsto} \mbox{logsum}_{i=1}^{N} ( \log
x_i )$.) In addition to the arithmetic precision issues, what are the other
computational advantages of the log version ?
\end{enumerate}


\subsubsection*{Answers\,:}
\expl{
\begin{enumerate}
\item Demonstration\,:
\vspace{-1em}
\[ a = e^{\log a} \;\;\;\;\;;\;\;\;\;\; b = e^{\log b} \]
\vspace{-2em}
\begin{eqnarray}
a+b & = & e^{\log a} + e^{\log b} \nonumber \\
    & = & e^{\log a} \left( 1 + e^{(\log b - \log a)} \right) \nonumber
\end{eqnarray}
\vspace{-1.5em}
\[ \log(a+b) = \log a + \log \left( 1 + e^{(\log b - \log a)} \right) \;\;\; \mbox{QED.}\]

The computation of the exponential overflows the double precision
arithmetics for big values ($\approx700$) earlier than for small
values. Similarly, the implementations of the exponential operation are
generally more precise for small values than for big values (since an error
on the input term is exponentially amplified). Hence, if $\log a > \log b$,
the first version ($\log(a+b) = \log a + \log \left( 1 + e^{(\log b - \log
a)} \right)$) is more precise since in this case $(\log b - \log a)$ is
small. If $\log a < \log b$, it is better to swap the terms (i.e. to use
the second version).
%
\item
\begin{enumerate}
\item {\bf Initialization}
%
	\[
		\alpha_1^{(log)}(i) = \log a_{1i} + \log b_i(x_1), \;\;\;\; 2 \leq i \leq N-1
	\]
%
\item {\bf Recursion}
	\[
		\alpha_{t+1}^{(log)}(j) = \left[ \mbox{logsum}_{i=2}^{N-1} \left(
							\alpha_{t}^{(log)}(i) + \log a_{ij}
						\right) \right] + \log b_j(x_{t+1}),
		\;\;\;\; \begin{array}{l} 1 \leq t \leq T \\ 2 \leq j \leq N-1 \end{array}
	\]
\item {\bf Termination}
	\[
		\log p(X|\Theta) = \left[ \mbox{logsum}_{i=2}^{N-1} \left(
			\alpha_{T}^{(log)}(i) + \log a_{iN} \right) \right]
	\]
%
\end{enumerate}
In addition to the precision issues, this version transforms the products
into sums, which is more computationally efficient. Furthermore, if the
emission probabilities are Gaussians, the computation of the
log-likelihoods $\log(b_j(x_t))$ eliminates the computation of the
Gaussians' exponential (see lab session 4).
\end{enumerate}

These two points just show you that once the theoretic barrier is crossed
in the study of a particular statistical model, the importance of the
implementation issues must not be neglected.  }

\vspace{-1em}
%%%%%%%%%
\subsection{Bayesian classification}
\label{sub:bayesian}
%%%%%%%%%
%\vspace{-0.5em}
\subsubsection*{Question\,:}
\vspace{-0.5em}
The forward recursion allows us to compute the likelihood of an observation
sequence with respect to a HMM. Hence, given a sequence of features, we are
able to find the most likely generative model in a Maximum Likelihood
sense. What additional quantities and assumptions do we need to perform a
true Bayesian classification rather than a Maximum Likelihood
classification of the sequences ?

Which additional condition makes the result of Bayesian classification
equivalent to the result of ML classification ?
\vspace{-0.5em}
\subsubsection*{Answer\,:}
\vspace{-0.5em}
\expl{
%
To perform a Bayesian classification, we need the prior probabilities
$P(\Theta_i|\Theta)$ of each model. In addition, we can assume that all the
observation sequences are equi-probable\,:
\begin{eqnarray}
P(\Theta_i|X,\Theta) & = & \frac{p(X|\Theta_i,\Theta)
					P(\Theta_i|\Theta)}{P(X|\Theta)} \nonumber \\
 & \propto & p(X|\Theta_i) P(\Theta_i) \nonumber
\end{eqnarray}
$P(\Theta_i)$ can be determined by counting the probability of occurrence
of each model (word or phoneme) in a database covering the vocabulary to
recognize (see lab session 4).

\tab If every model has the same prior probability, then Bayesian
classification becomes equivalent to ML classification.  }

\pagebreak
%%%%%%%%%
\subsection{Maximum Likelihood classification}
\label{sub:class}
%%%%%%%%%
In practice, for speech recognition, it is very often assumed that all the
model priors are equal (i.e. that the words or phonemes to recognize have
equal probabilities of occurring in the observed speech). Hence, the speech
recognition task consists mostly in performing the Maximum Likelihood
classification of acoustic feature sequences. For that purpose, we must
have of a set of HMMs that model the acoustic sequences corresponding to a
set of phonemes or a set of words. These models can be considered as
``stochastic templates''. Then, we associate a new sequence to the most
likely generative model. This part is called the {\em decoding} of the
acoustic feature sequences.

\subsubsection*{Experiment\,:}
Classify the sequences $X_1, X2, \cdots X_6$, given in the file
\com{data.mat}, in a maximum likelihood sense with respect to the six
Markov models given in table~\ref{tab:models}. Use the function
\com{logfwd} to compute the log-forward recursion expressed in the previous
section. Store the results in a matrix (they will be used in the next
section) and note them in the table below.

\subsubsection*{Example\,:}
\mat{plot(X1(:,1),X1(:,2));}
\mat{logProb(1,1) = logfwd(X1,hmm1)}
\mat{logProb(1,2) = logfwd(X1,hmm2)}
etc. \\
\mat{logProb(3,2) = logfwd(X3,hmm2)}
etc.

\medskip
\noindent Filling the \com{logProb} matrix can be done automatically with
the help of loops\,:
\begin{verbatim}
>> for i=1:6,
  for j=1:6,
    stri = num2str(i);
    strj = num2str(j);
    eval([ 'logProb(' , stri , ',' , strj , ')=logfwd(X' , stri , ',hmm' , strj , ');' ]); 
  end;
end;
>> logProb
\end{verbatim}

\noindent
\begin{center}
\renewcommand{\arraystretch}{1.5}
\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline
 Sequence &
\small $\log p(X|\Theta_1)$ & \small $\log p(X|\Theta_2)$ &
\small $\log p(X|\Theta_3)$ & \small $\log p(X|\Theta_4)$ &
\small $\log p(X|\Theta_5)$ & \small $\log p(X|\Theta_6)$ &
\parbox[c][3em][c]{11ex}{Most likely\\ model} \\ \hline
X1 & & & & & & & \\ \hline
X2 & & & & & & & \\ \hline
X3 & & & & & & & \\ \hline
X4 & & & & & & & \\ \hline
X5 & & & & & & & \\ \hline
X6 & & & & & & & \\ \hline
\end{tabular}
\end{center}

\subsubsection*{Answer\,:}
\expl{$X_1 \rightarrow HMM1$, $X_2 \rightarrow HMM3$, $X_3 \rightarrow
HMM5$, $X_4 \rightarrow HMM4$, $X_5 \rightarrow HMM6$, $X_6 \rightarrow
HMM2$.}


\pagebreak
%%%%%%%%%
%%%%%%%%%
\section{Optimal state sequence}
\label{sec:optimal}
%%%%%%%%%
%%%%%%%%%
\subsubsection*{Useful formulas and definitions\,:}
In speech recognition and several other pattern recognition applications,
it is useful to associate an ``optimal'' sequence of states to a sequence
of observations, given the parameters of a model. For instance, in the case
of speech recognition, knowing which frames of features ``belong'' to which
state allows to locate the word boundaries across time. This is called the
{\em alignment} of acoustic feature sequences.

A ``reasonable'' optimality criterion consists in choosing the state
sequence (or {\em path}) that brings a maximum likelihood with respect to a
given model.  This sequence can be determined recursively via the {\em
Viterbi algorithm}. This algorithm makes use of two variables\,:
\begin{itemize}
\item the {\em highest} likelihood $\delta_t(i)$ along a {\em single} path
among all the paths ending in state $i$ at time $t$\,:
\[
\delta_t(i) = \max_{q_1,q_2,\cdots,q_{t-1}}
p(q_1,q_2,\cdots,q_{t-1},q^t=q_i,x_1,x_2,\cdots x_t|\Theta)
\]
\item a variable $\psi_t(i)$ which allows to keep track of the ``best
path'' ending in state $i$ at time $t$\,:
\[
\psi_t(i) = \mbox{\hspace{2ex}arg}\hspace{-1ex}\max_{\hspace{-4.5ex}q_1,q_2,\cdots,q_{t-1}}
p(q_1,q_2,\cdots,q_{t-1},q^t=q_i,x_1,x_2,\cdots x_t|\Theta)
\]
\end{itemize}
Note that these variables are vectors of $(N-2)$ elements, $(N-2)$ being
the number of emitting states. With the help of these variables, the
algorithm takes the following steps\,: \\[1em]
\noindent \fbox{\bf The Viterbi Algorithm}
\begin{enumerate}
\item {\bf Initialization}
%
	\begin{eqnarray}
		\delta_1(i) & = & a_{1i} \cdot b_i(x_1), \;\;\;\; 2 \leq i \leq N-1 \nonumber \\
		\psi_1(i) & = & 0 \nonumber
	\end{eqnarray}
%
where, again, $a_{1i}$ are the transitions from the initial non-emitting
state to the emitting states with pdfs $b_{i,\,i = 2 \cdots N-1}(x)$, and
where $b_1(x)$ and $b_N{x}$ do not exist since they correspond to the
non-emitting initial and final states.
%
\item {\bf Recursion}