
% $Id:$


%The output data of each simulation run consists of a number of important
%variables, which are recorded during the simulation. 

%---------------------------------------------------------------------------
\subsection{AutoBayes}
%---------------------------------------------------------------------------
%\AB was then applied to discover more about the variables selected by \TAR.
%One way to explain {\AB} is to view it as an extension of the
%standard EM clustering algorithm.
%When working with unsupervised data,
A well-known method to find structure in large sets of data is to perform
clustering. Clustering is an unsupervised learning method that tries to
estimate class membership matrix and class parameters, only given the
data.
A variety of algorithms can be used, for example, the EM-algorithm 
\cite{dempster-laird-rubin-77-jrssb}.
A number of EM-implementations is available (e.g.,
Autoclass \cite{cheeseman-etal-88-ml,cheeseman-stutz-96-kdd}, EMMIX
\cite{mclachlan-etal-99-jss}, or MCLUST \cite{fraley-raftery-02-tr})
and could be used for this problem.
\comment{ Kevin Murphy's SW etc}

However, in order to refine the statistical model (e.g., by incorporating
other probability distributions for certain variables or to introduce
domain knowledge), the EM-algorithm needs to be modified substantially for
each problem variant, making experimentation a time-consuming and error-prone
undertaking.
Thus, we are using \AB tool to produce customized variants of the EM-algorithm.
\begin{figure*}[!t]
  \centerline{\includegraphics[width=0.65\textwidth]{figures/ab-architecture-latex}}
  \caption{\label{fig:arch}AutoBayes system architecture}
\end{figure*}


 \AB \cite{buntine-fischer-pressburger-99-kdd,fischer-schumann-03-jfp,
gray-etal-02-nips} is a fully automatic program synthesis system
that generates efficient and documented C/C++ code from abstract
statistical model specifications.
Developed at NASA Ames, \AB is 
implemented in approximately 90,000 lines of
documented SWI Prolog code.
From the outside, it looks similar
to a compiler for a very high-level programming language: it takes an 
abstract problem specification in the
form of a (Bayesian) statistical model and translates it into executable code.
On the inside, however, it works quite different: \AB first derives
a customized algorithm skeleton implementing the model and then
transforms it into optimized C/C++ code.
The program synthesis system uses a schema based approach to
translate the statistical problem into smaller problems and then,
after symbolically solving subproblems, to transform the instantiated
customized algorithm into efficient code.

\AB uses \emph{Bayesian Networks} (BNs) \cite{buntine-94-jair,pearl-88-mk}
as a compact internal representation of the statistical models.
Figure~\ref{fig:arch} shows the architecture of \AB. 
The input specification (see below) is first translated into its internal
BN representation. Then, the synthesis kernel applies synthesis schemas from
a domain-specific schema library to assemble the customized algorithm,
heavily supported by an elaborate symbolic subsystem and an efficient
rewriting engine. After optimization, C or C++ code is generated for
various platforms (e.g., embedded systems, Matlab, Simulink, or Octave).
For our experiment, we used AutoBayes to generate code that can be
called from Matlab (a MEX function).

%
%\AB has been developed at NASA Ames and has been used,
%among others, for data analysis of planetary nebulae 
%\comment{\cite{Maxent}}
%and for estimating Mercer kernels on galaxy data
%\comment{\cite{Ashok}}.
%

\begin{figure}
{\footnotesize \addtolength{\baselineskip}{-0.40mm}
\begin{alltt}
\ul{model} mog \ul{as} 'Multivar. Mix of Gaussians';

\ul{const} \ul{int} D \ul{as} 'number of variables' 
\ul{const} \ul{int} N \ul{as} 'number of data points'
      \ul{with} 1 < N;
\ul{const} \ul{int} C  \ul{as} 'number of classes'
      \ul{with} 1 < C; \ul{with} C << N;

\ul{double} phi(1..C) \ul{as} 'class probabilities'
      \ul{with} 1 = sum(_i := 1..C, phi(_i));
\ul{double} mu(1..D, 1..C); 
\ul{double} sigma(1..D, 1..C);

\ul{output} \ul{int} c(1..N) \ul{as} 'latent variable'; 
c(_) ~ discrete(phi);

\ul{data} \ul{double} x(1..D, 1..N); 
x(_i,_j) ~ gauss(mu(_i,c(_j)),sigma(_i,c(_j)));

\ul{max} pr(x | \{phi, mu, sigma\}) 
    \ul{wrt} \{phi, mu, sigma\};
\end{alltt}
\addtolength{\baselineskip}{0.40mm} } \caption{AutoBayes specification for
Gaussian mixture model. Keywords are underlined. \hspace*{\fill} }
\label{tab:mog-spec}
\end{figure}

\begin{figure*}[htb]
\centerline{\includegraphics[width=0.75\textwidth]{figures/entry_2_1_target_norm}}
\caption{Landing spots for each of the 1000 simulation runs. Colors indicate
class membership.
% Members of class~1 (green) only show an extremely small
%$\sigma$ and are represented by an enlarged green dot.
}
\label{fig:target2D}
\end{figure*}

\begin{figure*}[htb]
\centerline{\includegraphics[width=0.75\textwidth]{figures/entry_2_1_inp_par_norm}}
\caption{Clustering results for two simulation input parameters
($p_1$, $p_2$).}
\label{fig:inp_par}
\end{figure*}



The statistical model, which comprises the input to \AB describes
the properties of the data in a fully declarative fashion: for each
problem variable of interest (i.e., observation or parameter),
properties and dependencies are specified via probability
distributions and constraints. Figure~\ref{tab:mog-spec} shows how
our Gaussian mixture model with diagonal covariance
matrices can be represented in \AB's specification language. The
model assumes that the data consists of \texttt{N} points in
\texttt{D} dimensions such that each point belongs to one of
\texttt{C} classes; the first few lines of the specification just
declare these symbolic constants and specify the constraints on
them. 
Each point \texttt{x(1..C,\_j)}
(where \texttt{..} corresponds to Matlab's subrange
operator \texttt{:}, and \texttt{\_i}, \texttt{\_j} are index variables)
is drawn independently from a univariate Gaussian
with mean {\tt mu(\_i,c(\_j))} and standard deviation {\tt
sigma(\_i,c(\_j))}.
%%% \texttt{x(I,J)\~{}gauss(mu(I,c(J)),sigma(I,c(J)))}.
The unknown distribution parameters can be different for each class
and each dimension; hence, we declare them as matrices. The unknown
assignment of the points to the distributions (i.e., classes) is
represented by the latent variable {\tt c}; since we are interested
in the classification results as well (and not only the distribution
parameters), \texttt{c} is declared as \texttt{output}. \texttt{c}
is distributed as a discrete distribution with the relative class
frequencies given by the also unknown vector \texttt{phi}. Since
each point must belong to a class, the sum of the probabilities must
be equal to one.
%Additional constraints (lines~4, 6, 7) express further basic assumptions
%of the model.
Finally, we specify the goal inference task, maximizing the
conditional probability\\
\centerline{\tt pr(x | \{phi, mu, sigma\})}
with
respect to the parameters of interest, {\tt phi}, {\tt mu}, and {\tt
sigma}.
%%% Say something about MLE and MAP
This means that in this case, we are interested in a maximum likelihood 
estimate (MLE) of the model parameters.
More complicated models (e.g., maximum aposteriori estimates, MAP)
can be easily specified by adding priors to the model 
(see Figure~\ref{fig:priors}). 
Note that the model
is completely declarative and does not require the user to prescribe
any algorithmic aspects of the estimation program.

\subsection{\AB: Results}

We have conducted several clustering experiments with data from our
simulation example.
The data produced by the simulation are actually time series data over
a large number of variables (e.g., altitude, G-force, longitude, latitude).
\AB processed these data in order to produce a single vector of data for
each simulation run. Data dimensions obviously include the landing
position, the sum of consumed fuel, maximal structural loads, as
well as a measure of the duration of extended time intervals where the
gravitational forces exceed a safe limit.
With this preprocessing, we obtained a data set with 10 dimensions.
All data were normalized.

These data then were clustered using the Matlab/C code as was generated
by \AB (790 lines of documented C code). We set the maximum number of 
EM-iterations to 30, the termination
error to $1.0\times 10^{-6}$. The generated data analysis algorithm determined
that with 6 classes, a good separation can be achieved%
\footnote{The \AB specification of Fig~\ref{tab:mog-spec} does parameter
estimation for a model with a fixed number of classes. To obtain a the
best number of classes, a simple iteration and comparison with respect to
the model log-likelihood is performed.}.

The results of clustering, projected on the landing spot is shown
in Figure~\ref{fig:target2D}. Different colors indicate into which
class a specific simulation run falls.
Simulation runs in Class~1 (which comprises about 28\% of the total runs) result in a 
landing on target spot with very little deviation.
Thus all runs belonging
to this class are displayed as a single (enlarged) green dot. 
shown as a large green dot, bottom right, of Figure~\ref{fig:target2D}). 
For many of the 
other runs,
however, strong deviations of the landing position from the target exist.

\TAR' speed and succinctness are enabled by taking a coarse-grained view of the data.
(i.e., the output dependent variables are binarized by \BORE and 
the independent variables are discretized via percentile 
divisions).
Figure~\ref{fig:inp_par} 
compares \TAR's  coarse-grained view with the more nuanced view provided by \AB. 
The six
clusters found by \AB are shown 
in Figure~\ref{fig:inp_par} as different colored dots
(green denotes best).
Note that: \bi
\item
the best zone is a green strip that runs top-left to bottom-right 
\item
\TAR's coarse-grain divisions misses the middle of that strip.
\ei




