% file = sato01a.tex

\documentclass[11pt]{article}
\usepackage{jair}
\usepackage{theapa}
\usepackage{graphicx}

\jairheading{15}{2001}{391-454}{6/18}{12/01}
\ShortHeadings{Parameter Learning of Logic Programs
             for Symbolic-statistical Modeling}{Sato \& Kameya}
\firstpageno{391}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{corollary}{Corollary}[section]

\newcommand{\ifff}{\mbox{iff\,}}
\newcommand{\defined}{\stackrel{\rm def}{=}}
\newcommand{\db}{{D\!B}}
\newcommand{\pdb}{P_{D\!B}}
\newcommand{\taudb}{\tau_{D\!B}}
\newcommand{\taudbstar}{\tau_{D\!B}^{\ast}}
\newcommand{\Stil}{\widetilde{S}}
\newcommand{\psidbtilde}{\widetilde{\psi}_{D\!B}}
\newcommand{\ddefined}{\stackrel{{\rm def}}{=}}
\newcommand{\smbf}[1]{{\small\bf {#1}}}
\newcommand{\btheta}{\mbox{\boldmath $\theta$}}
\newcommand{\bpi}{\mbox{\boldmath $\Pi$}}
\newcommand{\bu}{\mbox{\boldmath $u$}}
\newcommand{\ttname}[1]{\mbox{\tt {#1}}}
\newcommand{\pname}[1]{\mbox{\small\em {#1}}}
\newcommand{\ON}{\mbox{$O\!N$}}
\newcommand{\OFF}{\mbox{$O\!F\!F$}}

\newcounter{linenumber}
\newenvironment{listing}{%
	\begin{list}{%
		\footnotesize\arabic{linenumber}:}{%
			\usecounter{linenumber}%
			\setlength{\baselineskip}{13pt}%
			\setlength{\itemsep}{0pt}%
			\setlength{\parsep}{0pt}}}%
	{\end{list}}


\begin{document}

\title{Parameter Learning of Logic Programs
              for Symbolic-statistical Modeling}

\author{\name Taisuke Sato \email sato@mi.cs.titech.ac.jp \\
        \name Yoshitaka Kameya \email kame@mi.cs.titech.ac.jp \\
        \addr Dept. of Computer Science, Graduate School of Information\\
              Science and Engineering, Tokyo Institute of Technology\\
              2-12-1  Ookayama Meguro-ku Tokyo Japan 152-8552
		}

\maketitle

\begin{abstract}
We    propose  a      logical/mathematical framework   for
statistical  parameter learning   of parameterized   logic
programs,   i.e.\     definite clause  programs containing
probabilistic facts with a parameterized distribution.  It
extends the traditional least  Herbrand model semantics in
logic   programming to   {\em   distribution semantics\/},
possible world semantics  with  a probability distribution
which is   unconditionally applicable  to  arbitrary logic
programs   including ones for    HMMs, PCFGs  and Bayesian
networks.

We also propose a new EM  algorithm, the {\em graphical EM
algorithm}, that  runs for a  class of parameterized logic
programs  representing sequential decision processes where
each decision is exclusive  and independent.  It runs on a
new  data   structure  called    {\em  support   graph\/}s
describing the  logical relationship  between observations
and their explanations, and learns parameters by computing
inside  and  outside  probability  generalized  for  logic
programs.

The complexity analysis shows that when combined with OLDT
search for   all   explanations  for  observations,    the
graphical  EM  algorithm, despite  its generality, has the
same time complexity  as existing EM algorithms, i.e.\ the
Baum-Welch algorithm   for   HMMs,  the     Inside-Outside
algorithm  for  PCFGs, and  the  one for  singly connected
Bayesian  networks that  have been developed independently
in each research field.   Learning experiments with  PCFGs
using   two corpora  of   moderate size indicate  that the
graphical  EM algorithm  can  significantly outperform the
Inside-Outside algorithm.
\end{abstract}


\section{Introduction}
\label{sec:introduction}

Parameter learning is common in various fields from neural
networks to reinforcement  learning to statistics.   It is
used to tune  up systems for   their best performance,  be
they classifiers  or   statistical models.   Unlike  these
numerical    systems described   by mathematical  formulas
however, symbolic systems, typically programs, do not seem
amenable  to  any  kind of   parameter learning.  Actually
there has been little  literature on parameter learning of
programs.

This paper is an attempt to incorporate parameter learning
into   computer   programs.    The   reason  is   twofold.
Theoretically we  wish to add the  ability of  learning to
computer programs,   which    the authors believe   is   a
necessary   step   toward building    intelligent systems.
Practically  it   broadens    the  class  of   probability
distributions, beyond   traditionally used numerical ones,
which are available for modeling complex phenomena such as
gene   inheritance, consumer   behavior, natural  language
processing and so on.

The type of  learning  we  consider here is    statistical
parameter learning applied to logic programs.\footnote{
%
In   this  paper,  logic  programs  mean  definite  clause
programs.  A  {\em definite clause program\/}  is a set of
definite clauses.  A {\em definite  clause\/} is a  clause
of   the form  $A\leftarrow\!L_1,\ldots,L_n$ $(0\leq   n)$
where $A,L_1,\ldots,L_n$   are atoms.  $A$  is  called the
head, $L_1,\ldots,L_n$    the  body.   All   variables are
universally quantified.  It  reads  if $L_1$ and  $\cdots$
and  $L_n$ hold, then $A$  holds.   In case of $n=0$,  the
clause is called a   {\em  unit clause}.  A {\em   general
clause\/} is one whose body may contain  negated atoms.  A
program including  general  clauses is  sometimes called a
general program \cite{Lloyd84,Doets94}. }
%
We assume   that facts (unit  clauses)  in  a  program are
probabilistically    true  and    have  a    parameterized
distribution.\footnote{
%
Throughout this paper, for familiarity and readability, we
will  somewhat  loosely use ``distribution''  as a synonym
for ``probability measure''. }
%
Other clauses,  non-unit definite clauses, are always true
as they encode laws  such as ``if  one has a pair of blood
type genes  {\tt a} and {\tt  b}, one's blood type is {\tt
AB}''.  We  call  logic  programs  of  this type   a  {\em
parameterized  logic program\/}   and use for  statistical
modeling in which ground atoms\footnote{
%
In logic programming,  the  adjective ``ground'' means  no
variables contained. }
%
provable from  the program represent our observations such
as ``one's blood type is {\tt AB}''  and the parameters of
the    program are  inferred    by performing ML  (maximum
likelihood) estimation on the observed atoms.

The probabilistic  first-order framework sketched above is
termed {\em statistical abduction} \cite{Sato00a} as it is
an amalgamation  of  statistical  inference  and abduction
where   probabilistic  facts   play  the    role of   {\em
abducible\/}s, i.e.\  primitive hypotheses.\footnote{
%
Abduction  means inference  to the best  explanation for a
set  of observations.   Logically, it  is  formalized as a
search  for an {\em  explanation\/} $E$ such that $E, K\!B
\vdash G$    where  $G$  is   an  atom   representing  our
observation, $K\!B$ a knowledge base and $E$ a conjunction
of atoms chosen from {\em  abducible\/}s, i.e.  a class of
formulas   allowed  as          primitive       hypotheses
\cite{Ka92,Flach00}.   $E$    must  be    consistent  with
$K\!B$. }
%
Statistical abduction is powerful    in that it  not  only
subsumes  diverse symbolic-statistical  frameworks such as
HMMs   (hidden Markov  models,   \citeR{Rabiner89}), PCFGs
(probabilistic         context     free          grammars,
\citeR{Wetherell80,Manning99}) and  (discrete)    Bayesian
networks \cite{Pearl88,Castillo97} but gives us freedom of
using     arbitrarily     complex   logic   programs   for
modeling.\footnote{
%
Existing symbolic-statistical modeling    frameworks  have
restrictions  and  limitations of various  types  compared
with arbitrary logic programs (see  Section~\ref{sec:work}
for details). For example, Bayesian  networks do not allow
recursion.   HMMs  and PCFGs,  stochastic grammars,  allow
recursion  but   lack  variables and    data   structures.
Recursive  logic   programs  are allowed  in \citeS{Ngo97}
framework but they assume  domains are finite and function
symbols seem prohibited.  }

The  semantic basis for   statistical  abduction is   {\em
distribution semantics\/} introduced  by   \citeA{Sato95}.
It defines a    parameterized  distribution, actually    a
probability  measure,   over the   set  of possible  truth
assignments to ground atoms and enables us to derive a new
EM algorithm\footnote{
%
``EM   algorithm'' stands    for   a  class of   iterative
algorithms   for    ML estimation   with  incomplete  data
\cite{McLachlan97}.}
%
for ML estimation called the {\em graphical EM algorithm}
\cite{Kameya00}.

Parameter learning in statistical abduction is done in two
phases,  search  and EM   learning.  Given a parameterized
logic  program and observations,  the first phase searches
for all explanations for  the observations.  Redundancy in
the first  phase    is eliminated by  tabulating   partial
explanations            using       OLDT            search
\cite{Tamaki86,Warren92,Sagonas94,Ramakrishnan95,Shen01}.
It returns a  {\em  support graph\/}  which  is a  compact
representation  of  the  discovered explanations.   In the
second  phase, we  run the graphical   EM algorithm on the
support graph and learn the parameters of the distribution
associated with the  program.   Redundancy in  the  second
phase is removed by the introduction of inside and outside
probability for  logic programs  computed from the support
graph.

The graphical EM algorithm has accomplished, when combined
with OLDT  search   for all explanations,  the   same time
complexity as the specialized  ones, e.g.\  the Baum-Welch
algorithm for HMMs \cite{Rabiner89} and the Inside-Outside
algorithm  for   PCFGs   \cite{Baker79},   despite     its
generality.  What is surprising is that, when we conducted
learning experiments  with PCFGs   using real corpora,  it
outperformed  the  Inside-Outside  algorithm  by orders of
magnitudes in terms of time   for one iteration to  update
parameters.  These   experimental    results enhance   the
prospect  for      symbolic-statistical    modeling     by
parameterized logic programs  of even more complex systems
than stochastic grammars whose modeling has been difficult
simply because of the lack of an appropriate modeling tool
and their sheer complexities.   The contributions of  this
paper therefore are

\begin{itemize}
\item  distribution semantics for parameterized logic programs
       which     unifies  existing    symbolic-statistical
       frameworks,
\item  the graphical EM algorithm (combined with tabulated search),
       a general  yet efficient EM algorithm
       that runs on support graphs and
\item  the prospect suggested by the learning experiments 
       for modeling and learning complex symbolic-statistical
       phenomena.
\end{itemize}

The rest of this  paper  is organized  as follows.   After
preliminaries     in  Section~\ref{sec:preliminaries},   a
probability space   for  parameterized  logic programs  is
constructed   in   Section~\ref{sec:semantics}    as     a
mathematical basis for  the subsequent sections.   We then
propose a  new EM algorithm,   the graphical EM algorithm,
for    parameterized      logic         programs        in
Section~\ref{sec:graphical-EM}.  Complexity analysis    of
the     graphical      EM algorithm     is   presented  in
Section~\ref{sec:complexity} for HMMs, PCFGs, pseudo PCSGs
and sc-BNs.\footnote{
%
Pseudo   PCSGs (probabilistic context sensitive  grammars)
are a context-sensitive extension of PCFGs proposed by
\citeA{Charniak94}.   sc-BN  is a  shorthand  for a singly
connected Bayesian network \cite{Pearl88}.  }
%
Section~\ref{sec:experiments}    contains     experimental
results of parameter  learning with PCFGs by the graphical
EM   algorithm  using real   corpora that  demonstrate the
efficiency  of the    graphical  EM algorithm.  We   state
related work  in Section~\ref{sec:work}, followed       by
conclusion in Section~\ref{sec:conclusion}.  The reader is
assumed   to be   familiar  with  the    basics of   logic
programming  \cite{Lloyd84,Doets94},    probability theory
\cite{Chow97}, Bayesian networks \cite{Pearl88,Castillo97}
and stochastic grammars \cite{Rabiner89,Manning99}.


\section{Preliminaries}
\label{sec:preliminaries}

Since our subject   intersects  logic programming and   EM
learning which are quite different  in nature, we separate
preliminaries.

\subsection{Logic Programming and OLDT}

In logic programming, a program $\db$ is a set of definite
clauses\footnote{
%
We do not deal with general  logic programs in this paper.
}
%
and the execution is  search for an {\em SLD refutation\/}
of a given goal $\leftarrow\!G$.  The top-down interpreter
recursively selects  the    next goal   and   unfolds   it
\cite{Tamaki84} into subgoals using a nondeterministically
chosen clause.  The computed result by the SLD refutation,
i.e.\   a  {\em  solution},  is  an   answer  substitution
(variable  binding)   $\theta$    such  that   $\db\vdash$
$G\theta$.\footnote{
%
By a solution we    ambiguously  mean   both   the  answer
substitution $\theta$    itself  and   the    proved  atom
$G\theta$, as one gives the other.  }
%
Usually   there    is more    than   one   refutation  for
$\leftarrow\!G$, and  the search space for all refutations
is described by an {\em SLD tree\/}  which may be infinite
depending   on    the         program and   the       goal
\cite{Lloyd84,Doets94}.

More often than  not, applications require  all solutions.
In natural language processing for instance, a parser must
be able  to find all  possible  parse  trees for  a  given
sentence as  every one of  them is  syntactically correct.
Similarly in statistical abduction, we need to examine all
explanations to    determine the most    likely  one.  All
solutions are obtained  by searching the entire SLD  tree,
and there is a choice of the  search strategy.  In Prolog,
the standard  logic  programming language, backtracking is
used  to  search for all solutions  in  conjunction with a
fixed     search   order  for    goals   (textually   from
left-to-right)  and  clauses (textually top-to-bottom) due
to the ease and simplicity of implementation.

The  problem    with  backtracking  is   that  it  forgets
everything until   up  to the  previous  choice point, and
hence it is quite likely to prove the  same goal again and
again, resulting  in exponential  search time.  One answer
to avoid this problem   is to store computed  results  and
reuse them whenever  necessary.  OLDT is  such an instance
of                  memoizing                       scheme
\cite{Tamaki86,Warren92,Sagonas94,Ramakrishnan95,Shen01}.
Reuse of proved  subgoals in OLDT search often drastically
reduces search  time  for  all solutions, especially  when
refutations   of    the top   goal  include   many  common
sub-refutations.   Take as   an example  a  logic  program
coding  an  HMM.   For  a given  string   $s$, there exist
exponentially many transition paths that output $s$.  OLDT
search applied to the    program however only takes   time
linear in  the length of  $s$ to find  all of  them unlike
exponential  time by Prolog's  backtracking search.

What does OLDT  have   to do with statistical   abduction?
From  the   viewpoint of statistical  abduction,  reuse of
proved  subgoals, or   equivalently, structure  sharing of
sub-refutations    for  the   top-goal $G$  brings   about
structure sharing of explanations for  $G$, in addition to
the  reduction of   search  time mentioned above,  thereby
producing  a  highly    compact   representation   of  all
explanations for $G$.

\subsection{EM Learning}
\label{sub:EM}

Parameterized  distributions    such  as   the multinomial
distribution    and  the    normal  distribution   provide
convenient modeling   devices  in  statistics.   Suppose a
random  sample  $x_1,\ldots,x_T$ of size   $T$ on a random
variable $X$ drawn from a distribution $P(X=x\mid \theta)$
parameterized by unknown $\theta$, is observed.  The value
of  $\theta$  is determined  by  ML  estimation as the MLE
(maximum  likelihood estimate)  of  $\theta$, i.e.\ as the
maximizer  of  the   likelihood   $\prod_{1\leq  i\leq  T}
P(x_i\mid \theta)$.

Things get  much  more  difficult   when  data  are   {\em
incomplete}.    Think of  a     probabilistic relationship
between non-observable cause $X$ and observable effect $Y$
such as one between  diseases and symptoms in medicine and
assume that $Y$ does not uniquely determine the cause $X$.
Then  $Y$  is incomplete in  the  sense that  $Y$ does not
carry enough information to completely determine $X$.  Let
$P(X=x,Y=y\mid  \theta)$    be   a  parameterized    joint
distribution over $X$ and $Y$.   Our task is to perform ML
estimation on $\theta$   under the condition that  $X$  is
non-observable   while   $Y$     is  observable.       Let
$y_1,\ldots,y_T$ be a random sample of size $T$ drawn from
the  marginal  distribution $P(Y=y\mid   \theta)  = \sum_x
P(X=x,Y=y\mid \theta)$.   The MLE  of $\theta$ is obtained
by  maximizing  the   likelihood   $\prod_{1\leq i\leq  T}
P(y_i\mid \theta)$ as a function of $\theta$.

While mathematical formulation looks  alike in both cases,
the  latter, ML estimation  with  incomplete data, is  far
more complicated  and direct maximization  is  practically
impossible  in many  cases.   People  therefore  looked to
indirect approaches to tackle the problem of ML estimation
with incomplete data to which  the EM algorithm has been a
standard solution \cite{Dempster77,McLachlan97}.  It is an
iterative algorithm   applicable  to  a   wide   class  of
parameterized   distributions  including the   multinomial
distribution and the normal distribution such that the MLE
computation is replaced by  the  iteration of two  easier,
more  tractable    steps. At $n$-th   iteration,  it first
calculates   the value  of $Q$   function introduced below
using        current     parameter   value  $\theta^{(n)}$
(E-step)\footnote{
%
$Q$ function is  related to ML  estimation as follows.  We
assume  here   only one  data,    $y$, is observed.   From
Jensen's  inequality  \cite{Chow97} and the  concavity  of
$\ln$ function, it follows that
%
\[
   \sum_x P(x\mid y,\theta^{(n)})  \ln P(x\mid y,\theta)
    - \sum_x  P(x\mid y,\theta^{(n)})\ln P(x\mid y,\theta^{(n)})
   \leq 0
\]
and hence that
% 
\begin{eqnarray*}
   \lefteqn{
     Q(\theta \mid \theta^{(n)}) - Q(\theta^{(n)} \mid \theta^{(n)})
   } \\
   & = &
     \sum_x P(x\mid y, \theta^{(n)}) \ln P(x\mid y,\theta)
       - \sum_x P(x\mid y, \theta^{(n)})\ln P(x\mid y,\theta^{(n)})
       + \ln P(y\mid \theta) - \ln P(y\mid \theta^{(n)}) \\
   & \leq &
       \ln P(y\mid  \theta) - \ln P(y\mid  \theta^{(n)}).
\end{eqnarray*}
%
Consequently, we have
\[ Q(\theta \mid \theta^{(n)}) \geq Q(\theta^{(n)} \mid \theta^{(n)})
   \;\Rightarrow\;
   \ln p(y\mid  \theta) \geq \ln p(y\mid  \theta^{(n)})
   \;\Rightarrow\;
   p(y\mid  \theta) \geq p(y\mid  \theta^{(n)}).
\]
}
%
:
%
\begin{equation}
Q(\theta \mid \theta^{(n)})\ddefined
   \sum_x P(x\mid y,\theta^{(n)}) \ln P(x, y \mid \theta).
\label{eq:Q1}
\end{equation}
%
Next,  it maximizes $Q(\theta  \mid  \theta^{(n)})$ as a
function of $\theta$ and updates $\theta^{(n)}$ (M-step):
%
\begin{equation}
\theta^{(n+1)} = \mbox{argmax}_{\theta}Q(\theta \mid \theta^{(n)}).
\label{eq:Q2}
\end{equation}
%
Since the old value  $\theta^{(n)}$ and the  updated value
$\theta^{(n+1)}$ do not  necessarily coincide, the E-steps
and M-steps are  iterated until  convergence, during which
the     (log)  likelihood   is    assured    to   increase
monotonically \cite{McLachlan97}.

Although   the   EM    algorithm   merely   performs local
maximization, it is used  in a variety  of settings due to
its simplicity and relatively  good performance.  One must
notice however that the EM algorithm is just a class name,
taking  different   form depending  on   distributions and
applications.  The  development of a concrete EM algorithm
such as the Baum-Welch algorithm for HMMs \cite{Rabiner89}
and  the Inside-Outside algorithm for PCFGs \cite{Baker79}
requires individual effort for each case.


\section{Distribution Semantics}
\label{sec:semantics}

In this  section, we  introduce {\em  parameterized  logic
programs\/} and  define their declarative  semantics.  The
basic  idea is as  follows.  We  start with  a set  $F$ of
probabilistic   facts (atoms) and   a  set $R$ of non-unit
definite clauses.  Sampling from $F$ determines a set $F'$
of true atoms, and the least Herbrand model of $F' \cup R$
determines the truth value of {\em every\/} atom in $\db =
F \cup R$.  Hence every atom can be considered as a random
variable,  taking on $1$ (true)  or  $0$ (false).  In what
follows,  we  formalize   this process   and construct the
underlying probability space for the denotation of $\db$.

\subsection{Basic Distribution $P_F$}
\label{subsec:P_F}

Let  $\db = F\cup   R$ be a definite  clause  program in a
first-order  language   $\cal   L$  with   countably  many
variables, function symbols  and  predicate symbols  where
$F$ is a set  of unit clauses ({\em facts})  and $R$ a set
of non-unit clauses  ({\em rules}).  In the sequel, unless
otherwise stated, we consider for  simplicity $\db$ as the
set of all  ground instances of  the clauses in $\db$, and
assume that  $F$  and $R$  consist  of  countably infinite
ground clauses (the finite case is similarly treated).  We
then construct a probability space for $\db$ in two steps.
First we introduce   a  probability space over  the   {\em
Herbrand interpretation\/}s\footnote{
%
A  {\em  Herbrand interpretation\/} interprets  a function
symbol uniquely as a  function on ground terms and assigns
truth values to ground atoms.  Since the interpretation of
function    symbols    is   common    to   all    Herbrand
interpretations, given ${\cal  L}$, they have a one-to-one
correspondence with truth assignments  to ground  atoms in
${\cal L}$.  So we do not distinguish them.  }
%
of $F$ i.e.\ the truth assignments to ground atoms in $F$.
Next we extend it to a probability space over the Herbrand
interpretations of {\em  all\/} ground atoms in ${\cal L}$
by using the least model semantics \cite{Lloyd84,Doets94}.

Let $A_1,A_2,\ldots$ be   a fixed enumeration  of atoms in
$F$.   We regard  an  infinite  vector  $\omega =  \langle
x_1,x_2,\ldots\rangle$  of $0$s  and   $1$s as a  Herbrand
interpretation   of $F$  in   such  a  way   that for   $i
=1,2,\ldots$ $A_i$ is  true (resp.\ false)  if and only if
$x_i = 1$ (resp.\ $x_i = 0$).  Under this isomorphism, the
set  of  all  possible   Herbrand  interpretations  of $F$
coincides with the Cartesian product:
%
\[
\Omega_F      \ddefined \prod_{i=1}^{\infty}{\{0,1\}}_i.
\]
%

We construct  a probability measure  $P_F$ over the sample
space $\Omega_F$\footnote{
%
We   regard  $\Omega_F$ as a  topological   space with the
product topology such that each $\{0,1\}$ is equipped with
the discrete topology.}
%
from   a  collection    of    finite  joint  distributions
$P^{(n)}_F(A_1=x_1,\ldots,A_n=x_n)$  $(n=1,2,\ldots,   x_i
\in \{0, 1\}, 1\leq i\leq n)$ such that
%
\begin{equation}
\left\{ \hspace{1em}
\begin{array}{l}
  0\leq P^{(n)}_F(A_1=x_1,\ldots,A_n=x_n)  \leq  1 \\
  \sum_{x_1,\ldots,x_n}P^{(n)}_F(A_1=x_1,\ldots,A_n=x_n)  =  1 \\
  \sum_{x_{n+1}}P^{(n+1)}_F(A_1=x_1,\ldots,A_{n+1}=x_{n+1})
                =  P^{(n)}_F(A_1=x_1,\ldots,A_n=x_n).
\end{array}
\right.
\label{eq:P_F}
\end{equation}
%
The  last   equation   is called  the {\em   compatibility
condition}.   It can   be  proved \cite{Chow97}  from  the
compatibility  condition  that there exists  a probability
space $(\Omega_F, {\cal F},     P_F)$ where $P_F$   is   a
probability  measure on  ${\cal F}$,  the minimal $\sigma$
algebra containing open sets  of $\Omega_F$, such that for
any $n$,
%
\[
P_F(A_1=x_1,\ldots,A_n=x_n) = P^{(n)}_F(A_1=x_1,\ldots,A_n=x_n).
\]
%
We call  $P_F$ a {\em  basic  distribution}.\footnote{
%
This  naming of  $P_F$,  despite its  being a  probability
measure, partly reflects  the observation  that it behaves
like    an        infinite         joint      distribution
$P_F(A_1=x_1,A_2=x_2,\ldots)$  for  an   infinite   random
vector     $\langle  A_1,A_2,\ldots  \rangle$  of    which
$P^{(n)}_F(A_1=x_1,\ldots,A_n=x_n)$ $(n=1,2,\ldots)$   are
marginal distributions.   Another reason is intuitiveness.
These considerations apply to  $\pdb$ defined in the  next
subsection as well.}

The choice  of $P_F^{(n)}$   is   free  as long   as   the
compatibility  condition    is met.    If   we  want   all
interpretations   to  be    equiprobable,  we should   set
$P_F^{(n)}(A_1=x_1,\ldots,A_n=x_n)   = 1/{2^n}$ for  every
$\langle  x_1,\ldots,x_n \rangle$.  The resulting $P_F$ is
a uniform distribution over  $\Omega_F$ just like the  one
over the unit interval $[0,1]$.  If, on the other hand, we
stipulate no  interpretation  except  $\omega_0  = \langle
c_1,c_2,\ldots \rangle$   should be possible,  we put, for
each $n$,
%
\[
P_F^{(n)}(A_1=x_1,\ldots,A_n=x_n)  = 
\left\{
   \begin{array}{lll}
     1    & \mbox{if} & \forall i\; x_i = c_i\;  (1\leq i\leq n)  \\
     0    & \mbox{o.w.}
   \end{array}
\right.
\]
%
Then  $P_F$ places all probability  mass on $\omega_0$ and
gives probability $0$ to the rest.\\

Define a {\em parameterized logic program\/} as a definite
clause program\footnote{
%
Here clauses are not necessarily ground.
}
%
$\db = F\cup R$ where $F$ is a set of unit clauses, $R$ is
a set of non-unit clauses such that no  clause head in $R$
is unifiable with a unit clause in $F$ and a parameterized
basic   distribution $P_F$  is  associated  with   $F$.  A
parameterized  $P_F$  is  obtained  from a   collection of
parameterized   joint      distributions  satisfying   the
compatibility   condition.   Generally,  the more  complex
$P_F^{(n)}$'s are, the more flexible  $P_F$ is, but at the
cost of tractability.  The  choice of parameterized finite
distributions made by \citeA{Sato95} was simple:
%
\begin{eqnarray}
\lefteqn{P_F^{(2n)}(\ON_1=x_1,\OFF_2=x_2,\ldots,\OFF_{2n}=x_{2n}
                 \mid  \theta_1,\ldots, \theta_n )} \nonumber \\
 & = &
 \prod_{i=1}^{n}
     P_{bs}(\ON_{2i-1}=x_{2i-1}, \OFF_{2i}=x_{2i}\mid \theta_i) \nonumber \\
\mbox{where} \nonumber \\
\lefteqn{P_{bs}(\ON_{2i-1}=x_{2i-1}, \OFF_{2i}=x_{2i} \mid \theta_i)}
                   \nonumber \\
 & = & \left\{
            \begin{array}{cl}
            0          & \mbox{if}\; x_{2i-1} = x_{2i} \\
            \theta_i     & \mbox{if}\; x_{2i-1} = 1,\; x_{2i} = 0 \\
            1 - \theta_i & \mbox{if}\; x_{2i-1} = 0,\; x_{2i} = 1.
            \end{array}
       \right.  \label{dist:bs}
\end{eqnarray}
%
$P_{bs}(\ON_{2i-1}=x_{2i-1},\OFF_{2i}    =      x_{2i}\mid
\theta_i)$ $(1\leq i\leq   n)$ represents a  probabilistic
binary switch,  i.e.\   a   Bernoulli trial,  using    two
exclusive atoms $\ON_{2i-1}$ and $\OFF_{2i}$ in such a way
that either one  of them is true  on each  trial but never
both.    $\theta_i$    is  a  parameter    specifying  the
probability   that the switch   $i$ is  on.  The resulting
$P_F$ is a probability measure  over the infinite  product
of independent binary outcomes.   It might look too simple
but expressive enough for Bayesian networks, Markov chains
and HMMs \cite{Sato95,Sato97}.


\subsection{Extending $P_F$ to $\pdb$ }
\label{T_P}

In   this  subsection, we  extend   $P_F$ to a probability
measure $\pdb$ over the {\em possible world\/}s for ${\cal
L}$, i.e.\ the  set of all  possible truth  assignments to
ground  atoms in  ${\cal L}$  through   the least Herbrand
model \cite{Lloyd84,Doets94}.  Before  proceeding however,
we need a  couple of notations.  For  an atom  $A$, define
$A^x$ by
%
\[
\left\{
\begin{array}{lll}
A^x = A      & \mbox{if} & x=1 \\
A^x = \neg A & \mbox{if} & x=0.
\end{array}
\right.
\]
%
Next take a  Herbrand interpretation $\nu \in \Omega_F$ of
$F$.  It  makes some atoms in  $F$ true and  others false.
Let   $F_{\nu}$ be the set of    atoms made true by $\nu$.
Then  imagine  a definite clause  program  ${\db}' = R\cup
F_{\nu}$  and   its   least   Herbrand model    $M_{\db'}$
\cite{Lloyd84,Doets94}.  $M_{\db'}$   is  characterized as
the least fixed point of a mapping $T_{\db'}(\cdot)$ below
%
\[
 T_{\db'}(I) \ddefined
    \left\{ A \left|
               \begin{array}{l}
                  \mbox{there is some}\;
                  A\leftarrow B_1,\ldots,B_k \in  \db'  \;(0\leq k) \\
                  \mbox{such that}\; \{ B_1,\ldots,B_k\} \subseteq I \\
               \end{array}
             \right.
    \right\}
\]
%
where   $I$  is a  set  of  ground  atoms.\footnote{
%
$I$ defines, mutually, a Herbrand interpretation such that
a ground atom $A$ is true if and only if $A \in I$. A {\em
Herbrand  model\/}  of    a    program  is    a   Herbrand
interpretation that makes every  ground instance of  every
clause in the program true. }
%
Or equivalently, it is inductively defined by
%
\begin{eqnarray*}
I_0     & = & \emptyset \\
I_{n+1} & = &  T_{\db'}(I_n)\\
M_{\db'} & = & \bigcup_n I_n.
\end{eqnarray*}
%
Taking into account the fact that $M_{\db'}$ is a function
of $\nu\in \Omega_F$, we   henceforth employ a  functional
notation $M_{\db}(\nu)$ to denote $M_{\db'}$.

Turning  back,    let    $A_1,A_2,\ldots$  be   again   an
enumeration, but of   {\em all\/} ground   atoms in ${\cal
L}$.\footnote{
%
Note that this  enumeration enumerates ground atoms in $F$
as well.}
%
Form $\Omega_{\db}$,   similarly to  $\Omega_{F}$, as  the
Cartesian  product  of  denumerably many   $\{0,1\}$'s and
identify it    with the  set   of  all  possible  Herbrand
interpretations  of the  ground atoms  $A_1,A_2,\ldots$ in
${\cal  L}$, i.e.\ the {\em  possible world\/}s for ${\cal
L}$.  Then  extend $P_F$ to  a probability  measure $\pdb$
over  $\Omega_{\db}$  as follows.   Introduce  a series of
finite      joint      distributions        $P^{(n)}_{\db}
(A_1=x_1,\ldots,A_n=x_n)$ for $n=1,2,\ldots$ by
%
\[
\begin{array}{rcl}
[A_1^{x_1} \wedge\cdots\wedge A_n^{x_n}]_F &
   \ddefined &
      \{\nu \in \Omega_F \mid M_{\db}(\nu)\models
                A_1^{x_1}\wedge\cdots\wedge A_n^{x_n}\}\\
P^{(n)}_{\db}(A_1=x_1,\ldots,A_n=x_n) &
   \ddefined &
       P_F([A_1^{x_1} \wedge\cdots\wedge A_n^{x_n}]_F).
\end{array}
\]
%
Note   that  the  set  $[A_1^{x_1}  \wedge  \cdots  \wedge
A_n^{x_n}]_F$ is  $P_F$-measurable and   by    definition,
$P^{(n)}_{\db}$'s satisfy the compatibility condition
%
\[
\sum_{x_{n+1}} P^{(n+1)}_{\db}(A_1=x_1,\ldots,A_{n+1}=x_{n+1})
   = P^{(n)}_{\db}(A_1=x_1,\ldots,A_n=x_n).
\]
%
Hence   there exists a    probability measure  $\pdb$ over
$\Omega_{\db}$ which is an extension of $P_F$ such that
%
\[
\pdb(A_1=x_1,\ldots,A_{n}=x_n) = P_{F}(A_1=x_1,\ldots,A_n=x_n)
\]
%
for any finite atoms $A_1,\ldots,A_n$ in $F$ and for every
binary   vector  $\langle x_1,\ldots,x_n\rangle$ $(x_i \in
\{0,1\},\; 1\leq i\leq n)$.  Define the {\em denotation\/}
of  the  program $\db   = F  \cup R$   w.r.t.\ $P_F$ to be
$\pdb$.  The denotational semantics of parameterized logic
programs   defined above    is called  {\em   distribution
semantics}.   As remarked  before,  we regard  $\pdb$ as a
kind    of     infinite          joint        distribution
$\pdb(A_1=x_1,A_2=x_2,\ldots)$.  Mathematical   properties
of $\pdb$  are listed in  Appendix~\ref{appendix:A}  where
our semantics is proved to be an extension of the standard
least model  semantics  in logic  programming  to possible
world semantics with a probability measure.

\subsection{Programs as Distributions}

Distribution  semantics views parameterized logic programs
as  expressing distributions.  Traditionally distributions
have been expressed by using mathematical formulas but the
use of  programs as (discrete)  distributions gives us far
more freedom and flexibility than mathematical formulas in
the construction of  distributions because  they have {\em
recursion\/} and  arbitrary  composition.  In particular a
program  can  contain infinitely many random  variables as
probabilistic  atoms   through   recursion, and hence  can
describe  stochastic   processes that  potentially involve
infinitely many random variables such as Markov chains and
derivations in PCFGs \cite{Manning99}.\footnote{
% 
An infinite derivation can occur  in PCFGs.  Take a simple
PCFG $\{ p:S  \rightarrow a, q:S  \rightarrow SS \}$ where
$S$ is a start symbol, $a$ a terminal symbol, $p+q =1$ and
$p,q>0$.  In this PCFG,  $S$  is rewritten either  to  $a$
with probability $p$ or to $SS$ with probability $q$.  The
probability of the occurrence of an infinite derivation is
calculated   as   $\mbox{max}\;\{0,   1-(p/q)\}$  which is
non-zero when $q>p$ \cite{Chi98}.  }

Programs  also enable us   to {\em procedurally\/} express
complicated constraints on distributions such as ``the sum
of  occurrences  of  alphabets {\tt  a} or  {\tt  b} in an
output string of an  HMM must  be  a multiple  of three''.
This feature, procedural expression of arbitrarily complex
(discrete)   distributions,      seems quite   helpful  in
symbolic-statistical modeling.

Finally,   providing mathematically  sound semantics   for
parameterized   logic    programs   is   one   thing,  and
implementing distribution  semantics in a tractable way is
another.   In the next  section, we investigate conditions
on  parameterized  logic programs  which  make probability
computation tractable, thereby  making  them usable as   a
means for large scale symbolic-statistical modeling.


\section{Graphical EM Algorithm}
\label{sec:graphical-EM}

According to the  preceding section, a parameterized logic
program $\db = F\cup R$  in a first-order language  ${\cal
L}$    with    a   parameterized     basic    distribution
$P_F(\cdot\!\mid\!\theta)$       over       the   Herbrand
interpretations    of  ground  atoms  in   $F$ specifies a
parameterized distribution  $\pdb(\cdot \mid \theta)$ over
the  Herbrand  interpretations for  ${\cal  L}$.  In  this
section,   we develop,  step   by step,   an efficient  EM
algorithm for   the parameter  learning  of  parameterized
logic programs  by interpreting  $\pdb$ as  a distribution
over the observable and non-observable events.  The new EM
algorithm is termed the {\em  graphical EM algorithm}.  It
is  applicable   to arbitrary   logic programs  satisfying
certain   conditions  described  later provided  the basic
distribution is   a direct  product  of   multi-ary random
switches,  which is  a slight  complication of the  binary
ones introduced in Section~\ref{subsec:P_F}.

From this  section on,  we  assume that $\db$  consists of
usual definite clauses containing (universally quantified)
variables.   Definitions and   changes  relating   to this
assumption are listed   below.   For a predicate   $p$, we
introduce $\ifff(p)$, the {\em iff definition of $p$\/} by
%
\begin{eqnarray*}
\ifff(p)
 & \defined &
   \forall x\;
     \left(\rule{0em}{0.8em} p(x)\leftrightarrow
       \exists y_1 (x=t_1 \wedge W_1) \vee \cdots \vee
       \exists y_n (x=t_n \wedge W_n) \right).
\end{eqnarray*}
%
Here $x$ is  a vector of new  variables of length equal to
the arity of $p$,  $p(t_i)\leftarrow W_i$ $(1\leq i\leq n,
0\leq n)$, an  enumeration of clauses  about $p$ in $\db$,
and  $y_i$,   a   vector  of    variables  occurring    in
$p(t_i)\leftarrow W_i$. Then define $comp(R)$ as follows.
%
\begin{eqnarray*}
head(R)
 & \defined &
     \{ B \mid B\; \mbox{is a ground instance of a clause head
                   appearing in}\; R \} \\
\ifff(R)
 & \defined &
    \{ \ifff(p) \mid p\;\mbox{appears in a clause head in}\; R \} \\
E_q
 & \defined &
        \{f(x)=f(y)\rightarrow x=y
                   \mid f\;\mbox{is a function symbol}\} \\
 &                   &
        \hspace{1em} \cup\; \{f(x)\neq g(y)
                   \mid \mbox{$f$ and $g$ are
                           different function symbols}\} \\
 &                   &
        \hspace{1em} \cup\; \{t\neq x
                    \mid t \;\mbox{is a term properly containing}\;x\}\\
comp(R)
 & \defined &
         \ifff(R)\;\cup\;E_q
\end{eqnarray*}
%
$E_q$, Clark's  equational    theory       \cite{Clark78},
deductively simulates unification. Likewise $comp(R)$ is a
first-order   theory   which    deductively  simulates SLD
refutation with  the help of  $E_q$ by  replacing a clause
head atom with the clause body \cite{Lloyd84,Doets94}.

We  here    introduce  some  definitions   which  will  be
frequently used.  Let $B$ be an atom.  An {\em explanation
for $B$\/} w.r.t.\  $\db = F\cup  R$ is a conjunction  $S$
such that $S,  R \vdash B$, and  as a set comprised of its
conjuncts, $S \subset F$ holds and no proper subset of $S$
satisfies this.  The set   of all explanations for $B$  is
called  the {\em support set for  $B$\/} and designated by
$\psi_{\db}(B)$.\footnote{
%
This definition of a support set differs from the one used
by \citeA{Sato95} and \citeA{Kameya00}.  }


\subsection{Motivating Example}
\label{subsec:DB-b}

First   of all,  we review   distribution  semantics by  a
concrete example.  Consider the following program $\db_b =
F_b   \cup  R_b$ in  Figure~\ref{program:ABO} modeling how
one's blood type     is determined by blood    type  genes
probabilistically inherited from the parents.\footnote{
%
When we implicitly   emphasize the  procedural  reading of
logic  programs,  Prolog  conventions       are   employed
\cite{Sterling86}. Thus, {\tt ;} stands for ``or'', {\tt ,}
``and''   {\tt :-} ``implied   by'' respectively.  Strings
beginning  with   a  capital    letter   are  (universally
quantified) variables, but  quoted ones such  as {\tt 'A'}
are  constants.  The underscore   {\tt \_} is an anonymous
variable.  }
%
\begin{figure}[htbp]
\[
\begin{array}{lcl}
R_b & = &
\left\{
\mbox{\tt
  \begin{tabular}{rcl}
  btype('A')	& :- 
	& (gtype(a,a) ; gtype(a,o) ; gtype(o,a)). \\
  btype('B')	& :- 
	& (gtype(b,b) ; gtype(b,o) ; gtype(o,b)). \\
  btype('O')	& :- 
	& gtype(o,o). \\
  btype('AB')	& :- 
	& (gtype(a,b) ; gtype(b,a)). \\
  gtype(X,Y)	& :- 
	& gene(father,X), gene(mother,Y).\\
  gene(P,G)		& :- 
	& msw(gene,P,G).\\
  \end{tabular}
}
\right.\\
\\
F_b & = &
   \{ \mbox{\tt msw(gene,father,a)}, \mbox{\tt msw(gene,father,b)},
      \mbox{\tt msw(gene,father,o)}, \\
	&	& \hspace*{2em}
     \mbox{\tt msw(gene,mother,a)}, \mbox{\tt msw(gene,mother,b)},
     \mbox{\tt msw(gene,mother,o)} \}

\end{array}
\]
\caption{ABO blood type program $\db_b$}
\label{program:ABO}
\end{figure}

The first four  clauses in  $R_b$  state a  blood type  is
determined by a genotype, i.e.\ a pair of blood type genes
{\tt   a}, {\tt  b}  and   {\tt o}.    For instance,  {\tt
btype('A'):-} {\tt (gtype(a,a) ; gtype(a,o) ; gtype(o,a))}
says  that  one's blood  type   is  {\tt A}  if his  (her)
genotype  is $\langle  {\tt  a},{\tt  a}\rangle$, $\langle
{\tt    a},{\tt  o}\rangle $    or  $\langle  {\tt o},{\tt
a}\rangle $.  These are propositional rules.

Succeeding clauses state general rules in terms of logical
variables.  The fifth  clause says that regardless of  the
values of {\tt  X}  and {\tt  Y}, event {\tt   gtype(X,Y)}
(one's having genotype $\langle {\tt X}, {\tt Y} \rangle$)
is caused by two  events, {\tt gene(father,X)} (inheriting
gene  {\tt X}  from  the father)  and {\tt gene(mother,Y)}
(inheriting  gene   {\tt  Y}   from  the   mother).   {\tt
gene(P,G):-}  {\tt msw(gene,P,G)}  is  a clause connecting
rules in $R_b$   with probabilistic facts   in $F_b$.   It
tells us that the gene {\tt G} is  inherited from a parent
{\tt    P}      if  a    choice     represented   by  {\tt
msw(gene,P,G)}\footnote{
%
{\tt  msw}   is   an  abbreviation of  ``multi-ary  random
switch'' and   ${\tt  msw}(\cdot,\cdot,\cdot)$ expresses a
probabilistic  choice from  finite  alternatives.   In the
framework of statistical  abduction, {\tt  msw} atoms  are
abducibles  from which  explanations  are constructed as a
conjunction.}
%
is made.  The genetic knowledge that the choice of {\tt G}
is by chance and made from $\{{\tt a}, {\tt b}, {\tt o}\}$
is expressed  by specifying a  joint distribution $F_b$ as
follows.
%
\[
P_{F_b}(\mbox{\tt msw(gene,{\it t},a)}=x,
		\mbox{\tt msw(gene,{\it t},b)}=y,
        \mbox{\tt msw(gene,{\it t},o)}=z 
       \mid \theta_a, \theta_b, \theta_o)
   \defined \theta^x_a \theta^y_b \theta^z_o
\]
%
where $x,y,z  \in   \{0,1\}$,  $x+y+z  =   1$,  $\theta_a,
\theta_b, \theta_o  \in [0,1]$,  $\theta_a  + \theta_b   +
\theta_o =1$  and {\it t}  is either {\tt father}  or {\tt
mother}.  Thus $\theta_a$ is the probability of inheriting
gene  {\tt a} from a  parent.  Statistical independence of
the  choice of gene, once from  {\tt father} and once from
{\tt mother}, is expressed by putting
%
\begin{eqnarray*}
\lefteqn{
   P_{F_b}(\,	\mbox{\tt msw(gene,father,a)}=x,
				\mbox{\tt msw(gene,father,b)}=y,
                \mbox{\tt msw(gene,father,o)}=z,
  		} \\
 & & \hspace{1em}
			\mbox{\tt msw(gene,mother,a)}=x',
			\mbox{\tt msw(gene,mother,b)}=y',
            \mbox{\tt msw(gene,mother,o)}=z' \\
 & & \hspace{1em}
            \mid \theta_a, \theta_b, \theta_o) \\
 & = &
    P_{F_b}(x,y,z    \mid \theta_a, \theta_b, \theta_o)
    P_{F_b}(x',y',z' \mid \theta_a, \theta_b, \theta_o).
\end{eqnarray*}

In this  setting, atoms  representing our  observation are
$obs(\db_b)=$   $\{\mbox{\tt  btype('A')},$     $\mbox{\tt
btype('B')},$   $\mbox{\tt    btype('O')},$     $\mbox{\tt
btype('AB')}\}$.   We observe  one of  them, say \mbox{\tt
btype('A')}, and infer a possible  explanation $S$, i.e.\ a
minimal   conjunction    of      abducibles     $\mbox{\tt
msw(gene,$\cdot$,$\cdot$)}$ such that
%
\[ S, R_b  \vdash \mbox{{\tt btype('A')}.} \]
%
$S$  is  obtained  by applying  a   special SLD refutation
procedure   to     the  goal    $\leftarrow\!    \mbox{\tt
btype('A')}$ which preserves {\tt msw} atoms resolved upon
in the refutation.  Three explanations are found.
%
\begin{eqnarray*}
S_1 & = &
	\mbox{\tt msw(gene,father,a)}\wedge
	\mbox{\tt  msw(gene,mother,a)} \\
S_2 & = &
	\mbox{\tt msw(gene,father,a)}\wedge
	\mbox{\tt  msw(gene,mother,o)} \\
S_3 & = &
	\mbox{\tt msw(gene,father,o)}\wedge
	\mbox{\tt  msw(gene,mother,a)}
\end{eqnarray*}
%
So $\psi_{{\db}_b}(\mbox{\tt  btype(a)})$, the support set
for    \mbox{\tt   btype(a)},  is $\{S_1,S_2,S_3\}$.   The
probability  of each explanation  is respectively computed
as  $P_{F_b}(S_1)  =    \theta^2_a$  and  $P_{F_b}(S_2)  =
P_{F_b}(S_3)     =            \theta_a\theta_o$.      From
Proposition~\ref{proposition:2}                         in
Appendix~\ref{appendix:A},        it     follows      that
$P_{{\db}_b}(\mbox{\tt btype('A')}) = P_{{\db}_b}(S_1 \vee
S_2 \vee S_3) = P_{F_b}(S_1 \vee S_2 \vee S_3)$ and that
%
\begin{eqnarray*}
P_{{\db}_b}(\mbox{\tt btype('A')} \mid \theta_a, \theta_b, \theta_o)
  & = & P_{F_b}(S_1)+P_{F_b}(S_2)+P_{F_b}(S_3)  \\
  & = & {\theta^2_a + 2\theta_a\theta_o}.
\end{eqnarray*}
%
Here  we used  the fact that   $S_1$, $S_2$  and $S_3$ are
mutually  exclusive  as the choice   of gene is exclusive.
Parameters, i.e.\ $\theta_a$, $\theta_b$ and $\theta_o$ are
determined  by ML estimation performed  on a random sample
such as   $\{\mbox{\tt btype('A')}, \mbox{\tt btype('O')},
\mbox{\tt btype('AB')} \}$ of {\tt btype} as follows.
%
\begin{eqnarray*}
\langle \theta_a,\theta_b,\theta_o \rangle
 & = &
  \mbox{argmax}_{\langle \theta_a,\theta_b,\theta_o \rangle}\;
		P_{{\db}_b}(\mbox{\tt btype('A')})
		P_{{\db}_b}(\mbox{\tt btype('O')})
		P_{{\db}_b}(\mbox{\tt btype('AB')}) \\
 & = &
  \mbox{argmax}_{\langle \theta_a,\theta_b,\theta_o \rangle}\;
        (\theta^2_a + 2\theta_a\theta_o)\theta^2_o\theta_a\theta_b
\end{eqnarray*}

This  program   contains neither   function   symbol   nor
recursion though our semantics allows  for them.  Later we
see an  example containing  both,  a  program for an   HMM
\cite{Rabiner93}.

\subsection{Four Simplifying Conditions}
\label{subsec:simplifying-conditions}

${\db}_b$   in  Figure~\ref{program:ABO}  is  simple   and
probability computation is easy. This is not generally the
case.  Since  our primary interest is learning, especially
efficient  parameter  learning  of  parameterized    logic
programs, we   hereafter  concentrate on  identifying what
property of a  program makes  probability computation easy
like ${\db}_b$, thereby makes efficient parameter learning
possible.

To  answer this question  precisely, let  us formulate the
whole   modeling     process.      Suppose   there   exist
symbolic-statistical phenomena   such as gene  inheritance
for which    we    hope  to construct    a   probabilistic
computational  model.  We   first  specify  a  {\em target
predicate\/} $p$  whose ground  atom $p(s)$ represents our
observation  of  the  phenomena.    Then  to  explain  the
empirical  distribution  of  $p$,     we write  down     a
parameterized   logic program $\db  =   F \cup R$ having a
basic distribution $P_F$ with  parameter $\theta$ that can
reproduce all  observable  patterns of  $p(s)$.   Finally,
observing a  random  sample  $p(s_1), \ldots,  p(s_T)$  of
ground atoms of $p$, we adjust  $\theta$ by ML estimation,
i.e.\   by  maximizing   the   likelihood    $L(\theta)  =
\prod_{t=1}^{T}   \pdb(p(s_t) \mid  \theta)$      so  that
$\pdb(p(\cdot)  \mid \theta)$  approximates  as closely to
the empirically observed distribution of $p$ as possible.

At first  sight,  this  formulation  looks right,  but  in
reality it is not.  Suppose two  events $p(s)$ and $p(s')$
$(s\neq s')$ are observed.  We put $L(\theta) = \pdb(p(s)
\mid  \theta)   \pdb(p(s') \mid  \theta)$.   But this {\em
cannot  be\/}  a  likelihood   at  all simply  because  in
distribution  semantics,  $p(s)$   and    $p(s')$ are  two
different random variables,   not two realizations of  the
same random variable.

A quick remedy is to  note that in  the case of blood type
program   $\db_b$   where  $obs(\db_b)    =    \{\mbox{\tt
btype('A')}, \mbox{\tt btype('B')}, \mbox{\tt btype('O')},
\mbox{\tt btype('AB')} \}$ are observable atoms, only one
of them is true for each observation, and  if some atom is
true, others  must be false. In  other words,  these atoms
collectively behave as a single random variable having the
distribution      $P_{{\db}_b}$    whose     values    are
$obs({\db}_b)$.

Keeping   this  in mind,     we  introduce the   following
condition.  Let $obs(\db)$ $(\subset head(R))$ be a set of
ground atoms  which represent observable events.   We call
them {\em observable atom\/}s.
%
\begin{description}
\item[Uniqueness condition:]\mbox{}\\
    $\pdb(G\wedge G') = 0\;\;
           \mbox{for any}\; G\neq G' \in obs(\db)$,
           \hspace{0.3em} and \hspace{0.3em} 
    $\sum_{G\in obs(\db)} \pdb(G) =1$.
\end{description}
%
The uniqueness condition   enables us to introduce   a new
random variable $Y_o$ representing  our  observation.  Fix
an  enumeration   $G_1,G_2,\ldots$ of observable  atoms in
$obs(\db)$ and define $Y_o$ by\footnote{
%
$\sum_{G\in obs(\db)} \pdb(G) =1$ only guarantees that the
measure of  $\{\,\omega \mid  \omega\models G_k\;\mbox{for
some}\;k\,   (\geq 1) \}$  is  one,  so  there can be some
$\omega$ satisfying  no  $G_k$'s.   In  such case,   we put
$Y_o(\omega) = 0$.  But values on a set of measure zero do
not affect any part of the discussion  that follows.  This
also applies to  the definition of $X_e$ in (\ref{eq:Xe}).}


\begin{eqnarray}
   Y_o(\omega) = k & \mbox{iff} & \omega \models G_k
                \;\; \mbox{for}\; \omega \in \Omega_{\db}\; (k \geq 1).
\label{eq:Yo}
\end{eqnarray}
%
Let $G_{k_1},G_{k_2},\ldots,G_{k_T}  \in  obs(\db)$  be  a
random  sample    of  size  $T$.     Then  $L(\theta)   =$
$\prod_{t=1}^{T}    \pdb(G_{k_t}\mid        \theta)     =$
$\prod_{t=1}^{T} \pdb(Y_o  = {k_t}\mid \theta)$  qualifies
for the likelihood function w.r.t.\ $Y_o$.

The second condition concerns the reduction of probability
computation  to  addition.   Take  again   the blood  type
exmaple.    The  computation  of    $P_{{\db}_b}(\mbox{\tt
btype('A')})$   is decomposed   into a  summation  because
explanations in the  support set are mutualy exclusive. So
we introduce
%
\begin{description}
\item[Exclusiveness condition:]\mbox{}\\
    For every  $G \in obs(\db)$ and  the support set
	$\psi_{\db}(G)$,
	$\pdb(S\wedge S') = 0$
		for any $S \neq S'  \in \psi_{\db}(G)$.
\end{description}
%
%
Using    the         exclusiveness       condition    (and
Proposition~\ref{proposition:2}                        in
Appendix~\ref{appendix:A}), we have
%
\[
          \pdb(G) = \sum_{S\in\psi_{\db}(G)}P_F(S).
\]
%
From a  modeling  point of  view,  it means  that while  a
single event,   or a  single observation,   $G$,  may have
several  (or even  infinite) explanations $\psi_{\db}(G)$,
only one of $\psi_{\db}(G)$ is allowed to be true for each
observation.

Now  introduce   $\Psi_{\db}$,  i.e.\  the    set  of  all
explanations relevant to $ obs(\db)$ by
\[
   \Psi_{\db} \defined
                 \bigcup_{G \in obs(\db)} \psi_{\db}(G)
\]
%
and fix  an enumeration $S_1,S_2,\ldots$ of explanations
in    $\Psi_{\db}$.      It       follows            from
Proposition~\ref{proposition:2},  the uniqueness condition
and the exclusiveness condition that
%
\begin{eqnarray*}
\pdb(S_i \wedge S_j)
 & = & 0 \;\; {\mbox{for}}\; i \neq j \hspace{2em}\mbox{and} \\
\sum_{S\in \Psi_{\db}} \pdb(S)
 & = & \sum_{G \in obs(\db)}
			\sum_{S \in \psi_{\db}(G)} \pdb(S) \\
 & = & \sum_{G \in obs(\db)} \pdb(G) \\
 & = & 1.
\end{eqnarray*}
%
So we are able to introduce under the uniqueness condition
and  the    exclusiveness  condition  yet   another random
variable  $X_e$, representing  an  explanation   for  $G$,
defined by
%
\begin{eqnarray}
X_e(\omega) = k & \mbox{iff} & 
  \omega \models S_k  \;\; {\mbox for}\; \omega \in \Omega_{\db}.
\label{eq:Xe}
\end{eqnarray}

The third condition concerns termination.
%
\begin{description}
\item[Finite support condition:]\mbox{}\\
   For every $G\in obs(\db)$  $\psi_{\db}(G)$ is finite.

\end{description}
%
$\pdb(G)$  is   then    computed  from  the   support  set
$\psi_{\db}(G) = \{S_1,\ldots,S_m\}$ $(0\leq m)$, with the
help of the exclusiveness condition, as a finite summation
$\sum_{i=1}^{m}   P_{F}(S_i)$.  This condition prevents an
infinite summation that is hardly computable.

The      fourth  condition   simplifies the    probability
computation to multiplication.  Recall that an explanation
$S$   for   $G\in   obs(\db)$    is  a   conjunction  $a_1
\wedge\cdots\wedge  a_m$     of         some    abducibles
$\{a_1,\ldots,a_m\}  \subset F$ $(1\leq  m)$.  In order to
reduce  the     computation      of     $P_{F}(S)        =
P_{F}(a_1\wedge\cdots\wedge a_m)$ to   the  multiplication
$P_{F}(a_1)\cdots P_{F}(a_m)$, we assume
%
\begin{description}
 \item[Distribution  condition:]\mbox{}\\   $F$  is  a set
     $F_{\tt msw}$ of  ground atoms  with a  parameterized
     distribution $P_{\tt msw}$ specified below.
\end{description}
%
Here atom {\tt msw($i$,$n$,$v$)} is intended to simulate a
multi-ary  random  switch  whose name  is   $i$ and  whose
outcome is  $v$ on trial $n$.   It is a  generalization of
primitive  probabilistic  events such as coin  tossing and
dice rolling.

\begin{enumerate}
\item

$F_{\tt   msw}$  consists  of    probabilistic  atoms {\tt
msw($i$,$n$,$v$)}.  The arguments   $i$, $n$ and $v$   are
ground terms called  {\em switch name}, {\em trial-id} and
a {\em  value\/}  (of the switch   $i$), respectively.  We
assume that a finite set $V_i$ of  ground terms called the
{\em value set of $i$\/}  is associated with each $i$, and
$v\in V_i$ holds.

\item

Write      $V_i$      as   $\{v_1,v_2,\ldots,v_m       \}$
($m=|V_i|$). Then, one   of   the ground  atoms  \{   {\tt
msw($i$,$n$,$v_1$)},  {\tt  msw($i$,$n$,$v_2$)},   \ldots,
{\tt msw($i$,$n$,$v_m$)}\} becomes exclusively true (takes
on value  1)  on  each  trial.    With  each $i$,  a  {\em
parameter\/} $\theta_{i,v}\in [0,1]$ such that $\sum_{v\in
V_i}\theta_{i,v}=1$ is  associated.  $\theta_{i,v}$ is the
probability of   {\tt msw($i$,$\cdot$,$v$)}   being   true
($v\in V_i$).

\item

For each ground terms $i$, $i'$, $n$, $n'$, $v\in V_i$ and
$v'\in  V_{i'}$, random variable {\tt msw($i$,$n$,$v$)} is
independent of {\tt  msw($i'$,$n'$,$v'$)} if $n\neq n'$ or
$i\neq i'$.

\end{enumerate}
%
In  other words,  we introduce  a family of  parameterized
finite distributions $P_{(i,n)}$ such that
%
\begin{eqnarray}
P_{(i,n)}(
	\mbox{\tt msw($i$,$n$,$v_1$)}=x_1,\ldots,
	\mbox{\tt msw($i$,$n$,$v_m$)}=x_m
                   \mid \theta_{i,v_1},\ldots,\theta_{i,v_m})
          \nonumber \\
 \defined
 \left\{
  \begin{array}{lcl}
     \theta_{i,v_1}^{x_1}\cdots \theta_{i,v_m}^{x_m}
         & \mbox{if} & \sum_{k=1}^{m} x_k=1 \\
     0   & \mbox{o.w.}
   \end{array}
   \right.
\label{eq:PF}
\end{eqnarray}
%
where $m=|V_i|$, $x_k  \in \{0,1\}$ $(1\leq k\leq m)$, and
define $P_{\tt msw}$ as their infinite product
%
\[
P_{\tt msw} \defined \prod_{i,n}P_{(i,n)}.
\]
%
Under this condition, we can compute $P_{\tt msw}(S)$, the
probability  of an explanation    $S$, as the product   of
parameters.   Suppose $\mbox{\tt msw($i_j$,$n$,$v$)}$  and
$\mbox{\tt   msw($i_{j'}$,$n'$,$v'$)}$    are    different
conjuncts    in     an explanation        $S=$  $\mbox{\tt
msw($i_1$,$n_1$,$v_1$)}$ $\wedge \cdots \wedge$ $\mbox{\tt
msw($i_k$,$n_k$,$v_k$)}$.  If either $j\neq j'$ or $n \neq
n'$ holds, they are independent  by construction.  Else if
$j=j'$ and $n=n'$ but $v\neq v'$, they are not independent
but  $P_{\tt msw}(S) =  0$ by  construction.  As a result,
whichever condition may hold, $P_{\tt msw}(S)$ is computed
from the parameters.

\subsection{Modeling Principle}
\label{subsec:modeling-principle}

Up to this point, we  have introduced four conditions, the
uniqueness condition,  the  exclusiveness condition,   the
finite  support condition  and the distribution condition,
to simplify probability computation.  The last one is easy
to satisfy.  We   just adopt $F_{\tt msw}$  together  with
$P_{\tt  msw}$.  So, from here on,   we always assume that
$F_{\tt  msw}$  has a parameterized  distribution  $P_{\tt
msw}$ introduced      in        the previous   subsection.
Unfortunately the   rest are not  satisfied automatically.
According to our  modeling experiences however, it is only
mildly difficult to satisfy  the uniqueness  condition and
the  exclusiveness  condition as    long  as we  obey  the
following modeling principle.

\begin{quotation}
{\bf  Modeling   principle:} $\db =   F_{\tt  msw} \cup R$
describes  a sequential decision process (modulo auxiliary
computations) that uniquely produces an observable atom $G
\in obs(\db)$ where each decision is expressed by some {\tt
msw} atom.\footnote{
%
Decisions  made in  the process  are   a finite subset  of
$F_{\tt msw}$. }
%
\end{quotation}
%

Translated into  programming level, it   says that we must
take care when   writing a program so  that   {\em for any
sample $F'$ from $P_{\tt  msw}$, there must uniquely exist
goal  $\leftarrow\!  G$ $(G  \in   obs(\db))$ which has  a
successful  refutation  from   $\db'=F'\cup R$}.   We  can
confirm  the principle by the  blood type program $\db_b =
F_b   \cup R_b$.    It   describes   a  process   of  gene
inheritance, and   for  an arbitrary  sample  $F_b'$  from
$P_{\tt     msw}$,   say     $F_b'    =    \{    \mbox{\tt
msw(gene,father,a)},\mbox{\tt  msw(gene,mother,o)}    \}$,
there exists   a  unique goal,   $\leftarrow\!   \mbox{\tt
btype('A')}$  in  this  case, that  has  a  successful SLD
refutation from $F_b' \cup R_b$.

The idea behind this  principle is that a decision process
always produces  some  result  (an observable  atom),  and
different decision processes must differ at some {\tt msw}
thereby entailing  mutually exclusive observable atoms. So
the  uniqueness condition and  the exclusiveness condition
will be automatically satisfied.

Satisfying the finite  support condition is more difficult
as it is virtually equivalent  to writing a program  $\db$
for which all solution  search for $\leftarrow\!G$ $(G \in
obs(\db))$  always   terminates.  Apparently   we  have no
general solution to  this problem, but  as far as specific
models  such  as  HMMs,  PCFGs  and Bayesian  networks are
concerned, it can be met.   All programs for these  models
satisfy the finite support condition (and other conditions
as well).

\subsection{Four Conditions Revisited}
\label{subsec:revisited}

In  this  subsection, we   discuss how  to relax the  four
simplifying            conditions          introduced   in
Subsection~\ref{subsec:simplifying-conditions}   for   the
purpose of   flexible  modeling.   We  first   examine the
uniqueness  condition considering its  crucial role in the
adaptation of the EM algorithm to our semantics.

The  uniqueness condition  guarantees that there  exists a
(many-to-one) mapping from explanations to observations so
that the EM algorithm is applicable \cite{Dempster77}.  It
is possible,   however, to relax  the uniqueness condition
while justifying the application  of the EM algorithm.  We
assume the  {\em  MAR  (missing  at  random)  condition\/}
introduced  by \citeA{Rubin76}  which  is  a   statistical
condition on how a complete  data (explanation) becomes an
incomplete data (observation),  and is customarily assumed
implicitly  or     explicitly    in   statistics      (see
Appendix~\ref{appendix:B}).       By assuming  the     MAR
condition, we can apply  our EM algorithm to non-exclusive
observations  $O$ such that  $\sum_{O} P(O)  \geq 1$ where
the uniqueness condition is seemingly destroyed.

Let  us see  the  MAR  condition in action   with a simple
example.  Imagine we walk along a road in front of a lawn.
We occasionally observe their state  such as ``the road is
dry but  the lawn  is   wet''.  Assume  that the  lawn  is
watered by a   sprinkler running probabilistically.    The
program ${\db}_{\tt rl}$ $=  R_{\tt rl} \cup F_{\tt rl}  $
in  Figure~\ref{prog:wet-dry}  describes   a    sequential
process        which     outputs   an   observation   {\tt
observed(road($x$),lawn($y$))} (``the  road is $x$ and the
lawn is $y$'') where $x, y \in \{{\tt wet,dry}\}$.

\begin{figure}[hbt]
\hspace*{3em}
\begin{tabular}{rcl}
$R_{\tt  rl}$	& = 
	& {\verb! { observed(road(X),lawn(Y)):-!} \\
&	& \hspace{3em}{\verb!msw(rain,once,A),!} \\
&	& \hspace{3em}{\verb!( A = yes, X = wet, Y = wet!} \\
&	& \hspace{3em}{\verb!; A = no,  msw(sprinkler,once,B),!} \\
&	& \hspace{3em}{\verb!     ( B = on, X = dry, Y = wet!} \\
&	& \hspace{3em}{\verb!     ; B = off, X = dry, Y = dry ) ).  }!} \\
%
$F_{\tt  rl}$	& = 
	& {\verb! { msw(rain,once,yes), msw(rain,once,no),!} \\
&	& {\verb!      msw(sprinkler,once,on), msw(sprinkler,once,off) }!}
\end{tabular}

\caption{${\db}_{\tt rl}$}
\label{prog:wet-dry}
\end{figure}


The basic distribution over $F_{\tt rl}$ is specified like
$P_{F_b}(\cdot)$ in  Subsection~\ref{subsec:DB-b},   so we
omit it.  {\tt msw(rain,once,A)} in the program determines
whether it rains (${\tt A} =  {\tt yes}$) or not (${\tt A}
=   {\tt  no}$),   whereas   {\tt   msw(sprinkler,once,B)}
determines  whether  the sprinkler  works  fine  (${\tt B}
={\tt on}$) or not (${\tt B} ={\tt off}$).  Since for each
sampled  values of ${\tt A}  = a$ $(a\in \{ {\tt yes},{\tt
no}  \})$ and  ${\tt  B}=b$ $(b\in  \{  {\tt on},{\tt off}
\})$,   there   uniquely   exists   an    observation {\tt
observed(road($x$),lawn($y$))}      ($x,y\in$      $\{{\tt
wet,dry}\}$), there    is  a  many-to-one  mapping  $\chi:
\chi(a,b) = \langle x,y \rangle$.  In  other words, we can
apply  the  EM   algorithm   to the    observations   {\tt
observed(road($x$),lawn($y$))}   ($x,   y      \in  \{{\tt
wet,dry}\}$).  What would happen if we observe exclusively
either  a  state   of the  road   or that   of  the  lawn?
Logically, this means  we  observe $\exists  y\, \mbox{\tt
observed(road($x$),lawn($y$))}$ or $\exists x\,  \mbox{\tt
observed(road($x$),lawn($y$))}$.       Apparently      the
uniqueness condition is not   met, because $\exists   y \,
\mbox{\tt observed(road(wet),lawn($y$))}$ and $\exists x\,
\mbox{\tt observed(road($x$),lawn(wet))}$  are  compatible
(they   are   true   when   it    rains).   Despite    the
non-exclusiveness of the observations,  we can still apply
the EM algorithm to them under the MAR condition, which in
this case translates into  that we observe either the lawn
or the road {\em randomly regardless of their state}.\\

We now briefly check other conditions.  Basically they can
be relaxed at the  cost of increased computation.  Without
the exclusiveness condition for instance, we would need an
additional process   of    transforming  the support   set
$\psi_{\db}(G)$  for a  goal $G$ into   a set of exclusive
explanations.  For  instance, if $G$ has  explanations $\{
\mbox{\tt           msw($a$,$n$,$v$)},$         $\mbox{\tt
msw($b$,$m$,$w$)}\}$,  we  have    to   transform it  into
$\{\mbox{\tt    msw($a$,$n$,$v$)},$    $\neg     \mbox{\tt
msw($a$,$n$,$v$)} \wedge \mbox{\tt  msw($b$,$m$,$w$)}  \}$
and so on.\footnote{
%
$\neg  \mbox{\tt msw($a$,$n$,$v$)}$ is further transformed
to  a   disjunction  of exclusive  {\tt   msw} atoms  like
$\bigvee_{v'     \neq      v, v'\in        V_a}  \mbox{\tt
msw($a$,$n$,$v'$)}$.  }
%
Clearly, this  transformation is exponential in the number
of  {\tt msw}   atoms   and efficiency concern leads    to
assuming the exclusiveness condition.

The finite support  condition is in practice equivalent to
the  condition  that the  SLD  tree for $\leftarrow  G$ is
finite. So  relaxing this  condition might induce infinite
computation.

Relaxing  the   distribution   condition    and  accepting
probability distributions other  than $P_{\tt msw}$  serve
to   expand   the  horizon  of     the   applicability  of
parameterized    logic programs.    In particular      the
introduction   of parameterized    joint     distributions
$P(v_1,\ldots,v_k)$  like   Boltzmann  distributions  over
switches   ${\tt  msw}_1,  \ldots,    {\tt msw}_k$   where
$v_1,\ldots,v_k$ are values of   the switches, makes  them
{\em  correlated}.   Such distributions facilitate writing
parameterized  logic   programs  for complicated  decision
processes in  which  decisions  are not   independent  but
interdependent.  Obviously,  on  the  other   hand,   they
increase learning  time, and whether the added flexibility
of distributions  deserves the increased learning  time or
not is yet to be seen.


\subsection{Naive Approach to EM Learning}
\label{subsec:learning:naive}

In this subsection, we derive  a concrete EM algorithm for
parameterized logic programs  $\db  = F_{\tt   msw}\cup R$
assuming that  they satisfy  the uniqueness condition, the
exclusiveness condition and the finite support condition.

To start,   we    introduce  $Y_o$, a    random   variable
representing our observations according to   (\ref{eq:Yo})
based on  a fixed  enumeration   of  observable atoms   in
$obs(\db)$.  We   also  introduce another random  variable
$X_e$ representing   their   explanations   according   to
(\ref{eq:Xe})  based    on    some  fixed   enumeration of
explanations in $\Psi_{\db}$.   Our understanding  is that
$X_e$  is non-observable while  $Y_o$  is observable,  and
they have a joint distribution $\pdb(X_e = x, Y_o = y \mid
\theta)$  where $\theta$ denotes  relevant parameters.  It
is   then    immediate,  following     (\ref{eq:Q1})   and
(\ref{eq:Q2})     in Section~\ref{sec:preliminaries},   to
derive a concrete EM algorithm from the Q function defined
by
%
$ Q(\theta \mid \theta') \defined
    \sum_x  \pdb(x\mid y,\theta') \ln \pdb(x, y \mid \theta)$
%
whose input  is a  random  sample of  observable atoms and
whose output is the MLE of $\theta$.

In  the  following,   for the  sake   of  readability,  we
substitute  an observable atom  $G$ $(G \in obs(\db))$ for
$Y_o =y$ and  write    $\pdb(G \mid \theta)$   instead  of
$\pdb(Y_o=y  \mid   \theta)$.   Likewise we  substitute an
explanation $S$ $(S \in \Psi_{\db})$ for $X_e=x$ and write
$\pdb(S, G  \mid  \theta)$  instead of  $\pdb(X_e=x, Y_o=y
\mid  \theta)$.   Then  it  follows   from the  uniqueness
condition that
%
\begin{eqnarray*}
\pdb(S,G \mid \theta)
 & = & \left\{ \begin{array}{lcl}
                           0  & \mbox{if} & S \not\in \psi_{\db}(G) \\
   P_{\tt msw}(S \mid \theta) & \mbox{if} & S \in \psi_{\db}(G).
               \end{array}
       \right.
\end{eqnarray*}
%
We need yet another   notation here.  For  an  explanation
$S$,     define      the   {\em     count   of  $\mbox{\tt
msw($i$,$n$,$v$)}$ in $S$\/} by
%
\[
\sigma_{i,v}(S)\defined
   \left| \left\{\; n\; \big| \;\mbox{\tt msw($i$,$n$,$v$)} \in S
	\right\}\right|.
\]
%
We  have done all preparations  now.  Suppose we make some
observations ${\cal  G}=G_1,\ldots,G_T$   where $G_t   \in
obs(\db)$ $(1\leq t \leq T)$. Put
%
\begin{eqnarray*}
I & \defined &
    \{ i \mid 
           \mbox{\tt msw($i$,$n$,$v$)}
				\in S \in \psi_{\db}(G_t), 1\leq t \leq T
    \} \\
\btheta
 & \defined &
    \{ \theta_{i,v} \mid 
           \mbox{\tt msw($i$,$n$,$v$)}
				\in S\in \psi_{\db}(G_t), 1\leq t \leq T
    \}.
\end{eqnarray*}
%
$I$  is   a  set  of  switch  names  that  appear  in some
explanation for one of  the $G_t$'s and $\btheta$  denotes
parameters associated with  these  switches.  $\btheta$ is
finite due to the finite support condition.\\

Various  probabilities and the Q  function are computed by
using           Proposition~\ref{proposition:2}         in
Appendix~\ref{appendix:A} together with our assumptions as
follows.
%
\begin{eqnarray}
\pdb(G_t \mid \btheta)
 & = & \pdb\left( \left. \bigvee \psi_{\db}(G_t)\; \right|
                                                   \btheta  \right)
       \hspace{1em}
    =  \sum_{S\in\psi_{\db}(G_t)}P_{\tt msw}(S \mid  \btheta)
      \label{eq:PDB} \\
P_{\tt msw}(S \mid  \btheta)
 & = &  \prod_{i\in I, v\in V_i}
			         \theta_{i,v}^{\sigma_{i,v}(S)}    \nonumber \\
Q(\btheta \mid \btheta')
  & \defined &
   \sum_{t =1}^{T} \sum_{S\in \Psi_{\db}} 
      \pdb(S \mid G_t, \btheta') \ln \pdb(S, G_t \mid \btheta)
                                     \nonumber \\
  & = &
	\sum_{i\in I, v\in V_i}
		\eta(i,v,\btheta')\ln\theta_{i,v}
	\leq
	\sum_{i\in I, v\in V_i}
		\left(
			\eta(i,v,\btheta')\ln
				\frac{\eta(i,v,\btheta')}
						{\sum_{v'\in V_i}\eta(i,v',\btheta')}
		\right) \label{eq:eta}\\
\mbox{where}\hspace{2em} & \nonumber \\
\eta(i,v,\btheta)
  & \defined &
	\sum_{t=1}^T \frac{1}{\pdb(G_t \mid \btheta)}
		\sum_{S\in\psi_{\db}(G_t)}
			P_{\tt msw}(S \mid \btheta) \sigma_{i,v}(S)  \nonumber
\end{eqnarray}
%
Here we used Jensen's inequality to obtain (\ref{eq:eta}).
Note that  ${\pdb(G_t  \mid  \btheta)}^{-1}$  $\sum_{S \in
\psi_{\db}(G_t)}$    $P_{\tt     msw}(S   \mid    \btheta)
\sigma_{i,v}(S)$  is    the {\em  expected   count\/}   of
$\mbox{\tt  msw($i$,$\cdot$,$v$)}$ in an SLD refutation of
$G_t$.  Speaking of the likelihood function $L(\btheta) =
\prod_{t=1}^{T} \pdb(G_t  \mid \btheta)$,  it  is  already
shown   in  Subsection~\ref{sub:EM}   (footnote)      that
$Q(\btheta  \mid  {\btheta'})    \geq  Q({\btheta}'   \mid
{\btheta}')$    implies  $L(\btheta) \geq  L({\btheta}')$.
Hence   from  (\ref{eq:eta}),  we    reach  the  procedure
\pname{learn-naive($\db$,${\cal G}$)} below that finds the
MLE of  the parameters.    The array variable  $\eta[i,v]$
stores $\eta(i,v,\btheta)$ under the current $\btheta$.

\begin{listing}
\samepage
\item \smbf{procedure} \pname{learn-naive($\db,{\cal G}$)} \smbf{begin}
\item  Initialize $\btheta$ with appropriate values
       and $\varepsilon$ with a small positive number ;
\item
       $\lambda^{(0)}:=\sum_{t=1}^T\ln \pdb(G_t \mid \btheta)$;
              \hspace{3.5em}{\% Compute the log-likelihood.}
	\label{line:learn-naive:get-loglike1}
\item \smbf{repeat}
\item \hspace{1em} \smbf{foreach} $i\in I,v\in V_i$ \smbf{do}
\item \hspace{2em}
        $\displaystyle \eta[i,v]:=
		  \sum_{t=1}^T\frac{1}{\pdb(G_t \mid \btheta)}
		  \sum_{S\in \psi_{\db}(G_t)}
			P_{\tt msw}( S \mid \btheta)\sigma_{i,v}(S)$;
      \label{line:learn-naive:get-On}
%\halfempline
\item \hspace{1em} \smbf{foreach} $i\in I,v\in V_i$ \smbf{do}
\item \hspace{2em}
         $\displaystyle \theta_{i,v}
            \,:=\, \frac{\eta[i,v]}{\sum_{v'\in V_i}\eta[i,v']}$;
		\hspace{4em}{\% Update the parameters.}
%\halfempline
\item  \hspace{1em} $m := m+1$;
\item  \hspace{1em}
           $\lambda^{(m)} \,:=\, \sum_{t=1}^T\ln \pdb(G_t \mid \btheta)$
		\hspace{2em}{\% Compute the log-likelihood again.}
     \label{line:learn-naive:get-loglike2}
\item \smbf{until} $\lambda^{(m)}-\lambda^{(m-1)}<\varepsilon$
		\hspace{5.6em}{\% Terminate if converged.}
\item \smbf{end}
\end{listing}

This  EM algorithm is  simple and correctly calculates the
MLE of $\btheta$,  but the  calculation of $\pdb(G_t  \mid
\btheta)$                  and                  $\eta[i,v]
$(Line~\ref{line:learn-naive:get-loglike1},
\ref{line:learn-naive:get-On}                          and
\ref{line:learn-naive:get-loglike2})   may      suffer   a
combinatorial   explosion   of  explanations.    That  is,
$|\psi_{\db}(G_t)|$   often grows   exponentially in   the
complexity   of      the    model.    For        instance,
$|\psi_{\db}(G_t)|$ for   an  HMM   with  $N$  states   is
$O(N^L)$, exponential in the length $L$ of an input/output
string.  Nonetheless, suppressing the explosion to realize
efficient  computation in a  polynomial order is possible,
under  suitable    conditions,    by   avoiding   multiple
computations of the same subgoal as we see next.

\subsection{Inside Probability and Outside Probability for
 Logic Programs}
\label{subsec:io-prob}

In this subsection,  we  generalize the notion of   inside
probability and outside probability \cite{Baker79,Lari90}
to    logic     programs.      Major   computations     in
\pname{learn-naive($\db$,${\cal  G}$)}  are  those of  two
terms   in   Line~\ref{line:learn-naive:get-On}, $\pdb(G_t
\mid   \btheta)$ and $\sum_{S  \in \psi_{\db}(G_t)} P_{\tt
msw}(S  \mid    \btheta)  \sigma_{i,v}(S)$.  Computational
redundancy lurks in  the naive computation  of both terms.
We  show   it   by   an example.   Suppose    there  is  a
propositional program $\db_p = F_p \cup R_p$ where $F_p =
\{ {\tt a,b,c,d,m} \}$ and
%
\begin{equation}
R_p = \left\{
\begin{array}{l}
   {\tt f \leftarrow  a \wedge g} \\
   {\tt f \leftarrow  b \wedge g} \\
   {\tt g \leftarrow  c} \\
   {\tt g \leftarrow  d \wedge h} \\
   {\tt h \leftarrow  m}. \\
\end{array}
     \right.
              \label{prog:p}
\end{equation}
%
Here {\tt f}  is an observable  atom.  We assume that {\tt
a}, {\tt b}, {\tt c}, {\tt  d} and {\tt m} are independent
and  also that  $\{{\tt  a,b}\}$  and $\{{\tt  c,d}\}$ are
pair-wise exclusive.  Then the  support set for {\tt f} is
calculated as
%
\[
\psi_{{DB}_p}({\tt f}) =
      \{ {\tt a \wedge c},\; {\tt a \wedge d \wedge m},\;
          {\tt b \wedge c},\; {\tt b \wedge d \wedge m}\;
      \}.
\]
%
Hence, in   light   of  (\ref{eq:PDB}),   we   may compute
$P_{{DB}_p}({\tt f})$ as
%
\begin{equation}
P_{{DB}_p}({\tt f})
      =   P_{F_p}({\tt a \wedge c}) + 
          P_{F_p}({\tt a \wedge d \wedge m}) + 
          P_{F_p}({\tt b \wedge c}) +
          P_{F_p}({\tt b \wedge d \wedge m}).
                   \label{eq:PDBP}
\end{equation}
%
This  computation   requires  6 multiplications   (because
$P_{F_p}({\tt a \wedge c}) = P_{F_p}({\tt a}) P_{F_p}({\tt
c})$ etc.) and 3  additions.   On  the  other hand,  it  is
possible  to compute   $P_{{DB}_p}({\tt  f})$  much   more
efficiently by factoring out common computations.  Let $A$
be a ground  atom.  Define the  {\em inside probability\/}
$\beta(A)$ of $A$ as
%
\begin{equation}
 \beta(A) \defined \pdb(A \mid \btheta).
      \label{eq:inside-def}
\footnote{
%
Note that  if $A$  is  a fact  in $F$, $\beta(A)  = P_{\tt
msw}(A \mid \theta)$.
%
}
\end{equation}
%
Then       by   applying     Theorem~\ref{th:1}         in
Appendix~\ref{appendix:A} to
%
\begin{equation}
comp(R_p) \vdash {\tt  f \leftrightarrow (a \wedge g)\vee (b \wedge g)},\;
                 {\tt  g \leftrightarrow  c \vee (d \wedge h)},\;
                 {\tt  h \leftrightarrow  m}
      \label{eq:p-iff}
\end{equation}
%
which unconditionally holds in our semantics, and by using
the independent and  the exclusiveness assumption  made on
$F_p$,  the  following equations about  inside probability
are derived.
%
\begin{equation}
\left\{
\begin{array}{lll}
  \beta({\tt f})
     & = &  \beta({\tt a})\beta({\tt g})+ \beta({\tt b})\beta({\tt g}) \\
  \beta({\tt g})
     & = &  \beta({\tt c})+ \beta({\tt d})\beta({\tt h}) \\
  \beta({\tt h})             
     & = &  \beta({\tt m})
\end{array}
\right.
      \label{eq:inside-example}
\end{equation}
%
$P_{{DB}_p}({\tt f}) (  =\beta({\tt  f}))$ is obtained  by
solving (\ref{eq:inside-example}) about  $\beta({\tt f})$,
for which  only   3 multiplications and 2    additions are
required.

It       is quite     straightforward     to    generalize
(\ref{eq:inside-example}) but before proceeding, look at a
program $\db_q=$  $\{ {\tt  m}\}  \cup \{\mbox{\tt   g:-m}
\wedge {\tt m}, \; \mbox{\tt g:-m}\}$  where {\tt g} is an
observable atom and  {\tt m} the  only {\tt  msw} atom. We
have ${\tt g}\!\leftrightarrow ({\tt m}\wedge{\tt m}) \vee
{\tt  m}$ in our semantics, but  to  compute $P({\tt g}) =
P({\tt m})P({\tt m}) + P({\tt  m})$ is clearly wrong as it
ignores the fact  that clause  bodies  for {\tt g},  i.e.\
{\tt m$\wedge$m} and  {\tt m} are  not mutually exclusive,
and  atoms  in the clause  body  {\tt  m$\wedge$m} are not
independent      (here    $P(\cdot) =  P_{\db_q}(\cdot)$).
Similarly, if we set {\tt a} = {\tt b} = {\tt c} = {\tt d}
=  {\tt m}, the equation~(\ref{eq:inside-example}) will be
totally incorrect.

We therefore add, {\em  temporarily in this subsection\/},
two assumptions on top  of the exclusiveness condition and
the   finite support   condition  so   that equations like
(\ref{eq:inside-example})  become mathematically  correct.
The  first   assumption  is   that  ``clause'' bodies  are
mutually exclusive   i.e.\    if there are   two   clauses
$B\leftarrow  W$  and $B\leftarrow W'$, $\pdb(W  \wedge W'
\mid \btheta) = 0$, and the second assumption is that body
atoms  are independent, i.e.\ if  $A \leftarrow B_1 \wedge
\cdots \wedge B_k$   is a  rule, $\pdb(B_1 \wedge   \cdots
\wedge B_k \mid \btheta) =$ $\pdb(B_1 \mid\btheta)\cdots
\pdb$ $(B_k \mid \btheta)$ holds.

Please note that ``clause'' used in  this subsection has a
special  meaning. It is   intended  to mean $G  \leftarrow
\tau$ where  $G$ is  a goal and   $\tau$ is a  {\em tabled
explanation\/} for  $G$  obtained by  {\em  OLDT search\/}
both   of   which  will   be    explained  in  the    next
subsection.\footnote{
%
The logical  relationship (\ref{eq:p-iff}) corresponds  to
(\ref{eq:tabled-completion})  where {\tt  f}, {\tt  g} and
{\tt h} are table atoms.  }
%
In other   words, these   additional conditions   are  not
imposed   on a source program  but   on the result of OLDT
search.  So clauses for auxiliary computations do not need
to satisfy them.\\

Now suppose clauses about  $A$ occur in $\db$ like
%
\[
\begin{array}{l}
A \leftarrow B_{1,{1}} \wedge \cdots \wedge B_{1,{i_1}} \\
\hspace{2em}\cdots \\
A \leftarrow  B_{L,{1}} \wedge \cdots \wedge B_{L,{i_L}}
\end{array}
\]
%
where $B_{h,j}$ $(1\leq h \leq L,  1\leq j\leq i_h)$ is an
atom.  Theorem~\ref{th:1} in Appendix~\ref{appendix:A} and
the above assumptions ensure
%
\begin{equation}
\beta(A)  =  \prod_{j=1}^{i_1} \beta(B_{1,j}) +  \cdots + 
                     \prod_{j=1}^{i_L}  \beta(B_{L,j}).
\label{eq:inside-rec}
\end{equation}
%
(\ref{eq:inside-rec})  suggests that  $\beta(G_t)$ can  be
considered as a function  of $\beta(A)$ if these equations
about inside probabilities are hierarchically organized in
such a way that $\beta(G_t)$  belongs to the top layer and
any $\beta(A)$ appearing on the left hand side only refers
to $\beta(B)$'s   which belong to  the lower   layers.  We
refer to    this condition as  the    {\em acyclic support
condition}.   Under  the  acyclic  support      condition,
equations of the  form (\ref{eq:inside-rec}) have a unique
solution, and  the computation  of $\pdb(G \mid  \btheta)$
via inside  probabilities allows  us to take  advantage of
reusing intermediate results stored as $\beta(A)$, thereby
contributing  to  faster  computation  of   $\pdb(G_t \mid
\btheta)$.\\

Next we tackle  a more intricate problem, the  computation
of      $\sum_{S \in  \psi_{\db}(G_t)}P_{\tt  msw}(S \mid
\btheta) \sigma_{i,v}(S)$.  Since  the sum equals  $\sum_n
\sum_{\mbox{\footnotesize\tt  msw($i$,$n$,$v$)}  \in S \in
\psi_{\db}(G_t)}P_{\tt   msw}(S  \mid  \btheta)$,      we
concentrate on the computation of
%
\[
 \mu(G_t,{\tt m}) \defined \sum_{{\tt m} \in S \in \psi_{\db}(G_t)}
  P_{\tt msw}(S \mid \btheta)
\]
%
where ${\tt m}  = \mbox{\tt msw($i$,$n$,$v$)}$.  First  we
note that if an explanation $S$ contains {\tt m} like $S =
a_1 \wedge \cdots \wedge a_h \wedge {\tt m}$, then we have
$\beta(S)  = \beta(a_1) \cdots \beta(a_h) \beta({\tt m})$.
So $\mu(G_t,{\tt m})$ is expressed as
%
\begin{equation}
\mu(G_t,{\tt   m}) = \alpha(G_t,{\tt  m})\beta({\tt  m})
                                 \label{eq:outside-m}
\end{equation}
%
where  $\alpha(G_t,{\tt m}) = \frac{\partial  \mu(G_t,{\tt
m})} {\partial \beta({\tt  m})}$ and $\alpha(G_t,{\tt m})$
does not  depend  on $\beta({\tt m})$.  Generalizing  this
observation to  arbitrary  ground atoms,  we introduce the
{\em outside probability\/}  of  ground atom  $A$  w.r.t.\
$G_t$ by
%
\[
 \alpha(G_t,A)  \defined
            \frac{\partial \beta(G_t)}{\partial \beta(A)}
\]
%\label{eq:outside-def}
%
assuming the same conditions  as  inside probability.   In
view  of  (\ref{eq:outside-m}), the  problem  of computing
$\mu(G_t,{\tt m})$  is  now reduced to  that  of computing
$\alpha(G_t,{\tt m})$, which  is recursively computable as
follows.  Suppose $A$ occurs  in the ground program  $\db$
like
%
\[
\begin{array}{l}
B_1 \leftarrow  A \wedge W_{1,{1}},
         \cdots, B_1 \leftarrow  A \wedge W_{1,{i_1}} \\
\hspace{2em}\cdots \\
B_K \leftarrow  A \wedge W_{K,{1}},
         \cdots, B_K \leftarrow  A \wedge W_{K,{i_K}}.
\end{array}
\]
%
As   $\beta(G_t)$  is a function  of  $\beta(B_1), \ldots,
\beta(B_K)$  by  our   assumption,  the   chain   rule  of
derivatives leads to
\[
\alpha(G_t,A)   =
  \left(\frac{\partial \beta(G_t)}{\partial \beta(B_1)}\right)
  \left(\frac{\partial \beta(A\wedge W_{1,1})}
                       {\partial \beta(A)}\right)
  + \cdots +
  \left(\frac{\partial \beta(G_t)}{\partial \beta(B_K)}\right)
  \left(\frac{\partial \beta(A\wedge 	W_{K,i_K})}
             {\partial \beta(A)}\right)
\]
and hence to\footnote{
%
Because of  the independence  assumption  on body   atoms,
$W_{h,j}$ $(1\leq h \leq K,  1\leq j\leq i_h)$ and $A$ are
independent.  Therefore
\[
\frac{\partial \beta(A\wedge W_{h,j})}{\partial \beta(A)}
=
\frac{\partial \beta(A)\beta(W_{h,j})}{\partial \beta(A)}
       = \beta(W_{h,j}).
\]
}
%
\begin{eqnarray}
\alpha(G_t,G_t) & = &  1
                  \label{eq:outside-ini} \\
\alpha(G_t,A) & = &  \alpha(G_t,B_1)\sum_{j=1}^{i_1}\beta(W_{1,j})
                     + \cdots +
                   \alpha(G_t,B_K)\sum_{j=1}^{i_K}\beta(W_{K,j}).
   \label{eq:outside-rec}
\end{eqnarray}
%
Therefore if  all inside  probabilities have already  been
computed, outside probabilities are recursively   computed
from    the      top       (\ref{eq:outside-ini})    using
(\ref{eq:outside-rec}) downward along the program  layers.
In the case of ${\db}_p$  with {\tt f}  and {\tt m}  being
chosen atoms, we compute
%
\begin{equation}
\left\{
\begin{array}{lll}
  \alpha({\tt f},{\tt f})
     & = &  1 \\
  \alpha({\tt f},{\tt g})
     & = &  \beta({\tt a})+ \beta({\tt b}) \\
  \alpha({\tt f},{\tt h})             
     & = &  \alpha({\tt f},{\tt g})\beta({\tt d}) \\
  \alpha({\tt f},{\tt m})             
     & = &  \alpha({\tt f},{\tt h}).
\end{array}
\right.
        \label{eq:outside-example}
\end{equation}
%
From (\ref{eq:outside-example}), the desired sum $\mu({\tt f},{\tt
m})$ is calculated as
%
\[ \mu({\tt f},{\tt m})  =  \alpha({\tt f},{\tt  m})\beta({\tt m})
                         =  (\beta({\tt a})+ \beta({\tt b}))
                             \beta({\tt d})\beta({\tt m})
\]
%
which requires only  two multiplications and one  addition
compared to four multiplications  and one addition  in the
naive computation.

Gains obtained by computing inside and outside probability
may be small for this case, but as the problem size grows,
they become enormous, and compensate enough for additional
restrictions imposed on the result of OLDT search.

\subsection{OLDT Search}
\label{subsec:OLDT}

To compute inside and outside probability recursively like
(\ref{eq:inside-rec})   or     (\ref{eq:outside-ini})  and
(\ref{eq:outside-rec}),   we need at  programming  level a
tabulation    mechanism  for structure-sharing of  partial
explanations  between subgoals.  We  henceforth deal  with
programs $\db$ in which a  set $table(\db)$ of {\em  table
predicate\/}s are   declared in  advance.   A ground  atom
containing a table predicate is called a {\em table atom}.
The purpose of table atoms is  to store their support sets
and eliminate the need of recomputation,  and by doing so,
to construct hierarchically organized explanations made up
of the table atoms and the {\tt msw} atoms.

Let ${\db} = F_{\tt msw}  \cup R$ be a parameterized logic
program which  satisfies the finite  support condition and
the    uniqueness            condition.         Also   let
$G_{1},G_{2},\ldots,G_{T}   $ be    a  random sample    of
observable  atoms in $obs(\db)$.    We make  the following
additional assumptions.\\

\noindent
{\bf Assumptions:}

For each $t$ $(1\leq t \leq T)$, there exists a finite set
$\{\tau_1^{t},  \ldots, \tau_{K_t}^{t}\}$  of  table atoms
associated with conjunctions $\Stil_{k,j}^{t}$ $(0 \leq k
\leq K_t, 1\leq j \leq m_k)$ such that
%
\begin{equation}
\begin{array}{lll}
comp(R)
& \vdash &	\left(G_t\leftrightarrow
							\Stil_{0,1}^{t}\lor\cdots\lor
							\Stil_{0,m_0}^{t}\right)\\
&        &	\wedge	\left(\tau_1^{t}\leftrightarrow
					\Stil_{1,1}^{t}\vee	\cdots \vee
					\Stil_{1,m_1}^{t}\right)
			\wedge
			\cdots
			\wedge  \left(\tau_{K_t}^{t}\leftrightarrow
					\Stil_{K_t,1}^{t} \vee \cdots \vee
					\Stil_{K_t,m_{K_t}}^{t}\right)
%	\nonumber
\end{array}
	\label{eq:tabled-completion}
\end{equation}
%
where
\begin{itemize}

\item

each  $\Stil_{k,j}^{t}$ $(0 \leq k \leq  K_t, 1\leq j \leq
m_k)$ is, as a set, a subset of $F_{\tt msw} \cup
\{\tau_{k+1}, \ldots, \tau_{K_t}\}$ {\bf (acyclic  support
condition)}.  As  a convention, we put  $\tau_0 = G_t$ and
call     respectively        $\taudb^{t}          \defined
\{\tau_0,\tau_1^{t},  \ldots, \tau_{K_t}^{t}\}$  the  {\em
set of table atoms for  $G_t$\/} and $\Stil_{k,j}^{t}$ $(k
\geq  0)$ a {\em t-explanation for $\tau_k^{t}$}.\footnote
{
%
Prefix ``t-'' is an abbreviation of ``tabled-''.}
%
The set  of all t-explanations  for $\tau_k$ is denoted by
$\psidbtilde(\tau_k^{t})$       and     we        consider
$\psidbtilde(\cdot)$ as a function of table atoms.

\item

t-explanations are mutually  exclusive, i.e.\ for each $k$
$(0   \leq   k  \leq  K_t)$,  $\pdb(\Stil_{k,j}^{t} \wedge
\Stil_{k,j'}^{t})  = 0$ $(1\leq j  \neq j' \leq m_k)$ {\bf
($t-$exclusiveness condition)}.

\item

$\Stil_{k,j}^{t}$ $(0 \leq k \leq K_t,  1\leq j \leq m_k)$
is  a conjunction of   independent atoms {\bf (independent
condition)}.\footnote{
%	
The independence  mentioned  here  only  concerns positive
propositions.   For  $B_1,B_2\in  head(\db)$, we say $B_1$
and $B_2$  are  independent if $\pdb(B_1   \wedge B_2 \mid
\theta)= \pdb(B_1  \mid  \theta)\pdb(B_2 \mid \theta)$ for
any $\theta$.}
%

\end{itemize}

These  assumptions are   aimed   at efficient  probability
computation.  Namely, the  acyclic support condition makes
dynamic  programming   possible,     the   t-exclusiveness
condition reduces   $\pdb(A \vee B)$  to $\pdb(A)+\pdb(B)$
and the independent condition  reduces $\pdb(A \wedge  B)$
to  $\pdb(A)\pdb(B)$.  There is  one more point concerning
efficiency however.  Note that  the computation in dynamic
programming proceeds following  the {\em partial  order\/}
on $\taudb^{t}$\footnote{
%
{\em $\tau_i$   precedes $\tau_j$\/}  if and only   if the
top-down execution   of  $\tau_i$  w.r.t.\   $\db$ invokes
$\tau_j$ directly or indirectly.  }
%
imposed by the acyclic support condition and access to the
table atoms  will be much simplified  if they are linearly
ordered.  We  therefore  topologically sort   $\taudb^{t}$
respecting the said partial order  and call the linearized
$\taudb^{t}$ satisfying the three assumptions (the acyclic
support condition, the  t-exclusiveness condition and  the
independent  condition)   a  {\em hierarchical  system  of
t-explanations for $G_t$}. We write it as $\taudb^{t} =
\langle      \tau_0^{t}, \tau_1^{t},\ldots, \tau_{K_t}^{t}
\rangle$ $(\tau_0 = G_t)$ assuming $\psidbtilde(\cdot)$ is
implicitly given.\footnote{
%
So now  it holds that  if $\tau_i$ precedes  $\tau_j$ then
$i<j$.}
%
Once a hierarchical  system of t-explanations for $G_t$ is
successfully built  from the source  program, equations on
inside  probability    and  outside  probability  such  as
(\ref{eq:inside-example})  and  (\ref{eq:outside-example})
are automatically derived and  solved in time proportional
to the size of the equations.   It plays a central role in
our approach to efficient EM learning.\\

One  way  to obtain  such t-explanations  is   to use OLDT
search~\cite{Tamaki86,Warren92},  a  complete   refutation
method  for logic programs.  In  OLDT search,  when a goal
$G$ is called for the first  time, we set  up an entry for
$G$  in a {\em   solution  table\/} and store  its  answer
substitutions $G\theta$ there.  When a call to an instance
$G'$ of $G$ occurs later, we stop solving $G'$ and instead
try to retrieve an answer substitution $G\theta$ stored in
the solution  table by unifying $G'$  with  $G\theta$.  To
record the  remaining   answer substitutions of    $G$, we
prepare a {\em lookup table\/} for $G'$ and hold a pointer
to them.

For self-containedness, we look  at details of OLDT search
using  a   sample  program  $\db_h   =   F_h\cup  R_h$  in
Figure~\ref{program:HMM}\footnote{
%
{\tt f1}, {\tt f2}, {\tt f3}, {\tt  h1}, {\tt h2} and {\tt
h3} are temporary marks, not part of the program.  }
%
which depicts an HMM\footnote{
%
An HMM defines a  probability distribution over strings in
the given  set of alphabets,  and  works  as a  stochastic
string     generator~\cite{Rabiner93}  such that an output
string is a sample from the defined distribution.  }
%
in Figure~\ref{fig:hmm1}.  This  HMM has two states  ${\tt
\{s0, s1\}}$.
%
\begin{figure}[hb]
\centerline{
	\includegraphics[scale=0.8]{./hmm1.eps}
}
\caption{Two state HMM}
\label{fig:hmm1}
%\end{figure}
%
%\begin{figure}[htb]
\[
\begin{array}{lcl}
F_h & = & \left\{
      \begin{small}
      \begin{tabular}{l}
      \verb|f1: values(init, [s0,s1]). | \\
      \verb|f2: values(out(_),[a,b]).  | \\
      \verb|f3: values(tr(_), [s0,s1]).| \\
      \end{tabular}
      \end{small}
      \right. \\
\\
R_h & = & \left\{
\begin{tabular}{ll}
\verb!h1: hmm(Cs):-               % To generate a string (chars) Cs,  ! \\
\verb!      msw(init,once,Si),    % Set initial state to Si, and then ! \\
\verb!      hmm(1,Si,Cs).         % Enter the loop with clock = 1.    ! \\
\verb!h2: hmm(T,S,[C|Cs]):- T=<3, % Loop:                             ! \\
\verb!      msw(out(S),T,C),      % Output C in state S.              ! \\
\verb!      msw(tr(S),T,NextS),   % Transit from S to NextS.          ! \\
\verb!      T1 is T+1,            % Put the clock ahead.              ! \\
\verb!      hmm(T1,NextS,Cs).     % Continue the loop (recursion).    ! \\
\verb!h3: hmm(T,_,[]):- T>3.      % Finish the loop if clock > 3.     ! \\
\end{tabular}
\right.
\end{array}
\]
\caption{Two state HMM program $\db_h$}
\label{program:HMM}
\end{figure}
%
At  a state  transition,  it probabilistically chooses the
next destination from  ${\tt  \{s0,  s1\}}$ and also    an
alphabet from   ${\tt \{a, b\}}$  to emit.   Note  that to
specify a fact  set $F_h$ and the associated  distribution
compactly, we introduce here   a new  notation  $\mbox{\tt
values($i$,[$v_1$,\ldots,$v_m$])}$. It declares that $F_h$
contains   {\tt  msw}   atoms   of  the  form   $\mbox{\tt
msw($i$,$n$,$v$)}$     ($v  \in \{v_1,\ldots,v_m\}$) whose
distribution   is  $P_{(i,n)}$  given by (\ref{eq:PF})  in
Subsection~\ref{subsec:simplifying-conditions}.        For
example, ({\tt f3}),   $\mbox{\tt values(tr(\_),[s0,s1])}$
introduces {\tt   msw(tr($t$),$n$,$v$)}  atoms   into  the
program  such  that $t$  can be  any  ground  term, $v \in
\{{\tt s0, s1}\}$ and for  a ground term  $n$, they have a
distribution
%
\[
P_{(\mbox{\footnotesize\tt tr($t$)},n)}
   (\mbox{\tt msw(tr($t$),$n$,s0)}=x,
    \mbox{\tt msw(tr($t$),$n$,s1)}=y
        \mid \theta_{i,s0}, \theta_{i,s1})
 = \theta^x_{i,s0}\theta^y_{i,s1}
\]
%
where $i = \mbox{\tt tr($t$)}$, $x,y \in \{0,1\}$ and $x+y
=1$.

This program runs like a Prolog program.  For a non-ground
top-goal  $\leftarrow  {\tt hmm(S)}$,   it functions as  a
stochastic string generator returning  a list of alphabets
such    as  {\tt [a,b,a]}  in   the   variable  {\tt S} as
follows. The  top-goal calls  clause ({\tt h1})  and ({\tt
h1}) selects the  initial state by executing  subgoal {\tt
msw(init,once,Si)}\footnote{
%
If {\tt  msw($i$,$n$,V)}   is called with ground  $i$  and
ground $n$, {\tt V},  a  logical variable, behaves  like a
random variable.  It is instantiated to some term $v$ with
probability $\theta_{i,v}$  selected  from the   value set
$V_i$ declared by a  {\tt values} atom.   If, on the other
hand,   {\tt V} is   a  ground term  $v$ when  called, the
procedural  semantics of  $\mbox{\tt msw($i$,$n$,$v$)}$ is
equal to that of $\mbox{\tt msw($i$,$n$,V)} \wedge {\tt V}
= v$.  }
%
which  returns in     ${\tt   Si}$   an initial      state
probabilistically chosen from \{{\tt s0}, {\tt s1}\}.  The
second clause ({\tt h2})  is called  from ({\tt h1})  with
ground  {\tt   S}  and  ground   {\tt  T}.   It  makes   a
probabilistic choice of   an output  alphabet {\tt  C}  by
asking {\tt   msw(out(S),T,C)}  and  then  determines {\tt
NextS},    the    next     state,     by     asking   {\tt
msw(tr(S),T,NextS)}.   ({\tt h3})  is  there to   stop the
transition.  For simplicity, the  length of output strings
is fixed to  three.  This  way of  execution  is termed as
{\em  sampling execution\/} because   it corresponds to  a
random  sampling from  $P_{\db_h}$.    If the  top-goal is
ground like   $\leftarrow\!\mbox{\tt   hmm([a,b,a])}$,  it
works as an acceptor,   i.e.\ returning success  (yes)  or
failure (no).\\

\begin{figure}[b]
\begin{verbatim}
 t1:  top_hmm(Cs,Ans):- tab_hmm(Cs,Ans,[]).
 t2:  tab_hmm(Cs,[hmm(Cs)|X],X):- hmm(Cs,_,[]).
 t3:  tab_hmm(T,S,Cs,[hmm(T,S,Cs)|X],X):- hmm(T,S,Cs,_,[]).
 t4:  e_msw(init,T,s0,[msw(init,T,s0)|X],X).
 t4': e_msw(init,T,s1,[msw(init,T,s1)|X],X).
   :
 t7:  hmm(Cs,X0,X1):- e_msw(init,once,Si,X0,X2), tab_hmm(1,Si,Cs,X2,X1).
 t8:  hmm(T,S,[C|Cs],X0,X1):-
          T=<3, e_msw(out(S),T,C,X0,X2), e_msw(tr(S),T,NextS,X2,X3),
          T1 is T+1, tab_hmm(T1,NextS,Cs,X3,X1).
 t9:  hmm(T,S,[],X,X):- T>3.
\end{verbatim}
\caption{Translated program of $\db_h$}
\label{program:translated-HMM}
\end{figure}

If  all explanations   for  $\mbox{\tt hmm([a,b,a])}$  are
sought for,  we keep  all  {\tt  msw} atoms  resolved upon
during the refutation as  a conjunction (explanation), and
repeat   this  process  by  backtracking   until  no  more
refutation is found.   If we need t-explanations  however,
backtracking must be abandoned  because sharing of partial
explanations  through     t-explanations,  the  purpose of
t-explanations  itself,  becomes impossible.  We therefore
instead use OLDT search for  all t-explanation search.  In
the   case  of the  HMM  program for  example,  to build a
hierarchical  system   of  t-explanations   for $\mbox{\tt
hmm([a,b,a])}$  by   OLDT search, we   first declare ${\tt
hmm}/1$ and ${\tt hmm}/3$ as table predicate.\footnote{
%
In general, ${\tt p}/n$  means  a predicate {\tt p}   with
arity $n$.  So  although ${\tt  hmm}/1$ and ${\tt  hmm}/3$
share   the predicate name  {\tt hmm},  they are different
predicates.}
%
So a t-explanation will be a  conjunction of ${\tt hmm}/1$
atoms, ${\tt hmm}/3$ atoms and  {\tt msw} atoms.  We  then
translate the  program    into  another  logic    program,
analogously to  the translation   of {\em definite  clause
grammars} (DCGs)  in Prolog~\cite{Sterling86}.  We add two
arguments  (which forms a  {\em D-list}) to each predicate
for the purpose of accumulating {\tt  msw} atoms and table
atoms as  conjuncts  in a  t-explanation.  The translation
applied      to   $\db_h$     yields    the    program  in
Figure~\ref{program:translated-HMM}.

In the translated program,  clause ({\tt t1})  corresponds
to  the top-goal $\leftarrow\!\mbox{\tt hmm($l$)}$ with an
input  string $l$, and a  t-explanation for the table atom
{\tt hmm($l$)} will be returned  in {\tt Ans}.  ({\tt t2})
and ({\tt  t3})   are auxiliary  clauses   to add  to  the
callee's  D-list a table atom of  the  form {\tt hmm($l$)}
and  {\tt  hmm($t$,$s$,$l$)} respectively ($t$: time step,
$s$: state).   In  general,  if ${\tt  p}/n$  is  a  table
predicate in the original program, ${\tt p}/(n+2)$ becomes
a  table   predicate  in the  translated    program and an
auxiliary predicate ${\tt   tab\_p}/(n+2)$ is inserted  to
signal the OLDT interpreter  to  check the solution  table
for ${\tt  p}/n$, i.e.\ to  check  if there  already exist
t-explanations   for ${\tt p}/n$.   Likewise clauses ({\tt
t4}) and ({\tt  t4'})  are a pair  corresponding  to ({\tt
f1}) which   insert  {\tt msw(init,T,$\cdot$)}   to    the
callee's D-list with  ${\tt T}={\tt once}$.   Clauses {\tt
(t7)}, {\tt (t8)}  and {\tt  (t9)} respectively correspond
to {\tt (h1)}, {\tt (h2)} and {\tt (h3)}.

\begin{figure}[bht]
\centerline{
	\includegraphics[scale=0.95]{./soltab.eps}
}
\caption{Solution table}
\label{tab:solution-tab}
\end{figure}

Then  after translation,    we   apply  OLDT  search    to
$\leftarrow\!    \mbox{\tt  top\_hmm([a,b,a],Ans)}$  while
noting  (i) the added D-list  does not  influence the OLDT
procedure,  and (ii) we  associate with each solution of a
table atom in the solution table a list of t-explanations.
The    resulting   solution     table   is     shown    in
Figure~\ref{tab:solution-tab}.  The first row reads that a
call to   {\tt   hmm([a,b,a])} occurred and   entered  the
solution table  and its  solution, {\tt hmm([a,b,a])}  (no
variable  binding  generated),  has  two   t-explanations,
$\mbox{\tt    msw(init,once,s0)}      \wedge     \mbox{\tt
hmm(1,s0,[a,b,a])}$     and $\mbox{\tt  msw(init,once,s1)}
\wedge  \mbox{\tt hmm(1,s1,[a,b,a])}$.  The remaining task
is  the topological sorting  of  the table atoms stored in
the solution   table     respecting the acyclic    support
condition.  This can  be done by using  depth-first search
(trace) of t-explanations from  the top-goal  for example.
Thus we obtain a hierarchical system of t-explanations for
$\mbox{\tt hmm([a,b,a])}$.

\subsection{Support Graphs}

Looking back, all we need  to  compute inside and  outside
probability is a   hierarchical system of  t-explanations,
which essentially  is a  boolean combination  of primitive
events ({\tt msw} atoms) and compound events (table atoms)
and  as such can   be more intuitively representable  as a
graph.  For  this  reason,  and  to help   visualizing our
learning  algorithm, we   introduce  a new  data-structure
termed {\em support graph\/}s, though the new EM algorithm
in  the next subsection itself  is described solely by the
hierarchical system of t-explanations.

\begin{figure}[htb]
\centerline{
	\includegraphics[scale=0.90]{./sgraph.eps}
}
\caption{A support graph (a) in general form, (b) for
    $G_t=\mbox{\tt hmm([a,b,a])}$ in the HMM program $\db_h$.
	A double-circled node refers to a table node.
	}
\label{fig:sgraph}
\end{figure}

As illustrated in Figure~\ref{fig:sgraph} (a), the support
graph for  $G_t$ is    a graphical
representation of       the   hierarchical   system     of
t-explanations   $\taudb^{t}       =   \langle \tau_0^{t},
\tau_1^{t},\ldots, \tau_{K_t}^{t} \rangle$  $(\tau_0^{t} =
G_t)$ for    $G_t$ in   (\ref{eq:tabled-completion}).   It
consists of  totally ordered disconnected  subgraphs, each
of which is  labeled   with the corresponding table   atom
$\tau_k^{t}$   in $\taudb^{t}$ ($0\leq   k  \leq K_t$).  A
subgraph labeled $\tau_k^{t}$  comprises two special nodes
(the {\em  start node\/} and  the {\em end node}) and {\em
explanation      graph\/}s,  each corresponding    to    a
t-explanation            $\Stil_{k,j}^{t}$              in
$\psidbtilde(\tau_k^{t})$ ($1 \leq j
\leq m_k$).

An  {\em  explanation graph\/}  of $\Stil_{k,j}^{t}$  is a
linear   graph in which a node   is labeled  either with a
table  atom   $\tau$ or     with   a  switch    $\mbox{\tt
msw($\cdot$,$\cdot$,$\cdot$)}$ in $\Stil_{k,j}^{t}$.  They
are called a  {\em table node\/} and a {\em switch node\/}
respectively.  Figure~\ref{fig:sgraph}  (b) is the support
graph  for $\mbox{\tt   hmm([a,b,a])}$  obtained  from the
solution   table   in Figure~\ref{tab:solution-tab}.  Each
table node labeled $\tau$  refers to the  subgraph labeled
$\tau$, so data-sharing  is achieved through the  distinct
table nodes referring to the same subgraph.


\subsection{Graphical EM Algorithm}
\label{sec:learning:GEM}

We describe here an efficient EM learning algorithm termed
the {\em  graphical EM  algorithm}  (Figure~\ref{alg:GEM})
introduced   by   \citeA{Kameya00}, that  runs on  support
graphs.  Suppose we have a random  sample ${\cal G} = G_1,
\ldots, G_T$  of observable  atoms.  Also suppose  support
graphs  for $G_t$ $(1\leq  t  \leq T)$, i.e.\ hierarchical
systems of  t-explanations  satisfying the acyclic support
condition,  the  t-exclusiveness     condition   and   the
independent  condition, have been successfully constructed
from a  parameterized logic  program $\db$  satisfying the
uniqueness condition and the finite support condition.

The    graphical   EM      algorithm               refines
\pname{learn-naive($\db$,${\cal G}$)} by introducing   two
subroutines, \pname{get-inside-probs($\db$, ${\cal G}$)\/}
to         compute  inside       probabilities         and
\pname{get-expectations($\db$,  ${\cal G}$)\/} to  compute
outside   probabilities.   They are  called  from the main
routine      \pname{learn-gEM($\db$,${\cal   G}$)}.   When
learning, we prepare four arrays   for each support  graph
for $G_t$ in ${\cal G}$:

\begin{itemize}
\item
${\cal P}[t,\tau]$ for  the {\em inside probability\/}  of
$\tau$, i.e.\ $\beta(\tau)  =  \pdb(\tau \mid  \btheta)$
(see (\ref{eq:inside-def}))

\item 
${\cal Q}[t,\tau]$ for  the {\em outside probability\/} of
$\tau$    w.r.t.\ $G_t$,    i.e.\  $\alpha(G_t,\tau)$ (see
(\ref{eq:outside-ini}) and (\ref{eq:outside-rec}))

\item
${\cal     R}[t,\tau,\Stil]$  for    the  {\em explanation
probability\/}  of $\Stil\, (\in\psidbtilde(\tau_k^{t}))$,
i.e.\ $\pdb(\Stil \mid \btheta)$

\item
$\eta[t,i,v]$ for the {\em expected count\/} of $\mbox{\tt
msw($i$,$\cdot$,$v$)}$,       i.e.\       $\sum_{S     \in
\psi_{\db}(G_t)}    P_{\tt    msw}(S   \mid       \btheta)
\sigma_{i,v}(S)$

\end{itemize}
%
and  call the   procedure    \pname{learn-gEM($\db$,${\cal
G}$)\/} in Figure~\ref{alg:GEM}.
%
\begin{figure}[tb]
\samepage
\begin{tabular}{|c|c|}
 	\hline
\begin{minipage}[t]{.46\textwidth}
\begin{listing}
\vspace{0.05em}
\item \smbf{procedure} {\em learn-gEM\/}($\db,{\cal G}$)
\item \smbf{begin}
\item \hspace{0.7em}
         Select some $\btheta$ as initial \\
         \hspace*{2em}parameters;
%		\hq\progcomment{Initialize the parameters appropriately.}
\item \hspace{0.7em}
        {\em get-inside-probs\/}($\db,{\cal G}$);
	\label{linenum:learn-gEM:get-init-iniside}
%		\hq\progcomment{Compute inside and local probs.}
\item \hspace{0.7em}
        $\lambda^{(0)} :=\, \sum_{t=1}^T\ln{\cal P}[t,G_t]$;
%		\hq\progcomment{Compute the log-likelihood of $\Gseq$.}
	\label{linenum:learn-gEM:get-init-likelihood}
\item \hspace{0.7em}
        \smbf{repeat}
\item \hspace{1.4em}
          {\em get-expectations\/}($\db,{\cal G}$);
%		\hq\progcomment{Expected occurence of
%						$\sym{msw($i$,$\cdot$,$v$)}$.}
	\label{lnum:learn-gEM:after-on-tiv}
\item \hspace{1.4em}
       \smbf{foreach} $i\in I,v\in V_i$  \smbf{do}
\item \hspace{2.1em} $\eta[i,v] \,:=$\\
      \hspace*{4em}
                  $\sum_{t=1}^T \eta[t,i,v]/{\cal P}[t,G_t]$;
	\label{list:learn-gEM:get-On}
\item \hspace{1.4em}
      \smbf{foreach} $i\in I,v\in V_i$  \smbf{do}
\item \hspace{2.1em}
         $\theta_{i,v} :=\, \eta[i,v]/\sum_{v'\in V_i}\eta[i,v']$;
	\label{lnum:learn-gEM:update}
%\halfempline
\item \hspace{1.4em}
         {\em get-inside-probs\/}$(\db,{\cal G})$;
	\label{lnum:learn-gEM:inside-before-loop}
\item \hspace{1.4em}
         $m := m+1$;
\item \hspace{1.4em}
         $\lambda^{(m)} :=\, \sum_{t=1}^T\ln {\cal P}[t,G_t]$
%		\hq\progcomment{Recompute the log-likelihood of $\Gseq$.}
%\halfempline
\item \hspace{0.7em}
        \smbf{until} $\lambda^{(m)}-\lambda^{(m-1)}<\varepsilon$
%		\hq\progcomment{Terminate if the loglikelihood converged.}
\item \smbf{end}.
\end{listing}
\end{minipage}&

\begin{minipage}[t]{.48\textwidth}
%%%%%%%% get_inside_prob %%%%%%%%
\begin{listing}
\vspace{0.05em}
\item \smbf{procedure} {\em get-inside-probs\/}$(\db,{\cal G})$
\item \smbf{begin}
\item \hspace{0.7em}
          \smbf{for} $t := 1$ \smbf{to} $T$ \smbf{do} \smbf{begin}
\item \hspace{1.4em} Let $\tau_0^{t} = G_t$;
    \label{line:get-inside-probs:top-layer}
\item \hspace{1.4em} \smbf{for} $k := K_t$  \smbf{downto} $0$
                     \smbf{do} \smbf{begin}
    \label{line:get-inside-probs:bottom-layer}
\item \hspace{2.1em} ${\cal P}[t,\tau_k^{t}] := 0$;
\item \hspace{2.1em}
         \smbf{foreach}
          $\Stil\in\psidbtilde(\tau_k^{t})$  \smbf{do} \smbf{begin}
    \label{line:get-inside-probs:Stil}
\item \hspace{2.8em}
          Let $\Stil=\{A_1,A_2,\ldots,A_{|\Stil|}\}$;
\item \hspace{2.8em}
          ${\cal R}[t,\tau_k^{t},\Stil] := 1$;
\item \hspace{2.8em}
           \smbf{for} $l := 1$  \smbf{to} $|\Stil|$  \smbf{do}
\item  \hspace{3.5em}
            \smbf{if} $A_l=\mbox{\tt msw($i$,$\cdot$,$v$)}$
			\smbf{then}
\item \hspace{4.2em}  ${\cal R}[t,\tau_k^{t},\Stil] \,\ast\!=\, \theta_{i,v}$
          \label{line:get-inside-probs:theta}
\item \hspace{3.5em}  \smbf{else}
            ${\cal R}[t,\tau_k^{t},\Stil] \,\ast\!=\, {\cal P}[t,A_l]$;
          \label{line:get-inside-probs:tau}
\item \hspace{2.8em}
	${\cal P}[t,\tau_k^{t}] \,+\!\!=\, {\cal R}[t,\tau_k^{t},\Stil]$
\item \hspace{2.1em} \smbf{end} {/* \smbf{foreach} $\Stil$ */}
\item \hspace{1.4em} \smbf{end} {/* \smbf{for} $k$ */}
\item \hspace{0.7em} \smbf{end} {/* \smbf{for} $t$ */}
\item \smbf{end}.
%\halfempline
\end{listing}
\vspace{0.1em}
	\end{minipage}\\
	\hline
	\multicolumn{2}{|c|}{
\begin{minipage}[t]{0.9\textwidth}
%%%%%%%% get_On %%%%%%%%%
\begin{listing}
\vspace{0.05em}
\item \smbf{procedure} {\em get-expectations\/}$(\db, {\cal G})$
      \smbf{begin}
\item \hspace{0.7em}
         \smbf{for} $t := 1$  \smbf{to} $T$  \smbf{do} \smbf{ begin}
\item \hspace{1.4em}
         \smbf{foreach} $i\in I,v\in V_i$  \smbf{do}
         $\eta[t,i,v] := 0$;
\item \hspace{1.4em}
         Let $\tau_0^{t} = G_t$;
         ${\cal Q}[t,\tau_0^{t}] \,:=1.0$;
    \label{line:get-expectations:ini-outside}
\item \hspace{1.4em}
         \smbf{for} $k := 1$  \smbf{to} $K_t$
         \smbf{do} ${\cal Q}[t,\tau_k^{t}] \,:=0$;
\item \hspace{1.4em}
         \smbf{for} $k := 0$  \smbf{to} $K_t$  \smbf{do}
    \label{line:get-expectations:loop-outside}
\item \hspace{2.1em}
         \smbf{foreach} $\Stil\in\psidbtilde(\tau_k^{t})$  \smbf{ do}
         \smbf{begin}
\item \hspace{2.8em} Let $\Stil=\{A_1,A_2,\ldots,A_{|\Stil|}\}$;
\item \hspace{2.8em} \smbf{for} $l := 1$  \smbf{to} $|\Stil|$  \smbf{do}
\item \hspace{3.5em}
             \smbf{if} $A_l=\mbox{\tt msw($i$,$\cdot$,$v$)}$
				\smbf{then}
                $\eta[t,i,v] \, +\!\!=\,
                 {\cal Q}[t,\tau_k^{t}] \cdot {\cal R}[t,\tau_k^{t},\Stil]$
      \label{line:get-expectations:counts}
\item \hspace{3.5em}
            \smbf{else}
                ${\cal Q}[t,A_l] \, +\!\!=\,
                    {\cal Q}[t,\tau_k^{t}]\cdot
                    {\cal R}[t,\tau_k^{t},\Stil]/{\cal P}[t,A_l]$
      \label{line:get-expectations:tau-outside}
\item \hspace{2.1em} \smbf{end} {/* \smbf{ foreach} $\Stil$} */
\item \hspace{0.7em} \smbf{end} {/* \smbf{ for} $t$ */}
\item \smbf{end}.
%\halfempline
\end{listing}
\vspace{0.1em}
\end{minipage}
  } \\ \hline
\end{tabular}

\caption{graphical EM algorithm.}
\label{alg:GEM}
\end{figure}
%
The   main  routine  \pname{learn-gEM($\db$,${\cal G}$)\/}
initially    computes all     inside         probabilities
(Line~\ref{linenum:learn-gEM:get-init-iniside}) and enters
a loop in which \pname{get-expectations($\db$,${\cal G}$)}
is   called  first   to   compute    the expected    count
$\eta[t,i,v]$  of  $\mbox{\tt  msw($i$,$\cdot$,$v$)}$  and
parameters are updated (Line~\ref{lnum:learn-gEM:update}).
Inside probabilities   are renewed  by using  the  updated
parameters   before   entering   the     next         loop
(Line~\ref{lnum:learn-gEM:inside-before-loop}).

The   subroutine      \pname{get-inside-probs($\db$,${\cal
G}$)\/} computes   the inside probability  $\beta(\tau)  =
\pdb(\tau    \mid  \btheta)$  (and  stores it  in   ${\cal
P}[t,\tau]$) of a table atom $\tau$  from the bottom layer
up     to    the  topmost    layer     $\tau_0   =    G_t$
(Line~\ref{line:get-inside-probs:top-layer})    of     the
hierarchical  system   of  t-explanations for   $G_t$ (see
(\ref{eq:tabled-completion})                            in
Subsection~\ref{subsec:io-prob}).     It       takes     a
t-explanation $\Stil$  in $\psidbtilde(\tau_k^{t})$ one by
one  (Line~\ref{line:get-inside-probs:Stil}),   decomposes
$\Stil$  into  conjuncts   and  multiplies    their inside
probabilities         which       are     either     known
(Line~\ref{line:get-inside-probs:theta})     or    already
computed (Line~\ref{line:get-inside-probs:tau}).

The other subroutine  \pname{get-expectations($\db$,${\cal
G}$)\/} computes   outside   probabilities  following  the
recursive      definitions    (\ref{eq:outside-ini})   and
(\ref{eq:outside-rec})  in Subsection~\ref{subsec:io-prob}
and stores the outside probability $\alpha(G_t,\tau)$ of a
table atom $\tau$  in  ${\cal Q}[t,\tau]$.  It first  sets
the outside probability of the top-goal  $\tau_0 = G_t$ to
$1.0$  (Line~\ref{line:get-expectations:ini-outside})  and
computes  the     rest  of        outside    probabilities
(Line~\ref{line:get-expectations:loop-outside}) going down
the  layers of the  t-explanation  for $G_t$  described by
(\ref{eq:tabled-completion})                            in
Subsection~\ref{subsec:io-prob}.
(Line~\ref{line:get-expectations:counts})    adds   ${\cal
Q}[t,\tau_k^{t}]\cdot     {\cal  R}[t,\tau_k^{t},\Stil] =$
$\alpha(G_t,\tau_k^{t})     \cdot    \beta(\Stil)$      to
$\eta[t,i,v]$,        the   expected    count    of   {\tt
msw($i$,$\cdot$,$v$)},   as    a   contribution  of   {\tt
msw($i$,$\cdot$,$v$)}  in $\Stil$ through $\tau_k^{t}$  to
$\eta[t,i,v]$.
(Line~\ref{line:get-expectations:tau-outside})  increments
the outside probability ${\cal Q}[t,A_l] =\alpha(G_t,A_l)$
of $A_l$ according to the equation~(\ref{eq:outside-rec}).
Notice  that  ${\cal  Q}[t,\tau_k^{t}]$  has  already been
computed    and   ${\cal R}[t,\tau_k^{t},\Stil]    / {\cal
P}[t,A_l] =\beta(W)$ for $\Stil = A_l \wedge W$.  As shown
in                 Subsection~\ref{subsec:learning:naive},
\pname{learn-naive($\db$,${\cal     G}$)\/}     is the MLE
procedure, hence the following theorem holds.

\begin{theorem}
\label{th:gEM}
\samepage
Let $\db$ be a parameterized  logic program, and ${\cal G}
=  G_1,\ldots,G_T$ a  random sample   of observable atoms.
Suppose the   five conditions (uniqueness,  finite support
(Subsection~\ref{subsec:simplifying-conditions}),  acyclic
support,        t-exclusiveness     and       independence
(Subsection~\ref{subsec:OLDT}))    are      met.      Then
\pname{learn-gEM  $(\db,{\cal  G})$\/}     finds  the  MLE
$\btheta^{\ast}$ which  (locally) maximizes the likelihood
$L({\cal G} \mid \btheta)  = \prod_{t=1}^{T} \pdb(G_t \mid
\btheta)$.
\end{theorem}
%
\noindent
(Proof) Sketch.\footnote{
%
A formal proof is  given by \citeA{Kame00b}.  It is proved
there   that  under   the  common    parameters  $\theta$,
$\eta[i,v]$   in  \pname{learn-naive($\db$,${\cal  G}$)\/}
coincides       with              $\eta[i,v]$           in
\pname{learn-gEM($\db$,${\cal G}$)\/}.  So, the parameters
are updated to the same  values.  Hence, starting with the
same initial values,  the parameters converge to  the same
values.  }
%
Since    the  main  routine  \pname{learn-gEM($\db$,${\cal
G}$)\/} is  the same  as   \pname{learn-naive($\db$,${\cal
G}$)\/} except    the   computation   of   $\eta[i,v]    =
\sum_{t=1}^T \eta[t,i,v]$,  we  show  that $\eta[t,i,v]  =
\sum_{S \in \psi_{\db}(G_t)} P_{\tt msw}(S \mid \btheta)
\sigma_{i,v}(S)$ $( =\sum_n \sum_{\mbox{\footnotesize\tt
msw($i$,$n$,$v$)} \in S \in \psi_{\db}(G_t)} P_{\tt msw}(S
\mid \btheta))$. However,

\begin{eqnarray*}
\eta[t,i,v]
 & = &
  \sum_{0\leq k\leq K_t} \sum_n
			\sum_{\mbox{\footnotesize\tt msw($i$,$n$,$v$)}
              \in  \Stil \in \psidbtilde(\tau_k^{t})}
        \alpha(G_t,\tau_k^{t}) \beta(\Stil)\\
 &   &
        \mbox{(see (Line~\ref{line:get-expectations:counts})
               in \pname{get-expectations}$(\db, {\cal G})$)}  \\
 & = & \sum_n \alpha(G_t,\mbox{\tt msw($i$,$n$,$v$)})
              \beta(\mbox{\tt msw($i$,$n$,$v$)})\\
 & = & \sum_n \mu(G_t,\mbox{\tt msw($i$,$n$,$v$)}) \;\;
          \mbox{(see the equation~(\ref{eq:outside-m}))} \\
 & = & \sum_n \sum_{\mbox{\footnotesize\tt   msw($i$,$n$,$v$)}
				\in S \in \psi_{\db}(G_t)}
          P_{\tt msw}(S \mid \btheta). \hspace{3em}\mbox{Q.E.D.}
\end{eqnarray*}

Here we used the fact that  if $\Stil$ contains $\mbox{\tt
msw($i$,$n$,$v$)}$  like $\Stil  = \Stil'\wedge  \mbox{\tt
msw($i$,$n$,$v$)}$,               $\beta(\Stil)          =
\beta(\Stil') \beta(\mbox{\tt   msw($i$,$n$,$v$)})$  holds,
and hence
%
\begin{eqnarray*}
\samepage
\lefteqn{\alpha(G_t,\tau_k^{t})\beta(\Stil)
  =  \alpha(G_t,\tau_k^{t})
		\beta(\Stil')\beta(\mbox{\tt msw($i$,$n$,$v$)})} \\
 & = & (\mbox{contribution of \mbox{\tt msw($i$,$n$,$v$)}
              in $\Stil$ through $\tau_k^{t}$ to}\;
			\alpha(G_t,\mbox{\tt msw($i$,$n$,$v$)}))
			\beta(\mbox{\tt msw($i$,$n$,$v$)}).
\end{eqnarray*}

The five conditions on the  applicability of the graphical
EM  algorithm  may    look  hard  to   satisfy   at  once.
Fortunately,     the       modeling         principle   in
Section~\ref{subsec:modeling-principle}  still stands, and
with due care  in modeling, it is likely  to lead us  to a
program that meets all of  them.  Actually, we will see in
the    next    section,     programs       for    standard
symbolic-statistical frameworks such as Bayesian networks,
HMMs and PCFGs all satisfy the five conditions.


\section{Complexity}
\label{sec:complexity}

In this section,  we analyze  the  time complexity of  the
graphical     EM     algorithm     applied  to     various
symbolic-statistical   frameworks  including HMMs,  PCFGs,
pseudo PCSGs and Bayesian networks.  The results show that
the  graphical  EM  algorithm   is competitive  with these
specialized EM algorithms  developed independently in each
research field.

\subsection{Basic Property}

Since  the  EM algorithm   is an iterative  algorithm and
since we are   unable to predict  when  it converges,  we
measure  time complexity  by    the time taken  for   one
iteration.  We therefore  estimate time per iteration  on
the \smbf{repeat}   loop  of {\em learn-gEM\/}($\db,{\cal
G}$) (${\cal G} = G_1,\ldots, G_T$).   We observe that in
one iteration, each support graph for $G_t$ ($1\leq t
\leq  T$)    is     scanned   twice, once      by    {\em
get-inside-probs\/}($\db,{\cal    G}$)   and once by {\em
get-expectations\/}($\db,{\cal G}$).  In      the   scan,
addition   is  performed  on   the   t-explanations,  and
multiplication (possibly  with  division) is performed on
the  {\tt msw} atoms  and table atoms once  for each.  So
time  spent for $G_t$  per iteration  by the graphical EM
algorithm is  linear   in the  {\em  size  of the support
graph}, i.e.\  the number of nodes  in the  support graph
for $G_t$.  Put
%
\begin{eqnarray*}
\widetilde{\Delta}_{\db}^{t}
   & \defined &
       \bigcup_{\tau\in\taudb^{t}}  \psidbtilde(\tau)\\
\xi_{\rm num}
   & \defined  &
        \max_{1\leq t\leq T} |\widetilde{\Delta}_{\db}^{t}| \\
\xi_{\rm maxsize}
   & \defined &
       \max_{1\leq t\leq T, \Stil\in\widetilde{\Delta}_{\db}^{t}} |\Stil|.
\end{eqnarray*}
%
Recall that $\taudb^{t}$  is the set  of table  atoms for
$G_t$, and  hence $\widetilde{\Delta}_{\db}^{t}$  is  the
set of  all  t-explanations appearing  in  the right hand
side     of         (\ref{eq:tabled-completion})       in
Subsection~\ref{subsec:OLDT}.  So $\xi_{\rm num}$ is  the
maximum number of t-explanations  in a support  graph for
the $G_t$'s and $\xi_{\rm maxsize}$ the maximum size of a
t-explanation for   the      $G_t$'s respectively.    The
following is obvious.

\begin{proposition}
\label{propositon:complexity}
The time  complexity of  the  graphical EM  algorithm per
iteration is linear in the  total size of support graphs,
$O(\xi_{\rm num}\xi_{\rm  maxsize} T)$ in notation, which
coincides with the space complexity because the graphical
EM algorithm runs on support graphs.
\end{proposition}

This is a rather general result, but  when we compare the
graphical EM algorithm  with other EM algorithms, we must
remember that the input to  the graphical EM algorithm is
support graphs (one  for  each  observed atom) and    our
actual total learning time is

\begin{center}
OLDT time + (the number of iterations)
       $\times O(\xi_{\rm num}\xi_{\rm maxsize} T)$
\end{center}
%
where ``OLDT time'' denotes time to construct all support
graphs  for ${\cal G}$.  It is  the sum  of time for OLDT
search and time for the topological  sorting of the table
atoms, but because the latter    is part of the    former
order-wise,\footnote{
%
Think of OLDT  search for a  top-goal $G_t$.  It searches
for {\tt msw} atoms and table atoms  to create a solution
table,   while    doing some     auxiliary  computations.
Therefore  its  time  complexity  is   never  less   than
$O(|\mbox{the  number of {\tt msw}  atoms and table atoms
in the support graph for  $G_t$}|)$, which coincides with
the time we need to topologically sort table atoms in the
solution table by depth-first search from $\tau_0 = G_t$.
}
%
we represent ``OLDT time'' by time for OLDT search.  Also
observe  that the total size of   support graphs does not
exceed time for OLDT  search  for ${\cal  G}$ order-wise.

To evaluate OLDT time for a specific class of models such
as  HMMs,  we need to  know   time for table  operations.
Observe that our OLDT search  in this paper is special in
the sense that table atoms are  always ground when called
and there is no resolution with solved goals. Accordingly
a solution table is used only

\begin{itemize}
\item
to  check if a  goal  $G$  already has  an  entry  in the
solution table, i.e.\ if it was called before, and

\item
to add  a new searched t-explanation for  $G$ to the list
of discovered t-explanations under $G$'s entry.

\end{itemize}

The time complexity of  these operations is equal to that
of table access which depends  both on the program and on
the implementation of the solution table.\footnote{
%
\citeA{Sagonas94}  and  \citeA{Ramakrishnan95}   discuss
about the implementation of OLDT.}
%
We first suppose programs are carefully written in such a
way that the arguments of   table atoms used as  indecies
for table   access are  integers.  Actually  all programs
used in the  subsequent complexity analysis (${\db}_h$ in
Subsection~\ref{subsec:OLDT}, ${\db}_g$ and  ${\db}_{g'}$
in Subsection~\ref{subsec:PCFGs},  ${\db}_{G^{\tau}}$  in
Subsection~\ref{subsec:Bayes})  satisfy   or  can satisfy
this  condition  by   replacing  non-integer  terms  with
appropriate integers.  We  also suppose that the solution
table is implemented using an  array so that table access
can be done in $O(1)$ time.\footnote{
%
If arrays   are not  available,  we  may be able   to use
balanced trees, giving $O(\log  n)$ access time where $n$
is the number  data in the solution table,  or  we may be
able to  use hashing, giving   average $O(1)$ time access
under a certain condition \cite{Cormen90}.  }

In what follows, we  present  a detailed analysis of  the
time complexity of  the graphical EM algorithm applied to
HMMs, PCFGs, pseudo PCSGs and Bayesian networks, assuming
$O(1)$  time access to the  solution table.  We remark by
the way  that  their space complexity  is  just the total
size of solution tables (support graphs).

\subsection{HMMs}
\label{subsec:HMMs}

The   standard  EM algorithm for  HMMs   is the Baum-Welch
algorithm  \cite{Rabiner89,Rabiner93}.  An  example of HMM
is         shown  in     Figure~\ref{fig:hmm1}          in
Subsection~\ref{subsec:OLDT}.\footnote{
%
We  treat here only  ``state-emission HMMs'' which emit a
symbol   depending   on    the  state.   Another    type,
``arc-emission HMMs'' in which the emitted symbol depends
on the transition arc, is treated similarly.  }
%
Given $T$  observations $w^1,\ldots,w^T$ of output string
of  length $L$, it   computes in $O(N^2LT)$  time in each
iteration   the  forward probability  $\alpha_m^t(q)  = $
$P(o^t_1  o^t_2  \cdots o^t_{m-1},\,q \mid  \btheta)$ and
the    backward   probability $\beta_m^t(q)  =$  $P(o^t_m
o^t_{m+1} \cdots  o^t_L \mid q,\,\btheta)$ for each state
$q \in Q$, time step $m$ $(1 \leq m \leq L)$ and a string
$w^t =  o^t_1  o^t_2 \cdots o^t_L$  $(1  \leq t \leq T)$,
where $Q$ is  the  set of states  and  $N$ the  number of
states.  The factor $N^2$ comes  from the fact that every
state  has  $N$  possible  destinations and   we  have to
compute the  forward  and backward probability  for every
destination   and   every state.   After    computing all
$\alpha_m^t(q)$'s   and  $\beta_m^t(q)$'s, parameters are
updated.  So,    the  total   computation time   in  each
iteration of  the  Baum-Welch  algorithm is estimated  as
$O(N^2LT)$ \cite{Rabiner93,Manning99}.

To compare this result with the graphical EM algorithm, we
use  the HMM program  $\db_h$ in  Figure~\ref{program:HMM}
with appropriate modifications   to $L$, the length  of  a
string, $Q$, the state set, and  declarations in $F_h$ for
the output  alphabets.  For a string  $w  = o_1 o_2 \cdots
o_L$, $\mbox{\tt hmm($n$,$q$,[$o_m,o_{m+1},\ldots,o_L$])}$
in $\db_h$ reads that  the HMM is  in state  $q \in Q$  at
time    $n$     and     has     to   output     $\mbox{\tt
[$o_m$,$o_{m+1}$,\ldots,$o_L$]}$   until it reaches    the
final  state.   After declaring   ${\tt hmm}/1$  and ${\tt
hmm}/3$ as  table     predicate   and translation     (see
Figure~\ref{program:translated-HMM}), we apply OLDT search
to       the            goal    $\leftarrow\!    \mbox{\tt
top\_hmm([$o_1$,\ldots,$o_L$],Ans)}$     w.r.t.\       the
translated program to  obtain all t-explanations  for {\tt
hmm([$o_1$,\ldots,$o_L$])}.  For a   complexity   argument
however, the translated program  and $\db_h$ are the same,
so we talk in terms of $\db_h$ for the sake of simplicity.
In   the   search, we fix    the search  strategy  to {\em
multi-stage  depth-first   strategy}~\cite{Tamaki86}.   We
assume   that the solution  table  is accessible in $O(1)$
time.\footnote{
%
$O(1)$ is possible   because in   the translated  program
$\db_h$ in  Section~\ref{subsec:OLDT},  we can identify a
goal               pattern          of               {\tt
hmm($\cdot$,$\cdot$,$\cdot$,$\cdot$,$\cdot$)}    by   the
first two arguments which are constants (integers).}
%
Since the  length of  the  list in the third  argument of
${\tt  hmm}/3$ decreases by  one  on each recursion,  and
there  are  only  finitely many    choices of  the  state
transition   and     the   output alphabet,   the  search
terminates, leaving finitely  many t-explanations in  the
solution  table like  Figure~\ref{tab:solution-tab}  that
satisfy the acyclic support condition respectively.  Also
the sampling execution of $\leftarrow\!\mbox{\tt hmm(L)}$
w.r.t.\ $\db_h$  is nothing  but  a sequential   decision
process  such that decisions made by  {\tt msw} atoms are
exclusive,  independent   and  generate a  unique string,
which  means   $\db_h$  satisfies   the   t-exclusiveness
condition,  the independence condition and the uniqueness
condition respectively.   So, the graphical  EM algorithm
is applicable  to   the set  of  hierarchical systems  of
t-explanations for {\tt hmm($w^t$)}   $(1\leq t \leq  T)$
produced   by   OLDT  search   for   $T$  observations  $
w^1,\ldots,w^T$ of     output    string.  Put    $w^t   =
o_1^to_2^t\cdots o_L^t$. It follows from

\[
\tau_{{\db}_h}^t 
    = 
    \{
      \mbox{\tt hmm($m$,$q$,[$o_m^t$,\ldots,$o_L^t$])}
		\mid 1 \leq m \leq L+1, q\in Q \}
     \cup
           \{ \mbox{\tt hmm([$o_1^t$,\ldots,$o_L^t$])} \}
\]
\[
\begin{array}[t]{c}
\widetilde{\psi}_{{\db}_h}(\mbox{\tt hmm($m$,$q$,[$o_m^t$,\ldots,$o_L^t$])}) \\
(1 \leq m \leq L)
\end{array}
= 
  \left\{  \left.
      \begin{array}{cc}
        \mbox{{\tt msw(out($q$),$m$,$o_m$)}},
        \mbox{{\tt msw(tr($q$),$m$,$q'$)}},\\
        \mbox{{\tt hmm($m+1$,$q'$,[$o_{m+1}^t$,\ldots,$o_L^t$])}}
      \end{array}
         \right|
        q' \in Q
 \right\}
\]
%
that  for  a top-goal {\tt hmm([$o_1^t$,\ldots,$o_L^t$])},
there  are   at  most $O(NL)$  calling  patterns  of ${\tt
hmm}/3$ and  each call causes  at most  $N$ calls to ${\tt
hmm}/3$, implying  there occur  $O(NL\cdot N)   = O(N^2L)$
calls to ${\tt hmm}/3$.   Since each call is computed once
due to  the  tabling mechanism, we  have  $\xi_{\rm num} =
O(N^2L)$.   Also     $\xi_{\rm maxsize}     =3$.  Applying
Proposition~\ref{propositon:complexity}, we reach

\begin{proposition}
\label{proposition:HMM}
Suppose we have $T$  strings of length $L$.  Also suppose
each table   operation in OLDT  search  is done in $O(1)$
time.  OLDT   time by   $\db_h$  is  $O(N^2LT)$   and the
graphical   EM algorithm    takes   $O(N^2LT)$  time  per
iteration where $N$ is the number of states.
\end{proposition}

$O(N^2LT)$ is   the  time complexity   of  the  Baum-Welch
algorithm.    So the   graphical   EM  algorithm  runs  as
efficiently as the Baum-Welch algorithm.\footnote{
%
Besides, the  Baum-Welch  algorithm and the  graphical EM
algorithm whose  input  are support  graphs generated  by
$\db_h$  update parameters to  the  same value if initial
values are the same.  }
\\

By     the     way,      the      Viterbi       algorithm
\cite{Rabiner89,Rabiner93} provides for HMMs an efficient
way of  finding  the most  likely  transition path for  a
given  input/output   string.  A  similar   algorithm for
parameterized  logic programs   that determines the  most
likely explanation for  a given goal  can be derived.  It
runs in time linear  in the  size  of the support  graph,
thereby  $O(N^2L)$   in  the  case   of  HMMs,  the  same
complexity as the Viterbi algorithm \cite{Sato00a}.


\subsection{PCFGs}
\label{subsec:PCFGs}

We  now compare the   graphical   EM algorithm with   the
Inside-Outside    algorithm   \cite{Baker79,Lari90}.  The
Inside-Outside algorithm is a well-known EM algorithm for
PCFGs \cite{Wetherell80,Manning99}.\footnote{
%
A PCFG (probabilistic context free grammar) is a backbone
CFG with  probabilities  (parameters)  assigned  to  each
production rule.    For  a   nonterminal $A$  having  $n$
production rules $\{A  \rightarrow \alpha_i \mid  1\leq i
\leq n  \}$,   a probability $p_i$    is assigned  to  $A
\rightarrow  \alpha_i$    $(1\leq  i   \leq    n)$  where
$\sum_{i=1}^n p_i =1$.  The probability of a sentence $s$
is the sum of probabilities of each (leftmost) derivation
of $s$.  The  latter is the  product of probabilities  of
rules used in the derivation.  }
%
It takes a  grammar in {\em Chomsky  normal form}.  Given
$N$ nonterminals, a  production rule in the grammar takes
the  form $i  \rightarrow j, k$  $(1 \leq  i,j,k \leq N)$
(nonterminals are named by   numbers from $1$ to  $N$ and
$1$ is  a starting symbol) or the  form $i \rightarrow w$
where  $1\leq  i\leq N$ and $w$   is a terminal.  In each
iteration,  it  computes  the  inside probability and the
outside probability of every   partial parse tree of  the
given sentence to  update parameters for these production
rules.    Time    complexity is   measured  by   time per
iteration, and  is   described  by  $N$, the   number  of
nonterminals,   and $L$, the    number of terminals in a
sentence.  It  is $O(N^3L^3T)$ for $T$ observed sentences
\cite{Lari90}.

To  compare  the graphical    EM   algorithm  with    the
Inside-Outside  algorithm, we start  from a propositional
program $\db_g = F_g  \cup  R_g$ below representing   the
largest   grammar  containing  all   possible   rules  $i
\rightarrow  j,k$ in $N$  nonterminals  where nonterminal
$1$ is a starting symbol, i.e.\ sentence.

\begin{figure}[htb]
\[
\begin{array}{lcl}
F_{g} & = &
       \begin{array}{l}
       \big\{ \mbox{\tt msw($i$,[$d$,$d'$],[$j$,$k$])}
                \mid 1 \leq i,j,k \leq N,
                       \mbox{$d,d'$ are numbers} \big\} \\
        \hspace{2em}\cup\;
        \big\{ \mbox{\tt msw($i$,$d$,$w$)}
                \mid 1 \leq i \leq N,
                     \mbox{$d$ is a number},
                     \mbox{$w$ is a terminal} \big\}
      \end{array}\\
\\
R_{g} & = &
     \left\{
     \left.     
        \begin{array}{lcl}
           \mbox{\tt q($i$,$d_0$,$d_2$)}
               & \mbox{\tt :-} & \mbox{\tt msw($i$,[$d_0$,$d_2$],[$j$,$k$]),} \\
               &          & \mbox{\tt q($j$,$d_0$,$d_1$),} \\
               &          & \mbox{\tt q($k$,$d_1$,$d_2$).}
        \end{array}
     \right|
        \begin{array}{l}
           1 \leq i, j, k \leq N, \\
           0 \leq d_0 < d_1 < d_2 \leq L
        \end{array}
    \right\} \\
\\
     &   &  \bigcup \;
     \left\{
     \left.
        \begin{array}{lcl}
          \mbox{\tt q($i$,$d$,$d\!+\!1$)}
             & \mbox{\tt :-} & \mbox{\tt msw($i$,$d$,$w_{d+1}$).}
        \end{array}
     \right|  \, 1\leq i\leq N, 0 \leq d \leq L-1
     \right\}
\end{array}
\]
\caption{PCFG program $\db_{g}$}
\label{program:PCFG}
\end{figure}

$\db_g$  is  an   artificial  parsing program  whose sole
purpose is to measure the size of an OLDT tree\footnote{
%
To be more precise, an OLDT structure,  but in this case,
it is a tree   because $\db_{g}$ contains  only constants
(Datalog program) and  there never   occurs the need   of
creating a new root node.}
%
created by the OLDT interpreter when it parses a sentence
$w_1w_2\cdots  w_L$.  So the input sentence $w_1w_2\cdots
w_L$   is embedded  in  the  program  separately  as {\tt
msw($i$,$d$,$w_{d+1}$)}  $(0\leq  d   \leq   L-1)$ in the
second clauses of  $R_g$ (this treatment  does not affect
the complexity argument).  {\tt q($i$,$d_0$,$d_1$)} reads
that the $i$-th nonterminal spans from position $d_0$ to
position  $d_1$,  i.e.\  the   substring $w_{d_0+1}\cdots
w_{d_1}$.  The first  clauses {\tt q($i$,$d_0$,$d_2$)  :-
msw($\cdot$,$\cdot$,$\cdot$),         q($j$,$d_0$,$d_1$),
q($k$,$d_1$,$d_2$)} are  supposed to be textually ordered
according to  the lexicographic order for tuples $\langle
i,j,k,d_0,d_2,d_1 \rangle$.  As a parser, the top-goal is
set to $\leftarrow \mbox{\tt q($1$,$0$,$L$)}$.\footnote{
%
$L$ here is not a Prolog variable but a constant denoting
the sentence length.  }
%
It  asks   the   parser  to  parse  the    whole sentence
$w_1w_2\cdots   w_L$  as the  syntactic category  ``$1$''
(sentence).

We  make  an exhaustive  search  for  this query by  OLDT
search.\footnote{
%
{\tt  q} is a table   predicate.
}
%
As before,  the multi-stage  depth-first search  strategy
and $O(1)$ time access to the solution table are assumed.
Then the time  complexity of OLDT  search is measured  by
the number of nodes in the OLDT tree.
%
\begin{figure}[tb]
\centerline{
	\includegraphics[scale=0.95]{./OLDTtree.eps}
}
\caption{OLDT tree for the query
           $\leftarrow\!\mbox{\tt q($1$,$d$,$L$)}$}
\label{fig:OLDTtree}
\end{figure}
%
Let      $T^{(k)}_d$    be   the     OLDT      tree   for
$\leftarrow\!\mbox{\tt                  q($k$,$d$,$L$)}$.
Figure~\ref{fig:OLDTtree} illustrates $T^{(1)}_d$ for $d$
($0\leq d\leq  L-3$)  where {\tt msw} atoms  are omitted.
As can be seen, the tree has many similar subtrees, so we
put    them   together       (see   {\bf    Note}      in
Figure~\ref{fig:OLDTtree}).   Due to    the   depth-first
strategy,  $T^{(1)}_d$  has   a  recursive structure  and
contains  $T^{(1)}_{d+1}$  as   a subtree.   Nodes  whose
leftmost atom is not underlined are solution nodes, i.e.\
they solve their leftmost atoms for the first time in the
entire   refutation  process.   The underlined atoms  are
already computed in  the subtrees to their left.\footnote{
%
It   can be   inductively   proved that   $T^{(1)}_{d+1}$
contains   every computed {\tt q($i$,$d'$,$d''$)} $(0\leq
d\leq  L-3,$ $d+1 \leq d'< d''  \leq L,$  $1\leq i\leq N,
d''-d' \leq L-d-2)$.}
%
They  only check  the  solution table if  there are their
entries  (= already computed) in  $O(1)$ time.  Since all
clauses are  ground,   such execution only   generates  a
single child node.

We   enumerate   $h^{(1)}_d$,  the  number   of  nodes in
$T^{(1)}_d$ but not in $T^{(k)}_{d+1}$ $(1\leq k\leq N)$.
From    Figure~\ref{fig:OLDTtree},        we          see
$h^{(1)}_d=O(N^3(L-d)^2)$.\footnote{
%
We here  focus on the subtree  $T'_d$.  $j$, $i'$ and $j'$
range from  $1$ to $N$, and $\bigl|\left\{(e,e')\;\left|\;
d+2\le  e'<e\le      L-1\right.\right\}\bigr|=O((L-d)^2)$.
Hence, the number of  nodes in $T'_d$  is $O(N^3(L-d)^2)$.
The  number   of   nodes in $T^{(1)}_d$    but  neither in
$T^{(1)}_{d+1}$  nor  in $T'_d$  is negligible,  therefore
$h^{(1)}_d=O(N^3(L-d)^2)$.  }
%
Let $h^{(k)}_d$ $(2\leq k \leq  N)$ be the number of nodes
in $T^{(k)}_{d+1}$ not contained in $T^{(1)}_{d+1}$. It is
estimated as $O(N^2(L-d-2))$.  Consequently, the number of
nodes  that  are    newly created   in   $T^{(1)}_{d}$  is
$h^{(1)}_d  + \sum^{N}_{k=2} h^{(k)}_d =   O(N^3(L-d)^2)$.
As a  result, total time  for OLDT  search is computed  as
$\sum_{d=0}^{L-3} h_d=O(N^3 L^3)$\footnote{
%
The number of nodes in $T^{(1)}_{L-1}$ and $T^{(1)}_{L-2}$
is negligible.}
%
which  is also  the size of the support graph.\\

We   now consider  a  non-propositional  parsing  program
$\db_{g'}      =     F_{g'}     \cup     R_{g'}$       in
Figure~\ref{program:PDCG}      whose   ground   instances
constitute    the   propositional   program    $\db_{g}$.
$\db_{g'}$ is  a  probabilistic variant  of  DCG  program
\cite{Pereira80} in which {\tt q'}/1, {\tt q'}/6 and {\tt
between}/3 are declared as table predicate.  Semantically
$\db_{g'}$  specifies a probability distribution over the
atoms of the form $\{\mbox{\tt q'($l$)} \mid l \;\mbox{is
a list of terminals} \}$.
%
\begin{figure}[htb]
\[
\begin{array}{lcl}
F_{g'} & = &
       \begin{array}{l}
       \big\{ \mbox{\tt msw($s_i$,$t$,[$s_j$,$s_k$])}
                \mid 1 \leq i,j,k \leq N,
                     \mbox{$t$ is a number} \big\} \\
        \hspace{2em}\cup\;
       \big\{ \mbox{\tt msw($s_i$,$t$,$w$)}
                \mid 1 \leq i \leq N, 
                     \mbox{$t$ is a number},
                     \mbox{$w$ is a terminal} \big\}
      \end{array}\\
\\
R_{g'} & = &
     \left\{
        \begin{array}{l}
           \begin{array}{lcl}
           \mbox{\tt q'(S)}
               & \mbox{\tt :-}&
                     \mbox{\tt length(S,D),}\;
                     \mbox{\tt q'($s_1$,0,D,0,\_,S-[]).}
           \end{array}\\
\\
           \begin{array}{lcl}
           \mbox{\tt q'(I,D0,D2,C0,C2,L0-L2)}
               & \mbox{\tt :-}&
                            \mbox{\tt between(D0,D1,D2),} \\
               &          & \mbox{\tt msw(I,C0,[J,K]),} \\
               &          & \mbox{\tt q'(J,D0,D1,s(C0),C1,L0-L1),} \\
               &          & \mbox{\tt q'(K,D1,D2,s(C1),C2,L1-L2).} \\
           \end{array}\\
\\
           \begin{array}{lcl}
           \mbox{\tt q'(I,D,s(D),C0,s(C0),[W|X]-X)}
              & \mbox{\tt :-} & \mbox{\tt msw(I,C0,W).}
           \end{array}
        \end{array}
     \right.
\end{array}
\]
\caption{Probabilistic DCG like program $\db_{g'}$}
\label{program:PDCG}
\end{figure}

The  top-goal to parse  a  sentence ${\tt  S} = \mbox{\tt
[$w_1,\ldots,w_L$]}$       is      $\leftarrow\!\mbox{\tt
q'([$w_1,\ldots,w_L$])}$.       It invokes      \mbox{\tt
q'($s_1$,0,D,0,\_,[$w_1$,$\ldots$,$w_L$]-[])}       after
measuring the length ${\tt D}$ of an input sentence ${\tt
S}$ by calling ${\tt length}/2$.
\footnote{
%
To make the program as simple as possible, we assume that
an   integer $n$ is   represented by a   ground term $s_n
\defined      \overbrace{\mbox{\tt    s($\cdots$s}}^{(n)}
\mbox{\tt (0)$\cdots$)}$.  We also  assume that when {\tt
D0}   and     {\tt   D2}   are    ground,      the   goal
$\leftarrow\!\mbox{\tt  between(D0, D1,  D2)}$ returns an
integer {\tt  D1} between  them  in time  proportional to
$|{\tt D1}- {\tt D0}|$.  }
%
\footnote{
%
We omit  an obvious  program for  {\tt length($l$,$s_n$)}
which computes the length $s_n$ of a list $l$ in $O(|l|)$
time.  }
%
In                 general,                     \mbox{\tt
q'($i$,$d_0$,$d_2$,$c_0$,$c_2$,$l_0$-$l_2$)}        works
identically to   \mbox{\tt q($i$,$d_0$,$d_2$)}  but three
arguments, $c_0$, $c_2$ and {\tt $l_0$-$l_2$}, are added.
$c_0$  supplies a unique   trial-id for {\tt  msw}s to be
used in the  body,   $c_2$  the latest trial-id    in the
current computation,   and {\tt  $l_0$-$l_2$}   a  D-list
holding a substring  between $d_0$ and $d_2$.  Since  the
added arguments  do not  affect the  shape of  the search
tree    in  Figure~\ref{fig:OLDTtree}    and   the  extra
computation caused by ${\tt  length}/2$ is $O(L)$ and the
one by the insertion of $\mbox{\tt between(D0,D1,D2)}$ is
$O(NL^3)$ respectively,\footnote{
%
$\mbox{\tt between(D0,D1,D2)}$ is   called  $O(N(L-d)^2)$
times in $T^{(1)}_d$.  So  it is called $\sum_{d=0}^{L-3}
O(N(L-d)^2) = O(NL^3)$ times in $T^{(1)}_0$.  }
%
OLDT time remains $O(N^3L^3)$, and hence so is the size of
the support graph.

To apply the graphical EM algorithm correctly, we need to
confirm the five conditions  on its applicability.  It is
rather  apparent however that  the OLDT refutation of any
top-goal   of     the  form     $\leftarrow     \mbox{\tt
q'([$w_1$,$\ldots$,$w_L$])}$      w.r.t.       $\db_{g'}$
terminates,  and  leaves a   support graph satisfying the
finite  support   condition  and  the   acyclic   support
condition.    The   t-exclusiveness condition    and  the
independent condition also  hold  because the  refutation
process   faithfully simulates the leftmost    stochastic
derivation  of $w_1\cdots w_L$  in  which the choice of a
production rule          made            by          {\tt
msw($s_i$,$s_c$,[$s_j$,$s_k$])}   is   exclusive      and
independent  (trial-ids    are   different  on  different
choices).

What remains is the uniqueness  condition.  To confirm it,
let   us    consider   another program    $\db_{g''}$,   a
modification of   $\db_{g'}$  such  that   the first  goal
\mbox{\tt length(S,D)} in the body of the first clause and
the first goal  \mbox{\tt between(D0,D1,D2)} in the second
clause of $R_{g'}$ are moved to the last position in their
bodies    respectively.  $\db_{g''}$   and $\db_{g'}$  are
logically  equivalent, and semantically equivalent as well
from the viewpoint of  distribution semantics.  Then think
of the  sampling execution by  the  OLDT interpreter of  a
top-goal $\leftarrow \mbox{\tt q'(S)}$ w.r.t.\ $\db_{g''}$
where  {\tt  S}   is a variable,    using the  multi-stage
depth-first search strategy.  It is easy to see first that
the execution  never fails, and  second that when the OLDT
refutation  terminates,   a   sentence  {\tt [$w_1,\ldots,
w_L$]} is returned in  {\tt S}, and third that conversely,
the set of {\tt msw} atoms resolved upon in the refutation
uniquely  determines      the output     sentence     {\tt
[$w_1,\ldots,w_L$]}.\footnote{
%
Because the trial-ids used in the refutation record which
rule  is  used at  what  step  in the  derivation of $w_1
\cdots w_L$.}
%
Hence, if the sampling execution  is guaranteed to always
terminate,   every    sampling from   $P_{F_{g''}}$   $(=
P_{F_{g'}})$ uniquely generates a sentence, an observable
atom,  so   the uniqueness  condition   is  satisfied  by
$\db_{g''}$, and hence by $\db_{g'}$.

Then when is the sampling  execution guaranteed to always
terminate?  In other  words, when  does  the grammar only
generate finite sentences?  Giving a general answer seems
difficult,  but it is known  that if the parameter values
in a PCFG are obtained by learning from finite sentences,
the  stochastic  derivation by the  PCFG  terminates with
probability  one   \cite{Chi98}.    In  summary, assuming
appropriate   parameter values,  we   can   say that  the
parameterized  logic program $\db_{g'}$  for  the largest
PCFG  with   $N$   nonterminal  symbols   satisfies   all
applicability  conditions,  and   the OLDT  time   for  a
sentence of length $L$ is $O(N^3L^3)$\footnote{
%
In $\db_{g'}$, we represent integers by ground terms made
out of ${\tt 0}$  and $\mbox{\tt s($\cdot$)}$ to keep the
program short.  If   we  use integers instead of   ground
terms  however, the first  three  arguments of $\mbox{\tt
q'($\cdot$,$\cdot$,$\cdot$,$\cdot$,$\cdot$,$\cdot$)}$ are
enough to check whether the  goal is previously called or
not, and this check can be done in $O(1)$ time.  }
%
and this  is also  the size  of the support  graph. 
From Proposition~\ref{propositon:complexity}, we conclude
%
\begin{proposition}
\label{proposition:PCFG}
Let $\db$ be a parameterized logic program representing a
PCFG with $N$  nonterminals  in the  form  of $\db_g'$ in
Figure~\ref{program:PDCG},      and       ${\cal      G}=
G_1,G_2,\ldots,G_T$ be  the sampled    atoms representing
sentences of length $L$.  We suppose each table operation
in OLDT search is  done in $O(1)$  time. Then OLDT search
for ${\cal G}$ and  one iteration in {\em learn-gEM}  are
respectively done in $O(N^3 L^3 T)$ time.
\end{proposition}

$O(N^3  L^3 T)$  is  also   the  time complexity  of  the
Inside-Outside   algorithm  per iteration  \cite{Lari90},
hence our algorithm is as efficient as the Inside-Outside
algorithm.

\subsection{Pseudo PCSGs}

PCFGs   can      be    improved  by    making     choices
context-sensitive, and  one  of such  attempts  is pseudo
PCSGs   (pseudo probabilistic context sensitive grammars)
in which a rule is  chosen probabilistically depending on
both   the nonterminal  to  be  expanded  and its parent
nonterminal \cite{Charniak94}.

A  pseudo  PCSG   is   easily programmed.   We  add    one
extra-argument, {\tt N}, representing  the parent node, to
the predicate     \mbox{\tt   q'(I,D0,D2,C0,C2,L0-L2)}  in
Figure~\ref{program:PDCG}      and  replace      \mbox{\tt
msw(I,C0,[J,K])}   with  \mbox{\tt   msw([N,I],C0,[J,K])}.
Since  the (leftmost)   derivation  of a  sentence  from a
pseudo    PCSG is still    a  sequential decision  process
described by   the  modified program,    the graphical  EM
algorithm applied to the support graphs generated from the
modified program and observed sentences correctly performs
the ML estimation of parameters in the pseudo PCSG.

A pseudo PCSG is  thought to be a PCFG  with rules of the
form $[n,i] \rightarrow [i,j] [i,k]$ $(1\leq n,i,j,k \leq
N)$ where  $n$ is the parent nonterminal   of $i$, so the
arguments in the   previous subsection are  carried  over
with minor changes.  We therefore have (details omitted)

\begin{proposition}
\label{proposition:PCSG}
Let $\db$ be a  parameterized logic program for  a pseudo
PCSG with $N$  nonterminals as  shown above,  and  ${\cal
G}=G_1,G_2,\ldots,G_T$  the observed  atoms  representing
sampled sentences  of  length  $L$.   Suppose  each table
operation  in  OLDT search can  be done   in $O(1)$ time.
Then OLDT  search  for ${\cal  G}$  and each iteration in
{\em learn-gEM} is completed in $O(N^4 L^3 T)$ time.
\end{proposition}

\subsection{Bayesian Networks}
\label{subsec:Bayes}

A relationship between  cause  $C$ and its effect  $E$ is
often probabilistic such  as the one between diseases and
symptoms, and as such   it is mathematically  captured as
the conditional probability $P(E =e \mid C =c)$ of effect
$e$ given the cause $c$.  What we wish to know however is
the inverse, i.e.\  the probability of a  candidate cause
$c$ given evidence $e$, i.e.\ $P(C=c \mid E= e)$ which is
calculated by Bayes' theorem as  $P(E= e \mid C=c) P(C  =
c)/   \sum_{c'} P(E  =e  \mid  C  =c')P(C=c')$.  Bayesian
networks  are a representational/computational  framework
that fits best this type of probabilistic inference
\cite{Pearl88,Castillo97}.

A  Bayesian network is   a graphical representation of  a
joint  distribution  $P(X_1=x_1,$ $\ldots,$ $X_N=x_N)$ of
finitely many random  variables   $X_1,\ldots,X_N$.   The
graph is a  dag (directed acyclic graph)  such as ones in
Figure~\ref{fig:Bnet},   and each    node   is a   random
variable.\footnote{
%
We only deal with discrete cases.}


\begin{figure}[htb]
\centerline{
	\includegraphics[scale=1.0]{./bn-example.eps}
}
\caption{Bayesian networks}
\label{fig:Bnet}
\end{figure}

In the graph,    a conditional  probability table   (CPT)
representing  $P(X_i  = x_i  \mid {\bpi}_i  =  {\bu}_i )$
$(1\leq i  \leq  N)$ is associated   with each node $X_i$
where ${\bpi}_i$   represents  $X_i$'s parent  nodes  and
$\bu_i$ their values.  When $X_i$ has  no parent, i.e.\ a
topmost node  in the graph, the  table is just a marginal
distribution $P(X_i    =    x_i)$.   The   whole    joint
distribution   is   defined  as the   product   of  these
conditional distributions:
%
\begin{equation}
P(X_1=x_1,\ldots,X_N=x_N)\footnote{
%
Thanks to the  acyclicity of  the graph, without   losing
generality, we may assume that   if $X_i$ is an  ancestor
node of $X_j$, then $i<j$ holds. }
%
   =  \prod_{i=1}^N P(X_i = x_i \mid  {\bpi}_i = \bu_i).
\label{eq:bn:def}
\end{equation}
%
Thus the graph $G_1$ in Figure~\ref{fig:Bnet} defines
%
\[
P_{G_1}(a,b,c,d,e,f)
   = P_{G_1}(a)P_{G_1}(b)P_{G_1}(c\mid  a)P_{G_1}(d\mid   a, b)
                         P_{G_1}(e\mid d)P_{G_1}(f\mid d)
\]
%
where  $a$, $b$, $c$, $d$, $e$    and   $f$  are   values  of
corresponding   random variables $A$, $B$, $C$, $D$, $E$   and
$F$, respectively.\footnote{
%
For notational simplicity, we  shall omit random variables
when no confusion arises.}
%
As mentioned before, one  of the basic tasks of  Bayesian
networks  is  to  compute   marginal  probabilities.  For
example,   the marginal  distribution $P_{G_1}(c,d)$   is
computed      either      by    (\ref{eq:bn:naive})    or
(\ref{eq:bn:factor}) below.
%
\begin{eqnarray}
P_{G_1}(c,d)
     & = & \sum_{a,b,e,f}P_{G_1}(a)P_{G_1}(b)P_{G_1}(c\,|\, a)
                       P_{G_1}(d\,|\, a, b)P_{G_1}(e\,|\, d)P(f\,|\, d)
                \label{eq:bn:naive}  \\
     & = & \left( \sum_{a,b}P_{G_1}(a)P_{G_1}(b) P_{G_1}(c\,|\, a)
                        P_{G_1}(d\, |\, a, b) \right)
           \left( \sum_{e,f}P_{G_1}(e \,|\, d)P_{G_1}(f\, |\, d)\right)
                \label{eq:bn:factor}
\end{eqnarray}

(\ref{eq:bn:factor})  is   clearly    more efficient  than
(\ref{eq:bn:naive}).  Observe that  if the graph were like
$G_2$ in Figure~\ref{fig:Bnet}, there  would be no  way to
factorize computations  like (\ref{eq:bn:factor})  but  to
use   (\ref{eq:bn:naive})  requiring  exponentially   many
operations.   The problem   is   that computing   marginal
probabilities  is  NP-hard  in  general, and factorization
such  as (\ref{eq:bn:factor})  is  assured  only when  the
graph is {\em singly connected\/} like $G_1$, i.e.\ has no
loop when viewed as   an undirected graph.  In  such case,
the computation is possible in $O(|V|)$  time where $V$ is
the   set  of vertices     in  the graph   \cite{Pearl88}.
Otherwise, the graph  is  called {\em multiply-connected},
and might   need exponential    time to compute   marginal
probabilities.  In the sequel, we show the following.

\begin{itemize}

\item
For  any   discrete  Bayesian   network  $G$   defining a
distribution     $P_G(x_1,\ldots,x_N)$,   there   is    a
parameterized logic program  ${\db}_{G}$  for a predicate
${\tt  bn}(\cdot)$  such  that   $P_{{\db}_{G}}(\mbox{\tt
bn($x_1$,$\ldots$,$x_N$)})$ $ = P_G(x_1,\ldots,x_N)$.

\item
For arbitrary factorizations and their order to compute a
marginal distribution, there exists a tabled program that
accomplishes the same computation in the specified way.

\item
When the  graph is singly connected  and  evidence $e$ is
given, there  exists  a  tabled program  $\db_{G^{\tau}}$
such that    OLDT time for  $\leftarrow \!   \mbox{\tt
bn($e$)}$ is $O(|V|)$, and  hence the time complexity per
iteration of the  graphical  EM algorithm is $O(|V|)$  as
well.

\end{itemize}

Let $G$    be   a Bayesian   network   defining a   joint
distribution  $P_G(x_1,\ldots,x_N)$  and $\{  P_G(X_i=x_i
\mid {\bpi}_i = \bu_i) \mid $ $1\leq i  \leq N,$ $x_i \in
val(X_i),$  $\bu_i  \in  val(\bpi_i) \}$  the conditional
probabilities associated with $G$ where $val(X_i)$ is the
set   of $X_i$'s   possible   values and  $val({\bpi}_i)$
denotes the set of  possible  values of the  parent nodes
${\bpi}_i$  as   a  random   vector,  respectively.    We
construct a parameterized  logic program that defines the
same distribution  $ P_G(x_1,\ldots,x_N)$.  Our   program
$\db_{G}   = F_{G}     \cup   R_{G}$   is    shown     in
Figure~\ref{program:bn}.
%
\begin{figure}[hbt]
\[
\begin{array}{lcl}
F_{G} & = &
       \big\{\; \mbox{\tt msw(par($i$,${\bu}_{i}$),once,$x_{i}$) }
                \mid 1 \leq i \leq N, {\bu}_{i} \in val({\bpi}_i),
                       x_i \in val(X_i) \; \big\}\\
\\
R_{G} & = &
       \big\{\; \mbox{\tt bn(${\tt X}_1$,$\ldots$,${\tt X}_N$)}
              \mbox{\tt :-} \bigwedge_{i=1}^N
                \mbox{\tt msw(par($i$,${\bpi}_i$),once,${\tt X}_i$).}
        \; \big\} 
\end{array}
\]
\caption{Bayesian network program $\db_{G}$}
\label{program:bn}
\end{figure}

\noindent
$F_{G}$ is comprised of {\tt msw}  atoms of the form {\tt
msw(par($i$,${\bu}_{i}$),once,$x_{i}$)} whose probability
is  exactly the conditional probability $P_G(X_i=x_i \mid
{\bpi}_i =  \bu_i)$.     When  $X_i$  has   no   parents,
${\bu}_{i}$  is the empty list  {\tt  []}.  $R_{G}$ is  a
singleton, containing  only one  clause whose body   is a
conjunction  of {\tt msw}  atoms which corresponds to the
product of conditional probabilities.   Note that we {\em
intentionally\/}      identify       random     variables
$X_1,\ldots,X_N$    with   logical   variables      ${\tt
X}_1,\ldots,{\tt X}_N$ for convenience.

\begin{proposition}
\label{proposition:Bayes}
${\db_{G}}$ denotes the same distributions as $G$.
\end{proposition}
%
(Proof) Let  $\langle  x_{1},\ldots,x_{N}  \rangle$ be   a
realization      of     the    random    vector   $\langle
X_{1},\ldots,X_{N} \rangle$.    It  holds by  construction
that
%
\begin{eqnarray*}
P_{{\db}_{G}}(\mbox{\tt bn($x_{1}$,$\ldots$,$x_{N}$)})
  & = &
     \prod_{h=1}^N
         P_{\tt msw}(\mbox{\tt msw(par($i$,${\bu}_{i}$),once,$x_{i}$)}) \\
 & = & \prod_{h=1}^N P_G(X_i  = x_{i}   \mid {\bpi}_i = \bu_i) \\
 & = & P_G(x_{1},\ldots,x_{N}). \hspace{3em} Q.E.D.
\end{eqnarray*}
%
In the case of $G_1$ in Figure~\ref{fig:Bnet}, the program
becomes\footnote{
%
${\tt 'A'}, {\tt 'B'},\ldots$ are Prolog constants used in
place of integers.}
%
\begin{center}
{\tt
\begin{tabular}{llll}
 bn(A,B,C,D,E,F)
   & :- &  msw(par('A',[]),once,A),  & msw(par('B',[]),once,B), \\
   &    &  msw(par('C',[A]),once,C), & msw(par('D',[A,B]),once,D), \\
   &    &  msw(par('E',[D]),once,E), & msw(par('F',[D]),once,F).
\end{tabular}
}
\end{center}
%
and  the  left-to-right sampling execution  gives a sample
realization of the random  vector $\langle\, {\tt A, B, C,
D, E, F}\, \rangle$.  A marginal distribution  is computed
from \mbox{\tt   bn($x_{1}$,$\ldots$,$x_{N}$)} by adding a
new clause   to  ${\db}_G$.    For  example,  to   compute
$P_{G_1}(c,d)$, we add {\verb!bn(C,D):- bn(A,B,C,D,E,F)!}
to ${\db_{G_1}}$  (let  the result be  ${\db}_{G_1}'$) and
then  compute  $P_{{\db}_{G_1}'}(\mbox{\tt  bn($c$,$d$)})$
which is equal to $P_{G_1}(c,d)$ because
%
\begin{eqnarray*}
P_{{\db}_{G_1}'}(\mbox{\tt bn($c$,$d$)})
	& = & 
      P_{{\db}_{G_1}}(\exists\;a,b,e,f\;
                         \mbox{\tt bn($a$,$b$,$c$,$d$,$e$,$f$)}) \\
  &  = &
     \sum_{a,b,e,f}P_{{\db}_{G_1}}(\mbox{\tt bn($a$,$b$,$c$,$d$,$e$,$f$)}) \\
  & = & P_{G_1}(c,d).
\end{eqnarray*}
%
Regrettably    this    computation    corresponds     to
(\ref{eq:bn:naive}),    not    to  the     factorization
(\ref{eq:bn:factor}).  Efficient probability computation
using factorization is   made possible by   carrying out
summations in a proper order.\\

We next sketch  by an example  how to carry  out specified
summations    in a  specified   order  by  introducing new
clauses.  Suppose we have a joint distribution $P(x,y,z,w)
=$   $\phi_1(x,y) \phi_2(y,z,w)  \phi_3(x,z,w)$ such  that
$\phi_1(x,y)$, $\phi_2(y,z,w)$ and   $\phi_3(x,z,w)$   are
respectively   computed  by   atoms     {\tt   \mbox{${\tt
p}_1$}(X,Y)},   {\tt \mbox{${\tt   p}_2$}(Y,Z,W)} and {\tt
\mbox{${\tt p}_3$}(X,Z,W)}.  Suppose also that we  hope to
compute the sum
%
\[
P(x) = \sum_{y}
   \phi_1(x,y)\left(\sum_{z,w} \phi_2(y,z,w)\phi_3(x,z,w)\right)
\]
%
in  which    we  first eliminate   $z,w$    and  then $y$.
Corresponding to each   elimination, we introduce two  new
predicates,   {\tt  q(X,Y)}   to compute  $\phi_4(x,y)  =$
$\sum_{z,w} \phi_2(y,z,w)\phi_3(x,z,w)$  and {\tt p(X)} to
compute  $P(x)    =$  $\sum_{y}\phi_1(x,y)\phi_4(x,y)$  as
follows.
%
\begin{center}
{\tt 
\begin{tabular}{rcl}
p(X)
  & :- & \mbox{${\tt p}_1$}(X,Y), q(X,Y). \\
q(X,Y)
  & :- & \mbox{${\tt p}_2$}(Y,Z,W), \mbox{${\tt p}_3$}(X,Z,W).
\end{tabular}
}
\end{center}
%
Note that the clause body  of  ${\tt q}/2$ contains  ${\tt
Z}$ and  ${\tt  W}$  as (existentially  quantified)  local
variables and the    clause  head {\tt  q(X,Y)}   contains
variables shared  with {\em other  atoms}.  In view of the
correspondence between $\sum$ and $\exists$, it is easy to
confirm  that   this   program    realizes  the   required
computation.  It is also  easy to see by generalizing this
example, though we do not  prove here, that there exists a
parameterized logic program  that  carries out the   given
summations in  the given order   for an arbitrary Bayesian
network, in   particular we  are     able to simulate   VE
(variable  elimination, \citeR{Zhang96,Dambrosi99}) in our
approach.\\

Efficient computation  of  marginal  distributions  is not
always  possible but  there    is a well-known  class   of
Bayesian networks, singly connected Bayesian networks, for
which  there  exists  an  efficient algorithm  to  compute
marginal  distributions   by         message       passing
\cite{Pearl88,Castillo97}.  We   here  show that  when the
graph is singly connected,  we can construct  an efficient
tabled Bayesian network program $\db_{G^{\tau}}$ assigning
a  table predicate to  each node.  To avoid complications,
we   explain  the  construction   procedure informally and
concentrate on the case  where we have only one interested
variable.  Let $G$ be a  singly connected graph.  First we
pick  up a node $U$ whose  probability $P_G(u)$ is what we
seek.  We construct a tree  $G^{\tau}$ with the root  node
$U$ from $G$, by letting   other nodes dangling from  $U$.
Figure~\ref{fig:bn:trans1} shows how $G_1$ is  transformed
to a tree when we select node B as the root node.

\begin{figure}[hbt]
\centerline{
	\includegraphics[scale=0.95]{./bn-trans1.eps}
}
\caption{Transforming $G_1$ to a tree}
\label{fig:bn:trans1}
\end{figure}

Then  we examine each node  in $G^{\tau}$ one  by one.  We
add for each node $X$  in the graph a corresponding clause
to  $\db_{G^{\tau}}$ whose purpose  is to  visit all nodes
connected to  $X$ except the one that  calls $X$.  Suppose
we    started     from  the     root      node  $U_1$   in
Figure~\ref{fig:bn:trans2} where  evidence $u$  is  given,
and  have  generated  clause (\ref{b-cl:top}).    Now   we
proceed to  an inner node  $X$ ($U_1$ calls $X$).   In the
original graph $G$,  $X$ has parent nodes  $\{U_1,U_2, U_3
\}$ and child nodes $\{  V_1,V_2 \}$.  $U_3$ is a  topmost
node in $G$.

\begin{figure}[htb]
\centerline{
	\includegraphics[scale=0.9]{./bn-trans2.eps}
}
\caption{General situation}
\label{fig:bn:trans2}
\end{figure}

For node  $X$ in Figure~\ref{fig:bn:trans2}, we add clause
(\ref{b-cl:in}).   When it is  called from the parent node
$U_1$ with ${\tt U}_1$ being   ground, we first   generate
possible    values   of     $U_2$     by  calling     {\tt
val\_${U_2}$(${\tt   U}_2$)},   and    then   call    {\tt
call\_$X$\_$U_2$(${\tt     U}_2$)}   to  visit   all nodes
connected to    $X$  through $U_2$.   $U_3$   is  similary
treated.  After visiting  all nodes  in $G$ connecting  to
$X$ through   the  parent nodes  $U_2$  and  $U_3$  (nodes
connected to $U_1$  have already been visited),  the value
of random variable $X$ is determined  by sampling the {\tt
msw} atom  jointly indexed by '{\tt  X}' and the values of
$U_1$, $U_2$  and $U_3$.   Then we visit   $X$'s children,
$V_1$ and $V_2$.  For a topmost node $U_3$ in the original
graph, we add clause (\ref{b-cl:leaf}).

\begin{eqnarray}
   \mbox{\tt  tbn(${\tt U}_1$)}
   & \mbox{\tt :-} &
           \mbox{\tt  msw(par('${\tt U}_1$',[]),once,${\tt U}_1$),\,\,}
           \mbox{\tt call\_$U_1$\_$X$(${\tt U}_1$).}
       \label{b-cl:top} \\
\nonumber \\
\mbox{\tt call\_$U_1$\_$X$(${\tt U}_1$)}
   & \mbox{\tt :-} &
       \mbox{\tt val\_${U_2}$(${\tt U}_2$), call\_$X$\_$U_2$(${\tt U}_2$),}
          \nonumber \\
   &           &
       \mbox{\tt val\_${U_3}$(${\tt U}_3$),  call\_$X$\_$U_3$(${\tt U}_3$),}
           \nonumber \\
   &          &
    \mbox{\tt  msw(par('X',[${\tt U}_1$,${\tt U}_2$,${\tt U}_3$]),once,X),}
            \nonumber  \\
   &           &
      \mbox{\tt call\_$X$\_{$V_1$}({\tt X}), call\_$X$\_{$V_2$}({\tt X}).} 
       \label{b-cl:in} \\
\nonumber \\
\mbox{\tt call\_$X$\_$U_3$(${\tt U}_3$)}
   & \mbox{\tt :-} &
      \mbox{\tt  msw(par('${\tt U}_3$',[]),once,${\tt U}_3$).}
       \label{b-cl:leaf}
\end{eqnarray}

Let   $\db_{G^{\tau}}$ be   the final program   containing
clauses  like  (\ref{b-cl:top}),   (\ref{b-cl:in})     and
(\ref{b-cl:leaf}).  Apparently $\db_{G^{\tau}}$  can    be
constructed  in time linear in the  number of nodes in the
network.      Also    note  that     successive  unfolding
\cite{Tamaki84} of      atoms   of   the     form     {\tt
call\_...($\cdot$)} in the clause  bodies that starts from
(\ref{b-cl:top}) yields  a  program $\db'_{G}$ similar  to
the one  in  Figure~\ref{program:bn} which   contains {\tt
msw}  atoms   but   no  {\tt   call\_...($\cdot$)}'s.   As
$\db_{G^{\tau}}$  and     $\db'_{G}$   define  the    same
distribution,\footnote{
%
Since  distribution semantics is based  on the least model
semantics, and because        unfold/fold   transformation
\cite{Tamaki84} preserves the  least Herbrand model of the
transformed program, unfold/fold transformation applied to
parameterized  logic programs  preserves the denotation of
the transformed program.  }
%
it can be proved from  Proposition~\ref{proposition:Bayes}
that   $P_G(u)    =    P_{\db'_{G}}({\tt     bn}(u))     =
P_{\db_{G^{\tau}}}({\tt tbn}(u))$ holds (details omitted).
By the   way, in Figure~\ref{fig:bn:trans2} we  assume the
construction  starts   from the  topmost  node ${\tt U}_1$
where  the evidence  $u$   is   given,  but this  is   not
necessary.  Suppose we change to start from the inner node
{\tt X}.  In that case, we replace clause (\ref{b-cl:top})
with \mbox{\tt call\_$X$\_$U_1$(${\tt U}_1$)} \mbox{\tt :-}
\mbox{\tt  msw(par('${\tt   U}_1$',[]),once,${\tt  U}_1$)}
just like (\ref{b-cl:leaf}).  At the  same time we replace
the head of  clause (\ref{b-cl:in}) with  \mbox{\tt tbn()}
and add  a  goal  \mbox{\tt call\_$X$\_$U_1$($u$)} to  the
body     and     so on.     For    the   changed   program
${\db''}_{G^{\tau}}$,  it  is rather  straightforward   to
prove that $P_{{\db''}_{G^{\tau}}}(\mbox{\tt   tbn()})   =
P_G(u)$  holds.  It is  true that the  construction of the
tabled program  $\db_{G^{\tau}}$ shown here is  very crude
and there  is a  lot   of room for  optimization, but   it
suffices to show that a  parameterized logic program for a
singly connected  Bayesian  network runs  in $O(|V|)$ time
where $V$ is the set of nodes.\\

To  estimate time   complexity  of  OLDT  search   w.r.t.\
$\db_{G^{\tau}}$, we declare {\tt tbn} and every predicate
of the form   {\tt call\_...($\cdot$)} as  table predicate
and verify the five conditions on the applicability of the
graphical EM algorithm (details omitted).  We now estimate
the   time   complexity  of  OLDT  search   for  the  goal
$\leftarrow\!     \mbox{\tt           tbn($u$)}$   w.r.t.\
$\db_{G^{\tau}}$.\footnote{
%
$\db_{G^{\tau}}$ is further   transformed  for  the  OLDT
interpreter  to collect {\tt  msw} atoms like the case of
the HMM program.  }
%
We notice that calls occur according to the pre-order scan
(parents -- the node --  children) of the tree $G^{\tau}$,
and  calls to {\tt call\_$Y$\_$X$($\cdot$)} occur $val(Y)$
times.  Each call to {\tt call\_$Y$\_$X$($\cdot$)} invokes
calls to the   rest of  nodes,  $X$'s  parents  and  $X$'s
children in the  graph $G^{\tau}$ except the  caller node,
with diffrent set of variable instantiations, but from the
second call on, every call only refers to solutions stored
in the solution table in $O(1)$ time.  Thus, the number of
added computation  steps in OLDT  search by $X$ is bounded
from  above, by constant  $O(val(U1)val(U2)val(U3)val(X))$
in the case   of Figure~\ref{fig:bn:trans2}.  As a  result
OLDT time is  proportional to the number  of nodes  in the
original  graph  $G$ (and  so  is the size of  the support
graph) provided  that the number  of edges connecting to a
node, and that of values of  a random variable are bounded
from above. So we have

\begin{proposition}
\label{proposition:Bnet}

Let $G$ be a  singly  connected Bayesian network  defining
distribution $P_G$,     $V$   the set     of   nodes,  and
$\db_{G^{\tau}}$ the   tabled   program derived  as above.
Suppose the number of edges connecting to a node, and that
of values of  a random variable  are bounded from above by
some constant.   Also suppose table access  can be done in
$O(1)$ time.  Then,  OLDT time for  computing $P_G(u)$ for
an observed value $u$ of a random variable $U$ by means of
$\db_{G^{\tau}}$ is $O(|V|)$  and so is time per iteration
required by the graphical EM algorithm.   If there are $T$
observations, time complexity is $O(|V|T)$.
\end{proposition}

$O(|V|)$  is the time    complexity required to compute  a
marignal  distribution for  a  singly  connected  Bayesian
network by a standard algorithm \cite{Pearl88,Castillo97},
and also  is   that of  the  EM algorithm   using  it.  We
therefore conclude that  the graphical EM algorithm  is as
efficient   as  a  specialzed   EM  algorithm  for  singly
connected Bayesian networks.\footnote{
%
When a marginal distribution  of $P_G$ for more  than one
variable  is required, we  can construct a similar tabled
program that  computes  marginal probabilities   still in
$O(|V|)$ time by adding extra-arguments that convey other
evidence  or   by   embedding   other evidnece   in   the
program.}
%
We  must also quickly add that  the graphical EM algorithm
is applicable to arbitrary Bayesian networks,\footnote{
%
We    check  the   five  conditions    with   $\db_{G}$ in
Figure~\ref{program:bn}.    The uniqueness  condition   is
obvious  as sampling always  uniquely generates  a sampled
value  for   each random  variable.   The   finite support
condition is  satisfied  because there are only   a finite
number of random variables and  their values.  The acyclic
support condition is  immediate because of  the acyclicity
of Bayesian  networks.  The  t-exclusiveness condition and
the independent condition are easy to verify.  }
%
and  what Proposition~\ref{proposition:Bnet}  says is that
an  explosion of the    support graph can  be  avoided  by
appropriate  programming in the  case  of singly connected
Bayesian networks.\\

To summarize, the graphical EM algorithm, a single generic
EM  algorithm, is proved to  have the same time complexity
as   specialized   EM algorithms,   i.e.\   the Baum-Welch
algorithm   for HMMs,  the    Inside-Outside algorithm for
PCFGs, and the one  for singly connected Bayesian networks
that have been developed   independently in each  research
field.

Table~\ref{table:complexity}   summarizes    the      time
complexity of  EM   learning using   OLDT search and   the
graphical EM algorithm in the case of one observation.  In
the  first  column, ``sc-BNs'' represents singly connected
Bayesian  networks.  The second  column shows a program to
use.        $\db_h$       is   an        HMM proram     in
Subsection~\ref{subsec:OLDT}, $\db_{g'}$ a PCFG program in
Subsection~\ref{subsec:PCFGs}    and     $\db_{G^\tau}$  a
transformed      Bayesian    network    program         in
Subsection~\ref{subsec:Bayes}, respectively.  OLDT time in
the third  column is time  for OLDT search to complete the
search of all t-explanations.  gEM in the fourth column is
time in one iteration taken by  the graphical EM algorithm
to update parameters.  We   use  $N$,  $M$, $L$ and    $V$
respectively for   the number of   states in an   HMM, the
number of nonterminals in a  PCFG, the length of an input
string and the number of nodes in a Bayesian network.  The
last  column is a  standard (specialized) EM algorithm for
each model.

\begin{table}[hbt]
\begin{center}
%\begin{footnotesize}
{\def\arraystretch{1.2}
\begin{tabular}{l|c|c|c|c}
\hline
Model	& Program			& OLDT time		& gEM	& Specialized EM\\
\hline
HMMs	& $\db_h$			& $O(N^2L)$		& $O(N^2L)$		& Baum-Welch \\
PCFGs	& $\db_{g'}$		& $O(M^3L^3)$	& $O(M^3L^3)$	& Inside-Outside \\
%PCSGs	& 					& $O(M^4L^3)$	& $O(M^4L^3)$	& \\
sc-BNs	& $\db_{G^\tau}$	& $O(|V|)$		& $O(|V|)$
							& \cite{Castillo97} \\
user model
		& 	& $O(|\mbox{OLDT tree}|)$
			& $O(|\mbox{support graph}|)$ \\
\hline
\end{tabular}
}
%\end{footnotesize}
\end{center}
\caption{Time complexity of EM learning by OLDT search
            and the graphical EM algorithm}
\label{table:complexity}
\end{table}
%

\subsection{Modeling Language PRISM}
\label{sec:sub:PRISM}

We have   been developing a  symbolic-statistical modeling
laguage     {\sc    PRISM}       since    1995  (URL     =
\verb+http://mi.cs.titech.ac.jp/prism/+)       as       an
implementation        of          distribution   semantics
\cite{Sato95,Sato97,Sato98}.  The language is intented for
modeling  complex  symbolic-statistical phenomena  such as
discourse interpretation    in natural language processing
and gene inheritance interacting  with social rules.  As a
programming language, it looks like an extension of Prolog
with new   built-in  predicates including  the   {\tt msw}
predicate and  other special predicates  for  manipulating
{\tt msw} atoms and their parameters.

A {\sc PRISM} program is comprised of three parts, one for
directives,  one for modeling and  one for utilities.  The
directive part contains declarations such as {\tt values},
telling the system what {\tt  msw} atoms  will be used  in
the  execution.  The modeling   part is a  set of non-unit
definite clauses that define the distribution (denotation)
of the program by  using {\tt msw}  atoms.  The last part,
the utility  part,  is an  arbitary  Prolog  program which
refers to predicates defined in the modeling part.  We can
use in the utility part {\tt  learn} built-in predicate to
carry out EM learning from observed atoms.

{\sc  PRISM} provides  three   modes  of  execution.   The
sampling execution  correponds to a random  sampling drawn
from the distribution defined  by the modeling part.   The
second one computes the probability  of a given atom.  The
third one returns the support set for a given goal.  These
execution modes are available through built-in predicates.

We  must report however that  while  the implementation of
the  graphical EM algorithm   with a simpified OLDT search
mechanism has been under way, it is not completed yet.  So
currently,   only    Prolog         search     and    {\em
learn-naive\/}($\db,{\cal            G}$)               in
Section~\ref{sec:graphical-EM}    are  available for    EM
learning though we realized, partially, structrure sharing
of   explanations     in   the  implemention     of   {\em
learn-naive\/}($\db,{\cal    G}$).   Putting computational
efficiecy aside however, there is no problem in expressing
and learning  HMMs, PCFGs, pseudo PCSGs, Bayesian networks
and other probailistic models by the current version.  The
learning experiments in the next section  used a parser as
a  substitute   for   the   OLDT    interpreter, and   the
independently implemented graphical EM algorithm.


\section{Learning Experiments}
\label{sec:experiments}

After complexity  analysis of the  graphical EM  algorithm
for  popular symbolic-probabilistic models in the previous
section, we look at an actual behavior of the graphical EM
algorithm with real   data in this section.  We  conducted
learning experiments with   PCFGs using two corpora  which
have  contrasting characters, and compared the performance
of  the graphical    EM algorithm  against  that   of  the
Inside-Outside algorithm in terms of time per iteration (=
time for updating parameters).   The results indicate that
the   graphical    EM   algorithm   can  outperform    the
Inside-Outside algorithm by orders  of magnigude.  Detalis
are  reported by   \citeA{Sato01b}.  Before proceeding,  we
review the Inside-Outside algorithm for completeness.

\subsection{The Inside-Outside Algorithm}
\label{subsec:I-O-algorithm}

The     Inside-Outside    algorithm  was    proposed    by
\citeA{Baker79}   as a  generalization   of the Baum-Welch
algorithm to PCFGs.  The algorithm is designed to estimate
parameters    for a  CFG  grammar  in  Chomsky normal form
containing rules expressed  by numbers like $i \rightarrow
j,k$ ($1\leq i,j,k \leq N$ for $N$ nonterminals, where $1$
is a   starting  symbol).   Suppose  an    input  sentence
$w_1,\ldots,w_L$ is given.   In  each iteration, it  first
computes   in  a bottom  up  manner  inside  probabilities
$e(s,t,i)$       $=    P(i       \stackrel{*}{\Rightarrow}
w_s,\ldots,w_t)$ and then computes outside   probabilities
$f(s,t,i)$          $=P(S        \stackrel{*}{\Rightarrow}
w_1,\ldots,w_{s-1}\;   i\;  w_{t+1},  \ldots,  w_L)$ in  a
top-down manner for every $s$,  $t$ and $i$ $(1\leq s \leq
t  \leq  L,  1\leq  i  \leq N)$.    After  computing  both
probabilities, parameters are  updated by using  them, and
this process iterates  until  some predetermined criterion
such  as  a  convergence  of the  likelihood of  the input
sentence is  achieved.  Although  Baker  did not give  any
analysis of  the Inside-Outside algorithm,  \citeA{Lari90}
showed that it takes $O(N^3L^3)$ time in one iteration and
\citeA{Lafferty93} proved that it is the EM algorithm.


While it  is true  that  the Inside-Outside algorithm  has
been  recognized as a standard  EM  algortihm for training
PCFGs, it is notoriously slow.  Although there is not much
literature  explicitly   stating time  required   by   the
Inside-Outside     algorithm      \cite{Carroll98,Beil99},
\citeA{Beil99} reported for example that when they trained
a  PCFG with 5,508  rules for  a  corpus of 450,526 German
subordinate   clauses  whose  average  ambiguity is  9,202
trees/clause      using   four     machines        (167MHz
Sun~UltraSPARC$\times          2$          and      296MHz
Sun~UltraSPARC-II$\times   2$),  it   took  2.5 hours   to
complete {\em one iteration}.   We  discuss later why  the
Inside-Outside algorithm is slow.

\subsection{Learning Experiments Using Two Corpora}
\label{subsec:ATR-EDR}

We report here  parameter learning of existing PCFGs using
two corpora of moderate size and  compare the graphical EM
algorithm against the Inside-Outside algorithm in terms of
time per  iteration.  As mentioned before, support graphs,
input to the garphical  EM algorithm, were generated  by a
parser, i.e.\ MSLR parser.\footnote{
%
MSLR parser is  a Tomita (Generalized LR) parser developed
by    Tanaka-Tokunaga Laboratory  in   Tokyo Institute  of
Technology \cite{Tanaka97}.  }
%
All  measurements were made  on a 296MHz Sun~UltraSPARC-II
with 2GB memory under Solaris~2.6 and the threshold for an
increase of the $\log$  likelihood of input sentences  was
set to  $10^{-6}$ as  a  stopping  criterion  for the   EM
algorithms.

In  the experiments, we  used  ATR  corpus and EDR  corpus
(each converted to a  POS (part of speech)-tagged corpus).
They are similar in size (about 10,000) but contrasting in
their characters, sentence  length and ambiguity  of their
grammars.  The first  experiment employed ATR corpus which
is a  Japanese-English corpus (we  used only  the Japanese
part) developed   by  ATR \cite{Uratani94}.   It  contains
10,995   short  conversational   sentences, whose  minimum
length, average length and maximum length are respectively
2,  9.97 and 49.  As   a skeleton of   PCFG, we employed a
context free grammar $G_{\mbox{atr}}$ comprising 860 rules
(172 nonterminals and  441 terminals) manually  developed
for  ATR   corpus    \cite{Tanaka97}  which    yields  958
parses/sentence.

Because the Inside-Outside algorithm only accepts a CFG in
Chomsky normal  form, we converted  $G_{\mbox{atr}}$  into
Chomsky  normal           form         $G_{\mbox{atr}}^*$.
$G_{\mbox{atr}}^*$ contains 2,105 rules (196  nonterminals
and   441  terminals).  We  then  divided  the corpus into
subgroups of similar length like $(L=1,2), (L=3,4),\ldots,
(L=25,26)$, each containing randomly chosen 100 sentences.
After these preparations,  we  compare at each length  the
graphical  EM   algorithm applied to  $G_{\mbox{atr}}$ and
$G_{\mbox{atr}}^*$  against the  Inside-Outside  algorithm
applied    to $G_{\mbox{atr}}^*$ in    terms  of  time per
iteration by running them until convergence.


\begin{figure*}[htb]
\centerline{
  \includegraphics[scale=0.45]{./atr-io.eps}
  \hspace{1em}
  \includegraphics[scale=0.45]{./atr-io-gem.eps}
  \hspace{1em}
  \includegraphics[scale=0.45]{./atr-gem-magnified.eps}
}
\caption{Time per iteration : I-O vs. gEM (ATR)}
\label{fig:ATR-I-O-gEM}
\end{figure*}

Curves in  Figure~\ref{fig:ATR-I-O-gEM} show the  learning
results  where an  x-axis is  the  length $L$ of an  input
sentence  and a  y-axis is  average  time taken by  the EM
algorithm   in one iteration    to  update all  parameters
contained in the  support graphs generated from the chosen
100  sentences (other  parameters  in the  grammar do  not
change).  In the left  graph, the Inside-Outside algorithm
plots a  cubic curve labeled ``I-O''.   We omitted a curve
drawn by the graphical EM algorithm as it drew the x-axis.
The  middle graph  magnifies the   left graph.  The  curve
labeled ``gEM (original)''  is plotted by the graphical EM
algorithm applied to the original grammar $G_{\mbox{atr}}$
whereas  the  one   labeled   ``gEM  (Chomsky NF)''   used
$G_{\mbox{atr}}^*$.   At  length 10,  the average sentence
length, it is measured that whichever grammar is employed,
the graphical EM   algorithm runs several   hundreds times
faster  (845 times faster in  the case of $G_{\mbox{atr}}$
and 720  times  faster in the  case of $G_{\mbox{atr}}^*$)
than  the  Inside-Outside  algorithm per  iteration.   The
right graph  shows (almost) linear  dependency of updating
time  by the graphical  EM  algorithm within the measuared
sentence length.\\

Although some  difference is anticipated in their learning
speed, the speed  gap between the Inside-Outside algorithm
and the graphical EM algorithm is unexpectedly large.  The
most conceivable reason is  that ATR corpus  only contains
short sentences and $G_{\mbox{atr}}$ is not much ambiguous
so  that parse trees  are   sparse and  generated  support
graphs are  small, which affects favorably the perforamnce
of the graphical EM algorithm.

We therefore  conducted the  same experiment  with another
corpus which contains  much longer sentences  using a more
ambiguous grammar  that generates dense  parse  trees.  We
used EDR Japanese  corpus \cite{edr95} containing  220,000
Japanese news article sentences.  It  is however under the
process of re-annotation,  and only  part of it  (randomly
sampled 9,900 sentences)  has recently been made available
as a labeled corpus.   Compared with ATR corpus, sentences
are much longer (the average length  of 9,900 sentences is
20, the minimum length 5, the maximum length 63) and a CFG
grammar $G_{\mbox{edr}}$   (2,687 rules,   converted    to
Chomsky normal  form grammar $G_{\mbox{edr}}^*$ containing
12,798 rules) developed for it is  very ambiguous (to keep
a coverage rate),  having $3.0\times 10^8$ parses/sentence
at length 20 and $6.7  \times 10^{19}$ parses/sentence  at
length 38.

\begin{figure*}[htb]
\centerline{
   \includegraphics[scale=0.45]{./edr-io.eps}
  \hspace{1em}
  \includegraphics[scale=0.45]{./edr-io-gem.eps}
  \hspace{1em}
  \includegraphics[scale=0.45]{./edr-gem-magnified.eps}
}
\caption{Time per iteration : I-O vs. gEM (EDR)}
\label{fig:EDR-I-O-gEM}
\end{figure*}

Figure~\ref{fig:EDR-I-O-gEM}  shows  the  obtained  curves
from the experiments  with  EDR corpus (the  graphical  EM
algorithm      applied to     $G_{\mbox{edr}}$   vs.   the
Inside-Outside   algorithm  applied to $G_{\mbox{edr}}^*$)
under the   same condition as  ATR  corpus, i.e.\ plotting
average time per iteration to process 100 sentences of the
designated length, except that   the plotted time for  the
Inside-Outside  algorithm is the  average of 20 iterations
whereas that for the graphical EM algorithm is the average
of 100  iterations.   As is  clear from the  middle graph,
this time again, the graphical EM algorithm runs orders of
magnitude faster than  the  Inside-Outside algorithm.   At
average sentence length 20, the  former takes 0.255 second
whereas the latter takes 339 seconds, giving a speed ratio
of  1,300 to 1.  At  sentence length 38,  the former takes
2.541 seconds but the latter takes 4,774 seconds, giving a
speed ratio  of  1,878 to  1.  Thus  the speed  ratio even
widens compared  to ATR corpus.  This  can be explained by
the mixed effects  of   $O(L^3)$, time complexity of   the
Inside-Outside algorithm,  and a  moderate increase in the
total size of support graphs w.r.t.\ $L$.  Notice that the
right  graph  shows how  the total  size of support graphs
grows  with sentence length $L$  as  time per iteration by
the graphical EM algorithm is linear  in the total size of
support graphs.\\

Since  we  implemented  the  Inside-Outside      algorithm
faithfully to  \citeA{Baker79,Lari90},  there is much room
for   improvement.   Actually    Kita   gave a     refined
Inside-Outside algorithm  \cite{Kita99}.  There is also an
implementation  by  Mark  Johnson   of  the Inside-Outside
algorithm                 down-loadable               from
{\verb!http://www.cog.brown.edu/%7Emj/!}.  The use of such
implementations  may  lead  to different  conclusions.  We
therefore conducted  learning  experiments  with the  {\em
entire\/} ATR corpus  using these  two implementations and
measured  updating time per  iteration \cite{Sato01b}.  It
turned out that both implementations run  twice as fast as
our naive implementation  and take  about 630 seconds  per
iteration while the  graphical  EM algorithm  takes  0.661
second per iteration, which  is still orders  of magnitude
faster than    the  former two.   Regrettably    a similar
comparison using  the {\em entire\/} EDR  corpus available
at the moment was abandoned  due to memory overflow during
parsing for the construction of support graphs.

Learning   experiments so  far   only   compared time  per
iteration which ignore   extra time  for search  (parsing)
required  by the graphical EM   algorithm.  So a  question
naturally   arises w.r.t.\  comparison   in terms of total
learning time.  Assuming  100 iterations for learning  ATR
corpus  however,  it is   estimated  that even considering
parsing time,  the  graphical  EM algorithm  combined with
MSLR parser runs orders of magnitude faster than the three
implementations    (ours, Kita's   and Johnson's)   of the
Inside-Outside algorithm  \cite{Sato01b}.  Of  course this
estimation does not directly apply   to the graphical   EM
algorithm     combined with OLDT     search,  as the  OLDT
interpreter will take more time than a parser and how much
more time is needed depends  on the implementaiton of OLDT
search.\footnote{
%
We  cannnot  answer this    question   right now  as   the
implementation of OLDT search is not completed.}
%
Conversely, however, we may be able to  take it as a rough
indication of how  far  our  approach,  the graphical   EM
algorithm combined  with  OLDT search via  support graphs,
can go in the domain of EM learning of PCFGs.

\subsection{Examing the Performance Gap}
\label{subsec:examine}

In the previous subsection, we compared the performance of
the  graphical  EM algorithm  against   the Inside-Outside
algorithm when   PCFGs are given,   using  two corpora and
three implementations of the Inside-Outside algorithm.  In
all experiments, the  graphical EM algorithm  considerably
outperformed the Inside-Outside algorithm despite the fact
that both have the same time complexity.  Now we look into
what causes such a performance gap.

Simply   put,   the  Inside-Outside   algorithm  is   slow
(primarily)  because it     lacks parsing.  Even   when  a
backbone CFG grammar is explicitly given, it does not take
any advantage of  the constraints imposed  by the grammar.
To see  it,  it might   help to   review how the    inside
probability $e(s,t,A)$, i.e.\ P(nonterminal $A$ spans from
$s$-th word to $t$-th word) $(s \leq t)$, is calculated by
the Inside-Outside algorithm for the given grammar.

\[
e(s,t,A) = 
    \sum_{B,C\; \mbox{s.t.}\; A \rightarrow BC\; \mbox{in the grammar}}
       \sum_{r=s}^{r=t-1} \mbox{P($A \rightarrow BC$)}e(s,r,B)e(r+1,t,C)
\]
%
Here P($A \rightarrow  BC$)  is a   probability associated
with a production rule $A \rightarrow  BC$.  Note that for
a fixed triplet   $(s,t,A)$, it  is  usual that  the  term
$\mbox{P($A\rightarrow     BC$)}e(s,r,B)e(r+1,t,C)$     is
non-zero only for a relatively small number of $(B,C,r)$'s
determined   from   successful  parses and   the  rest  of
combinations always give $0$ to the term.  Nonetheless the
Inside-Outside  algorithm attempts to  compute the term in
every iteration for all  possible combinations of $B$, $C$
and $r$ and this is repeated for every possible $(s,t,A)$,
resulting in   a  lot of redundancy.    The   same kind of
redundancy    occurs   in    the   computation of  outside
probability by the Inside-Outside algorithm.

The  graphical EM  algorithm is  free  of such  redundancy
because    it runs on   parse    trees  (a parse   forest)
represented by the support graph.\footnote{
%
We   emphasize   that      the   difference  between   the
Inside-Outside algorithm and the graphical EM algorithm is
solely computational efficiency,  and they converge to the
same parameter values when  starting from the same initial
values.     Linguistic    evaluations   of  the  estimated
parameters by the graphical EM algorithm are also reported
by \citeA{Sato01b}.  }
%
It must  be added, on  the other hand, that superiority in
learning speed of  the graphical EM algorithm  is realized
at the cost     of  space complexity  because   while  the
Inside-Outside  algorithm merely requires $O(NL^2)$  space
for  its array to  store  probabilities, the  graphical EM
algorithm needs $O(N^3L^3)$   space to store  the  support
graph where $N$  is the number  of nonterminals and $L$ is
the sentence length.  This  trade-off is understandable if
one notices that the graphical  EM algorithm applied to  a
PCFG   can be considered as    partial  evaluation of  the
Inside-Outside  algorithm  by   the    grammar (and    the
introduction    of  appropriate  data   structure  for the
output).

Finally we remark that the  use of parsing as a preprocess
for EM learning of PCFGs is not unique to the graphical EM
algorithm   \cite{Fujisaki89,Stolcke95}.  These approaches
however still seem to  contain redundancies compared  with
the  graphical      EM   algorithm.     For       instance
\citeA{Stolcke95} uses  an Earley  chart to compute inside
and   outside probability,   but   parses  are  implicitly
reconstructed  in each iteration  dynamically by combining
completed items.


\section{Related Work and Discussion}
\label{sec:work}

\subsection{Related Work}

The work  presented in this paper  is at the crossroads of
logic  programming and probability theory, and considering
an   enormous   body  of  work  done   in   these  fields,
incompleteness is unavoidable when reviewing related work.
Having said that,  we  look  at various attempts  made  to
integrate  probability with computational  logic or  logic
programming.\footnote{
%
We omit  literature  leaning strongly  toward  logic.  For
logic(s) concerning   uncertainty,  see  an    overview by
\citeA{Kyburg94}.  }
%
In reviewing, one  can immediately  notice there are   two
types of usage  of probability.  One type, {\em constraint
approach}, emphasizes    the   role of     probability  as
constraints and does  not necessarily  seek for a   unique
probability distribution over logical formulas.  The other
type, {\em distribution   approach}, explicitly defines  a
unique distribution by   model theoretical means  or proof
theoretical means,  to  compute  various  probabilities of
propositions.

A typical constraint approach is seen in the early work of
probabilistic logic  by \citeA{Nilsson86}.   His   central
problem,  ``probabilistic    entailment  problem'', is  to
compute the     upper and   lower  bound  of   probability
$\mbox{P}(\phi)$ of a target sentence $\phi$ in such a way
that the bounds are compatible with a given knowledge base
containing   logical   sentences (not  necessarily   logic
programs)   annotated    with    a  probability.     These
probabilities work as constraints on the possible range of
$\mbox{P}(\phi)$.    He used     the   linear  programming
technique to solve  this problem that  inevitably delimits
the applicability of his approach to finite domains.

Later \citeA{Lukasiewicz99} investigated the computational
complexity of   the probabilistic entailment  problem in a
slightly different  setting.  His knowledge base comprises
statements of  the form  $(H\mid G)[u_1,u_2]$ representing
$u_1\leq \mbox{P}(H\mid G)   \leq u_2$.   He  showed  that
inferring ``tight'' $u_1, u_2$ is  NP-hard in general, and
proposed a  tractable   class  of  knowledge  base  called
conditional constraint trees.

After  the  influential  work  of Nilsson, \citeA{Frish94}
introduced a deductive system for probabilistic logic that
remedies   ``drawbacks'' of  Nilsson's  approach,  that of
computational   intractability and  the  lack  of  a proof
system.  Their  system  deduces a probability range   of a
proposition   by  rules of probabilistic  inferences about
unconditional    and     conditional   probabilities.  For
instance, one of  the rules infers $P(\alpha  \mid \delta)
\in [0\; y]$ from $P(\alpha\vee\beta \mid \delta) \in [x\;
y]$ where $\alpha$,$\beta$  and $\delta$ are propositional
variables   and    $[x\;y]\;  (x\leq  y)$   designates   a
probability range.

Turning   to   logic  programming,   probabilistic   logic
programming        formalized   by   \citeA{Ng92}      and
\citeA{Dekhtyar97} was  also a constraint approach.  Their
program  is a set of annotated  clauses of the form $A:\mu
\leftarrow F_1:\mu_1,\ldots,F_n:\mu_n$  where $A$   is  an
atom,  $F_i$ ($1\leq  i\leq n$)  a basic  formula, i.e.\ a
conjunction or a disjunction of atoms, and $\mu_j$ ($0\leq
j\leq  n$)    a  sub-interval  of   $[0,1]$ indicating   a
probability   range.      A query   $\leftarrow  \exists\;
(F_1:\mu_1,\ldots,F_n:\mu_n)$ is  answered by an extension
of SLD refutation.  On formalization,  it is assumed  that
their language  contains only a  finite number of constant
and predicate symbols, and no function symbol is allowed.

A similar framework  was proposed  by \citeA{Lakshmanan94}
under   the  same   syntactic  restrictions (finitely many
constant and predicate symbols but no function symbols) in
a   different  uncertainty setting.    They used annotated
clauses    of   the    form   $A\stackrel{c}  {\leftarrow}
B_1,\ldots,B_n$ where $A$ and $B_i$ $(1\leq i \leq n)$ are
atoms            and          $c                  =\langle
[\alpha,\beta],[\gamma,\delta]\rangle$,     a   confidence
level,  represents   a   belief  interval $[\alpha,\beta]$
($0\leq\alpha\leq\beta\leq  1$)    and  a  doubt  interval
$[\gamma,\delta]$ ($0\leq\gamma\leq\delta\leq 1$),   which
an expert has in the clause.\\

As seen above,  defining a unique probability distribution
is of secondary or no  concern to the constraint approach.
This  is in sharp contrast   with Bayesian networks as the
whole discipline rests on  the ability of the  networks to
define        a     unique   probability      distribution
\cite{Pearl88,Castillo97}.   Researchers      in  Bayesian
networks have been seeking for  a  way of mixing  Bayesian
networks with a  logical representation to  increase their
inherently propositional expressive power.

\citeA{Breese92}  used logic programs to     automatically
build   a  Bayesian  network from  a   query.  In Breese's
approach, a  program is the  union  of a  definite  clause
program and a set of conditional  dependencies of the form
$\mbox{P}(P \mid Q_1\wedge\cdots\wedge Q_n)$ where $P$ and
$Q_i$s  are atoms.  Given a  query, a  Bayesian network is
constructed   dynamically  that   connects the query   and
relevant  atoms  in the program,  which  in turn defines a
{\em   local\/}   distribution  for  the  connected atoms.
Logical  variables  can appear  in  atoms  but no function
symbol is allowed.

\citeA{Ngo97} extended Breese's approach by  incorporating
a mechanism reflecting context.  They used a clause of the
form $\mbox{P}(A_0 \mid  A_1,\ldots,A_n)=\alpha \leftarrow
L_1,\dots,L_k$,   where    $A_i$'s   are   called  p-atoms
(probabilistic atoms)  whereas $L_j$'s  are context  atoms
disjoint from p-atoms,  and  computed by   another general
logic program (satisfying  certain restrictions).  Given a
query,  a set   of evidence and    context atoms, relevant
ground p-atoms are  identified by resolving context  atoms
away by SLDNF resolution, and a  local Bayesian network is
built to calculate  the probability  of  the query.   They
proved  the  soundness and   completeness  of their  query
evaluation procedure under the condition that programs are
acyclic\footnote {
%
The condition says  that every ground  atom $A$ must be
assigned    a  unique    integer    $n(A)$   such  that
$n(A)>n(B_1),\ldots,n(B_n)$ holds     for    any ground
instance  of   a   clause of the   form   $A \leftarrow
B_1,\ldots,B_n$.   Under this condition, when a program
includes  $p(X)\leftarrow  q(X,Y)$,     we cannot write
recursive    clauses  about $q$    such as  $q(X,[H|Y])
\leftarrow q(X,Y)$.}
%
and domains are finite.

Instead of defining  a local distribution for  each query,
\citeA{Poole93} defined a  {\em  global\/} distribution in
his   ``probabilistic    Horn abduction''.   His   program
consists  of   definite    clauses   and  {\em    disjoint
declarations\/}  of              the       form $\mbox{\tt
disjoint([$h_1$:$p_1$,...,$h_n$:$p_n$])}$  which specifies
a    probability  distribution     over    the  hypotheses
(abducibles)   $\{h_1,\ldots,h_n\}$.         He   assigned
probabilities to all  ground  atoms with the help  of  the
theory of logic programming,  and furthermore proved  that
Bayesian   networks are  representable  in  his framework.
Unlike previous approaches, his language contains function
symbols, but   the  acyclicity  condition imposed   on the
programs for his semantics to  be definable seems to be  a
severe  restriction.  Also,  probabilities are not defined
for quantified formulas.

\shortciteA{Bacchus96}  used  a    much    more   powerful
first-order  probabilistic language than clauses annotated
with probabilities.  Their language allows a statistically
quantified  term   such as  ${\parallel \phi(x)| \theta(x)
\parallel}_x$  to denote the   ratio  of individuals  in a
finite  domain  satisfying   $\phi(x)\wedge  \theta(x)$ to
those  satisfying  $\theta(x)$.  Assuming that every world
(interpretation  for their  language) is  equally likely,
they define the probability  of a sentence $\varphi$ under
the   given   knowledge   base     $K\!B$ as  the    limit
$\lim_{N\rightarrow         \infty}      \left(   \frac{\#
\mbox{worlds}^{\tau}_N(\varphi      \wedge   K\!B)}    {\#
\mbox{worlds}^{\tau}_N(K\!B)} \right)$       where     $\#
\mbox{worlds}^{\tau}_N(\chi)$   is the  number of possible
worlds containing $N$  individuals satisfying $\chi$,  and
$\tau$ parameters   used  in   judging     approximations.
Although   the limit does   not  necessarily exist and the
domain must be  finite, they showed  that their method can
cope with  difficulties arising from  ``direct inference''
and default reasoning.

In a more linguistic vein, \citeA[and others]{Muggleton96}
formulated SLPs (stochastic logic programs)  procedurally,
as an extension  of PCFGs to probabilistic logic programs.
So,     a      clause    $C$,       which      must     be
range-restricted,\footnote{
%
A syntactic property  that variables appearing in the head
also appear in the body of a clause. A unit clause must be
ground.}
%
is  annotated  with a  probability  $p$  like $p:C$.   The
probability  of a  goal $G$  is  the product of  such $p$s
appearing in its refutation  but with a  modification such
that  if a subgoal $g$ can  invoke  $n$ clauses, $p_i:C_i$
($1\leq    i \leq n$)    at    some refutation  step,  the
probability  of  choosing  $k$-th clause is  normalized to
$p_k/\sum_{i=1}^n p_i$.

More recently, \citeA{Cussens99,Cussens01}   enriched SLPs
by  introducing a special  class  of log-linear models for
SLD refutations w.r.t.   a  given  goal.  He for   example
considers all   possible  SLD   refutations for   the most
general  goal $\leftarrow\!  s(X)$ and defines probability
$\mbox{P}(R)$ of   a  refutation   $R$ as $\mbox{P}(R)   =
Z^{-1}\exp  \left(\sum_i \lambda_i \nu(R,i)\right)$.  Here
$\lambda_i$ is a number associated with a clause $C_i$ and
$\nu(R,i)$ is a  feature, i.e.\ the  number of occurrences
of $C_i$ in $R$.  $Z$  is the normalizing constant.  Then,
the  probability  assigned  to  $s(a)$   is   the sum   of
probabilities of refutation for $\leftarrow\!  s(a)$.


\subsection{Limitations and Potential Problems}

Approaches   described so far have   more  or less similar
limitations and   potential  problems.   Descriptive power
confined to finite domains is  the most common limitation,
which  is   due to  the  use   of  the linear  programming
technique  \cite{Nilsson86},   or due   to  the  syntactic
restrictions not allowing  for  infinitely many  constant,
function   or  predicate symbols \cite{Ng92,Lakshmanan94}.
Bayesian networks have the same limitation as well (only a
finite  number          of     random    variables     are
representable).\footnote{
%
However, RPMs (recursive  probability models)  proposed by
\citeA{Pfeffer00} as   an  extension of  Bayesian networks
allow for    infinitely many random  variables.   They are
organized  as  attributes  of classes   and a  probability
measure over attribute values is introduced.  }
%
Also there are  various semantic/syntactic restrictions on
logic programs.    For  instance the  acyclicity condition
imposed by  \citeA{Poole93} and \citeA{Ngo97} prevents the
unconditional use of clauses with local variables, and the
range-restrictedness   imposed  by \citeA{Muggleton96} and
\citeA{Cussens99}  excludes  programs  such  as the  usual
membership Prolog program.

There  is another  type  of  problem, the possibility   of
assigning   conflicting    probabilities   to    logically
equivalent  formulas.     In  SLPs,   $\mbox{P}(A)$    and
$\mbox{P}(A\wedge A)$ do  not necessarily coincide because
$A$ and   $A\wedge  A$  may  have different    refutations
\cite{Muggleton96,Cussens99,Cussens01}.    Consequently in
SLPs,    we would be in  trouble   if we naively interpret
$\mbox{P}(A)$ as the  probability   of $A$'s being   true.
Also   assigning probabilities    to arbitrary  quantified
formulas seems out of scope of both approaches to SLPs.

Last but not  least, there is a big  problem common to any
approach using  probabilities:  {\em where  do the numbers
come  from}?   Generally speaking,  if we   use $n$ binary
random variables in a   model, we have to determine  $2^n$
probabilities    to     completely specify  their    joint
distribution,     and  fulfilling  this   requirement with
reliable numbers quickly becomes  impossible as $n$ grows.
The  situation is even  worse when  there are unobservable
variables  in the  model   such as  possible  causes  of a
disease.  Apparently parameter learning from observed data
is  a   natural solution  to  this  problem, but parameter
learning of logic programs has not been well studied.

Distribution semantics  proposed by \citeA{Sato95}  was an
attempt to solve these problems    along the line of   the
global distribution approach.   It  defines a distribution
(probability measure) over the possible interpretations of
ground  atoms for an  arbitrary logic program in any first
order language and assigns consistent probabilities to all
closed formulas.  Also  distribution semantics enabled  us
to derive  an EM algorithm for   the parameter learning of
logic programs  for  the first time.   As  it  was a naive
algorithm   however,  dealing   with  large problems   was
difficult when there  are  exponentially many explanations
for  an   observation  like  HMMs.    We  believe that the
efficiency problem  is solved  to a  large extent   by the
graphical EM algorithm presented in this paper.

\subsection{EM Learning}

Since EM  learning is one  of the  central issues in  this
paper, we separately mention  work related to  EM learning
for symbolic  frameworks.  \citeA{Koller97a} used in their
approach to  KBMC  (knowledge-based model construction) EM
learning  to  estimate parameters  labeling clauses.  They
express   probabilistic   dependencies   among   events by
definite clauses   annotated with probabilities, similarly
to \citeS{Ngo97} approach,  and  locally build a  Bayesian
network  relevant to the  context  and evidence as well as
the   query.  Parameters are   learned by  applying to the
constructed   network the  specialized   EM algorithm  for
Bayesian networks \cite{Castillo97}.

Dealing with a  PCFG by a  statically constructed Bayesian
network  was   proposed  \citeA{Pynadath98},   and  it  is
possible to combine the EM  algorithm with their method to
estimate parameters    in the  PCFG.    Unfortunately, the
constructed  network  is  not singly  connected, and  time
complexity  of   probability computation  is   potentially
exponential in the length of an input sentence.

Closely related to our EM  learning is parameter  learning
of  log-linear models.  \citeA{Riezler98}  proposed the IM
algorithm   in  his approach  to  probabilistic constraint
programming.  The IM    algorithm is a  general  parameter
estimation algorithm from  incomplete data  for log-linear
models whose probability  function $\mbox{P}(x)$ takes the
form    $\mbox{P}(x) =   Z^{-1}\exp   \left(\sum_{i=1}^{n}
\lambda_i          \nu_i(x)    \right)p_0(x)$        where
$(\lambda_1,\ldots,\lambda_n)$   are  parameters  to    be
estimated,  $\nu_i(x)$ the  $i$-th feature of  an observed
object  $x$ and $Z$    the normalizing constant.   Since a
feature can  be any function of  $x$, the log-linear model
is highly  flexible and includes  our distribution $P_{\tt
msw}$ as a special case of $Z=1$.  There is a price to pay
however; the computational  cost   of $Z$. It requires   a
summation over exponentially many   terms.  To avoid   the
cost of  exact  computation, approximate computation  by a
Monte Carlo method is possible.   Whichever one may choose
however, learning    time increases  compared  to  the  EM
algorithm for $Z=1$.

The FAM (failure-adjusted maximization) algorithm proposed
by \citeA{Cussens01} is an EM algorithm applicable to pure
normalized  SLPs that may fail.  It   deals with a special
class of log-linear models  but is more efficient than the
IM algorithm.   Because the  statistical framework  of the
FAM  is   rather  different  from  distribution semantics,
comparison   with   the   graphical EM    algorithm  seems
difficult.

Being  slightly        tangential     to  EM     learning,
\shortciteA{Koller97b} developed  a    functional modeling
language defining a probability distribution over symbolic
structures in  which they  showed ``cashing'' of  computed
results  leads  to  efficient probability  computation  of
singly  connected   Bayesian  networks  and  PCFGs.  Their
cashing  corresponds   to  the  computation   of    inside
probability   in   the  Inside-Outside algorithm   and the
computation of outside probability is untouched.

\subsection{Future Directions}

Parameterized logic programs  are expected to be  a useful
modeling tool  for complex symbolic-statistical phenomena.
We  have   tried   various  types  of   modeling,  besides
stochastic  grammars and  Bayesian  networks, such  as the
modeling of   gene  inheritance  in   the  Kariera   tribe
\cite{White63}  where the rules of bi-lateral cross-cousin
marriage for four clans interact with the rules of genetic
inheritance    \cite{Sato98}.    The  model  was     quite
interdisciplinary,  but the flexibility  of combining {\tt
msw} atoms   by  means    of definite    clauses   greatly
facilitated the modeling process.

Although   satisfying the   five conditions          in
Section~\ref{sec:graphical-EM}
%
\begin{itemize}
\item the uniqueness condition
         (roughly, one cause yields one effect)
\item the finite   support condition
         (there  are  a   finite  number    of
          explanations for one observation)
\item the acyclic support condition
         (explanations must not be cyclic)
\item the t-exclusiveness condition
         (explanations must be mutually exclusive)
\item the independence condition
         (events in an explanation must be independent)
\end{itemize}
%
for the applicability of  the graphical EM algorithm seems
daunting, our modeling experiences so far tell us that the
modeling    principle in    Section~\ref{sec:graphical-EM}
effectively guides us to  successful modeling.  In return,
we can obtain a declarative model described compactly by a
high  level   language  whose  parameters  are efficiently
learnable by  the graphical EM algorithm  as shown  in the
preceding section.

One of the  future directions is  however to relax some of
the  applicability  conditions,  especially the uniqueness
condition that prohibits  a generative model from  failure
or from generating  multiple  observable events.  Although
we  pointed out in Section~\ref{subsec:revisited} that the
MAR  condition in Appendix~\ref{appendix:B} adapted to our
semantics   can    replace  the uniqueness  condition  and
validates the use of the graphical  EM algorithm even when
a complete data  does not uniquely determine the  observed
data just like the case of ``partially bracketed corpora''
\cite{Pereira92}, we feel the need  to do more research on
this topic.  Also investigating the role of the acyclicity
condition   seems      theoretically  interesting   as the
acyclicity  is often related  to   the learning  of  logic
programs \cite{Arimura97,Reddy98}.

In this paper we only  scratched the surface of individual
research fields such as HMMs, PCFGs and Bayesian networks.
Therefore, there remains  much to be done about clarifying
how experiences in each research   field are reflected  in
the  framework  of   parameterized  logic  programs.   For
example,   we need  to  clarify   the relationship between
symbolic  approaches  to   Bayesian  networks such  as SPI
\cite{Li94} and our approach.  Also   it is unclear how  a
compiled approach  using the  junction tree  algorithm for
Bayesian networks can be  incorporated into our  approach.
Aside   from  exact   methods,   approximate  methods   of
probability  computation specialized  for    parameterized
logic programs must also be developed.

There is also a direction of improving learning ability by
introducing  priors instead of ML  estimation to cope with
data sparseness.  The  introduction of basic distributions
that  make  probabilistic switches  correlated seems worth
trying in  the near future. It  is  also important to take
advantage of the logical  nature of our approach to handle
uncertainty.  For  example,  it    is already   shown   by
\citeA{Sato01c} that we can learn parameters from negative
examples such as ``the grass is {\em  not\/} wet'' but the
treatment   of  negative examples   in parameterized logic
programs is still in its infancy.

Concerning developing complex  statistical models based on
the  ``programs as  distributions''   scheme,   stochastic
natural  language  processing   which  exploits   semantic
information   seems         promising.     For   instance,
unification-based grammars such as HPSGs    \cite{Abney97}
may be a good target beyond PCFGs because they use feature
structures  logically describable,  and the  ambiguity  of
feature  values seems to  be  expressible by a probability
distribution.

Also building a mathematical basis for logic programs with
continuous random  variables   is a  challenging  research
topic.


\section{Conclusion}
\label{sec:conclusion}

We  have proposed   a logical/mathematical  framework  for
statistical  parameter  learning  of parameterized   logic
programs,  i.e.\ definite   clause programs     containing
probabilistic  facts   with  a parameterized   probability
distribution.   It extends  the traditional least Herbrand
model semantics in logic  programming to {\em distribution
semantics\/}, possible world  semantics with a probability
distribution  over  possible      worlds         (Herbrand
interpretations) which  is unconditionally  applicable  to
arbitrary logic programs   including ones for HMMs,  PCFGs
and Bayesian networks.

We   also have presented  a   new  EM algorithm, the  {\em
graphical            EM         algorithm\/}            in
Section~\ref{sec:graphical-EM}, which  learns  statistical
parameters from  observations for a class of parameterized
logic programs representing  a sequential decision process
in which  each decision is  exclusive and independent.  It
works on {\em support  graph\/}s,   a new data   structure
specifying a logical relationship between an observed goal
and   its   explanations,  and  estimates   parameters  by
computing  inside and  outside probability generalized for
logic programs.

The complexity  analysis   in Section~\ref{sec:complexity}
showed that when OLDT search, a complete tabled refutation
method for  logic  programs, is employed  for  the support
graph  construction and  table  access is  done  in $O(1)$
time, the graphical EM  algorithm, despite its generality,
has  the same time   complexity as existing EM algorithms,
i.e.\   the     Baum-Welch   algorithm  for     HMMs,  the
Inside-Outside algorithm for PCFGs  and the one for singly
connected Bayesian   networks  that   have  been developed
independently in  each research  field.  In addition,  for
pseudo  probabilistic context  sensitive grammars with $N$
nonterminals, we  showed  that the  graphical EM algorithm
runs in time $O(N^4 L^3)$ for a sentence of length $L$.

To   compare actual    performance  of the   graphical  EM
algorithm     against  the Inside-Outside  algorithm,   we
conducted     learning  experiments    with   PCFGs     in
Section~\ref{sec:experiments}  using two real corpora with
contrasting characters.   One  is   ATR  corpus containing
short  sentences for    which  the grammar  is not    much
ambiguous  (958   parses/sentence), and the  other  is EDR
corpus containing long  sentences for which the grammar is
rather  ambiguous  ($3.0\times  10^8$  at average sentence
length  20).  In  both cases,  the  graphical EM algorithm
outperformed the  Inside-Outside  algorithm by  orders  of
magnitude in terms of time  per iteration, which  suggests
the effectiveness of our approach   to EM learning by  the
graphical EM algorithm.

Since our  semantics is not limited  to finite  domains or
finitely many random variables but applicable to any logic
programs  of   arbitrary   complexity,  the   graphical EM
algorithm is  expected  to give  a  general  yet efficient
method  of parameter   learning   for  models of   complex
symbolic-statistical     phenomena governed  by  rules and
probabilities.


\acks{
%
The authors  wish to thank   three anonymous referees  for
their comments   and  suggestions.  Special thanks  go  to
Takashi Mori  and Shigeru Abe  for stimulating discussions
and learning   experiments, and  also  to  Tanaka-Tokunaga
Laboratory for kindly allowing them to use MSLR parser and
the linguistic data.}


\appendix

\section{Properties of $\pdb$}
\label{appendix:A}

In  this appendix,   we  list some   properties  of $\pdb$
defined by a parameterized logic program $\db  = F \cup R$
in a countable first-order language ${\cal L}$.\footnote{
%
For    definitions  of $\Omega_F$, $P_F$,  $M_{\db}(\nu)$,
$\Omega_{\db}$,  $\pdb$  and    others   used  below,  see
Section~\ref{sec:semantics}.}
%
First        of  all,   $\pdb$     assigns      consistent
probabilities\footnote{
%
By consistent, we  mean probabilities  assigned to logical
formulas respect the laws of  probability such as $0  \leq
P(A) \leq  1$, $P(\neg A) =1  -  P(A)$ and $P(A  \vee B) =
P(A)+P(B)-P(A\wedge B)$.  }
%
to every closed formula $\phi$ in ${\cal L}$ by
%
\[
\pdb(\phi)
 \ddefined
    \pdb(\{\omega \in \Omega_{\db} \mid \omega \models \phi\})
\]
%
while guaranteeing continuity in the sense that
\[
\begin{array}{rcl}
\lim_{n\rightarrow\infty}\pdb(\phi(t_1)\wedge\cdots\wedge \phi(t_n))
    & = & \pdb(\forall x \phi(x))\\
\lim_{n\rightarrow\infty}\pdb(\phi(t_1)\vee\cdots\vee \phi(t_n))
    & = & \pdb(\exists x \phi(x))
\end{array}
\]
%
where  $t_1,t_2,\ldots$ is an  enumeration of ground terms
in ${\cal L}$.\\

The    next proposition,  Proposition~\ref{proposition:1},
relates $\pdb$ to  the Herbrand  model.   To prove it,  we
need  some  terminology.   A {\em  factor\/}   is a closed
formula in prenex disjunctive  normal form $Q_1\cdots  Q_n
M$ where $Q_i$ $(1\leq  i\leq n)$ is either an existential
quantification or a  universal  quantification and  $M$  a
matrix.  The length of quantifications  $n$ is called  the
{\em  rank\/} of  the factor.  Define  $\Phi$  as a set of
formulas   made  out     of  factors,  conjunctions    and
disjunctions.   Associate   with  each  formula $\phi$  in
$\Phi$ a multi-set $r(\phi)$ of ranks by
%
\[
r(\phi) = \left\{
          \begin{array}{ccl}
             \emptyset & \mbox{if} &
                          \phi\;\mbox{is a factor with no quantification} \\
             \{ n \}   & \mbox{if} &
                          \phi\;\mbox{is a factor with rank}\; n \\
             r(\phi_1) \uplus r(\phi_2)
                       & \mbox{if} &
                          \phi=\phi_1 \vee \phi_2\;\mbox{or}\;
                          \phi=\phi_1 \wedge \phi_2.
         \end{array}
         \right.
\]
%
Here $\uplus$ stands for the union of two multi-sets.  For
instance $\{1,2,3\}\uplus \{2,3,4\}  =   \{1,2,2,3,3,4\}$.
We  use     the multi-set  ordering     in the    proof of
Proposition~\ref{proposition:1}   because   the      usual
induction on the complexity of formulas does not work.

\begin{lemma}
\label{lemma:1}
Let $\phi$ be a  boolean formula made  out of ground atoms
in $\cal   L$.   $\pdb(\phi)  =  P_F(\{\nu\in\Omega_F \mid
M_{\db}(\nu) \models \phi\})$.
\end{lemma}
%
(Proof)  We  have  only   to   prove the  lemma  about   a
conjunction of atoms of  the form $D^{x_1}_1 \wedge \cdots
\wedge D^{x_n}_n$ $(x_i\in\{0,1\}, 1\leq i\leq n)$.
%
\begin{eqnarray*}
\pdb(D^{x_1}_1 \wedge \cdots  \wedge D^{x_n}_n)
  & = & \pdb(\{\omega \in \Omega_{\db} \mid \omega \models
                D_1^{x_1}\wedge\cdots\wedge D_n^{x_n}\}) \\
  & = & \pdb(D_1=x_1,\ldots,D_n=x_n) \\
  & = & P_{F}(\{\nu \in \Omega_F \mid M_{\db}(\nu)\models
                D_1^{x_1}\wedge\cdots\wedge D_n^{x_n}\})
                                  \hspace{2em}\mbox{Q.E.D.}\\
\end{eqnarray*}
%
\begin{proposition}
\label{proposition:1}
Let $\phi$ be a closed formula in $\cal L$.  $\pdb(\phi)
= P_F(\{\nu\in\Omega_F \mid M_{\db}(\nu) \models \phi\})$.
\end{proposition}
%
(Proof) Recall that   a closed formula   has an equivalent
prenex disjunctive normal form that belongs to $\Phi$.  We
prove   the proposition for   formulas  in $\Phi$ by using
induction on the multi-set ordering over $\{ r(\phi) \mid
\phi\in \Phi\}$.  If  $r(\phi) =\emptyset$,  $\phi$ has no
quantification.     So  the  proposition   is  correct  by
Lemma~\ref{lemma:1}.  Suppose   otherwise.   Write $\phi =
G[Q_1Q_2\cdots Q_n    F]$   where  $Q_1Q_2\cdots  Q_n   F$
indicates   a  {\em single\/}  occurrence  of  a factor in
$G$.\footnote{
%
For an expression $E$, $E[\gamma]$ means that $\gamma$ may
occur  in   the    specified   positions  of    $E$.    If
$\gamma_1\vee\gamma_2$ in  $E[\gamma_1    \vee  \gamma_2]$
indicates    a    {\em      single\/}     occurrence    of
$\gamma_1\vee\gamma_2$ in a  positive boolean formula $E$,
$E[\gamma_1\vee   \gamma_2] = E[\gamma_1]\vee E[\gamma_2]$
holds.  }
%
We assume  $Q_1=\exists x$  ($Q_1=\forall x$ is  similarly
treated).  We also assume that bound variables are renamed
to  avoid name clash.  Then $G[\exists x Q_2\cdots Q_n F]$
is equivalent to $\exists  x G[Q_2\cdots Q_n F]$ in  light
of the  validity of  $(\exists x A)\wedge  B =  \exists  x
(A\wedge B)$ and $(\exists  x A)\vee B  = \exists x (A\vee
B)$ when $B$ contains no free $x$.
%
\begin{eqnarray*}
\pdb(\phi)
 & = & \pdb(G[Q_1Q_2\cdots Q_n  F]) \\
 & = & \pdb(\exists  x G[\,Q_2\cdots Q_n   F[x] \,]) \\
 & = & \lim_{k\rightarrow\infty}
               \pdb(G[\, Q_2\cdots Q_n F[t_1]]\vee\cdots\vee
                                         G[Q_2\cdots Q_n F[t_k] \,] )\\
 & = & \lim_{k\rightarrow\infty}
               \pdb(G [\, Q_2\cdots Q_n F[t_1]\vee\cdots\vee
                          Q_2\cdots Q_n F[t_k] \,] )\\
 & = & \lim_{k\rightarrow\infty}
          P_{F}(\{\nu\in\Omega_F \mid
                  M_{\db}(\nu) \models
                      G[\, Q_2\cdots Q_n F[t_1]\vee\cdots\vee
                          Q_2\cdots Q_n F[t_k] \,]\})\\
 &   & (\mbox{by induction hypothesis})\\
 & = & P_{F}(\{\nu\in\Omega_F \mid
                  M_{\db}(\nu) \models
                    \exists x G[Q_2\cdots Q_n F[x]] \})\\
 & = & P_{F}(\{\nu\in\Omega_F \mid
                  M_{\db}(\nu) \models \phi \}) \hspace{3em}\mbox{Q.E.D.}
\end{eqnarray*}

We next prove  a theorem on  the iff definition introduced
in Section~\ref{sec:graphical-EM}.  Distribution semantics
considers  the program  $\db  =  F  \cup R$  as  a set  of
infinitely many ground definite clauses such that $F$ is a
set of facts (with a  probability measure $P_F$) and $R$ a
set of rules, and  no clause head  in $R$ appears  in $F$.
Put
\[
   head(R) \ddefined \{ B \mid B\;\mbox{appears
           in}\;    R\;  \mbox{as  a  clause  head}\}.
\]
%
For    $B     \in    head(R)$,  let    $B\leftarrow   W_i$
$(i=1,2,\ldots)$ be an enumeration of clauses about $B$ in
$R$.   Define  $\ifff(B)$,  {\em the iff  (if-and-only-if)
form of rules about $B$} in $\db$\footnote{
%
This definition    is   different  from    the   usual one
\cite{Lloyd84,Doets94} as we  are  here talking at  ground
level.  $W_1 \vee W_2 \vee \cdots$ is true  if and only if
one of the disjuncts is true.}
%
by
\[
  \ifff(B) \ddefined
         B\leftrightarrow W_1 \vee  W_2 \vee \cdots
\]
%
Since  $M_{\db}(\nu)$  is   a least   Herbrand model,   the
following is obvious.
\begin{lemma}
\label{lemma:2}
For $B$ in  $head(R)$ and $\nu \in \Omega_F$, $M_{\db}(\nu)
\models \ifff(B)$.
\end{lemma}

Theorem~\ref{th:1}  below  is about $\ifff(B)$.  It states
that at general level, both  sides  of the iff  definition
$p(x)\leftrightarrow \exists y_1 (x=t_1 \wedge W_1) \vee
\cdots \vee \exists y_n  (x=t_n \wedge W_n)$ of $p(\cdot)$
coincide as random variables  whenever $x$ is instantiated
to a ground term.

\begin{theorem}
\label{th:1}

Let $\ifff(B) =  B \leftrightarrow W_1\vee W_2\vee \cdots$
be  the   iff   form  of rules  about     $B \in head(R)$.
$\pdb(\ifff(B))=1$   and    $\pdb(B)= \pdb(W_1\vee
W_2\vee\cdots)$.
\end{theorem}
%
(Proof)
\begin{eqnarray*}
\pdb(\ifff(B))
& = &
      \pdb(\{ \omega\in\Omega_{\db} \mid
         \omega\models B \wedge (W_1\vee W_2\vee \cdots) \}) \\
&   & \hspace{3em}
        + \pdb(\{ \omega\in\Omega_{\db} \mid
         \omega\models \neg B \wedge \neg (W_1\vee W_2\vee \cdots) \}) \\
& = &
      \lim_{k\rightarrow \infty}
        \pdb(\{ \omega\in\Omega_{\db} \mid
         \omega\models B \wedge \bigvee^{k}_{i=1} W_i \}) \\
&   & \hspace{3em}
        + \lim_{k\rightarrow \infty}
             \pdb(\{ \omega\in\Omega_{\db} \mid
             \omega\models \neg B \wedge \neg\bigvee^{k}_{i=1} W_i \}) \\
& = &
      \lim_{k\rightarrow \infty}
        P_{F}(\{ \nu\in\Omega_{F} \mid
         M_{\db}(\nu)\models B \wedge \bigvee^{k}_{i=1} W_i \}) \\
&   & \hspace{3em}
      + \lim_{k\rightarrow \infty}
           P_{F}(\{ \nu\in\Omega_{F} \mid
            M_{\db}(\nu)\models \neg B \wedge \neg\bigvee^{k}_{i=1} W_i \})\\
&   & (\mbox{Lemma~\ref{lemma:1}}) \\
& = &
       P_{F}(\{ \nu\in\Omega_{F} \mid M_{\db}(\nu)\models \ifff(B) \})
              \\
& = & P_{F}(\Omega_{F})  \;\; (\mbox{Lemma~\ref{lemma:2}})\\
& = & 1
\end{eqnarray*}
It follows from $\pdb(\ifff(B)) = 1$ that
\[
\pdb(B) =  \pdb(B\wedge\ifff(B))
		=  \pdb(W_1\vee W_2\vee \cdots).
             \hspace{4em}\mbox{Q.E.D.}
\]

We  then prove  a   proposition   useful   in  probability
computation.  Let  $\psi_{\db}(B)$ be the  support set for
an   atom $B$ introduced in Section~\ref{sec:graphical-EM}
(it is the set   of all explanations   for $B$).   In  the
sequel, $B$ is    a ground atom.  Write   $\psi_{\db}(B) =
\{S_1,S_2,\ldots \}$ and $\bigvee \psi_{\db}(B) = S_1 \vee
S_2 \vee\cdots$\footnote{
%
For a set $K = \{E_1,E_2, \ldots \}$ of formulas, $\bigvee
K$ denotes a (-n infinite) disjunction $E_1 \vee E_2 \vee
\cdots$}
%
Define a set $\Lambda_B$ by
\[
\Lambda_B
      \ddefined
       \{ \omega \in \Omega_{\db} \mid
           \omega \models  B  \leftrightarrow  \bigvee \psi_{\db}(B) \}.
\]

\begin{proposition}
\label{proposition:2}

For every $B\in head(R)$,
$\pdb(\Lambda_B) = 1$ and $\pdb(B) =\pdb(\bigvee \psi_{\db}(B))$.

\end{proposition}
%
(Proof) We first  prove $\pdb(\Lambda_B) =1$ but the proof
exactly parallels  that of  Theorem~\ref{th:1} except that
$W_1\vee W_2\vee\cdots$ is replaced by  $S_1 \vee S_2 \vee
\cdots$  using the fact  that  $B \leftrightarrow S_1 \vee
S_2 \vee \cdots$ is true in every  least Herbrand model of
the form $M_{\db}(\nu)$.  Then from $\pdb(\Lambda_B) = 1$,
we have
\begin{eqnarray*}
\pdb(B)
 & = & \pdb(B \wedge  (B\leftrightarrow \bigvee \psi_{\db}(B))) \\
 & = & \pdb(\bigvee \psi_{\db}(B)). \hspace{4em}\mbox{Q.E.D.}
\end{eqnarray*}

Finally,  we   show that   distribution  semantics   is  a
probabilistic extension  of the traditional least Herbrand
model   semantics  in    logic programming    by   proving
Theorem~\ref{th:2}. It says  that the probability mass  is
distributed  exclusively  over possible    least  Herbrand
models.

Define $\Lambda$ as    the set of  least  Herbrand  models
generated by fixing $R$ and varying a subset of $F$ in the
program $\db=F\cup R$. In symbols,

\[
\Lambda \ddefined
   \{ \omega \in \Omega_{\db} \mid \omega = M_{\db}(\nu)
                   \;\mbox{for some}\; \nu \in \Omega_F
   \}.
\]
%
Note   that   as  $\Lambda$  is     merely  a  subset   of
$\Omega_{\db}$,  we  cannot conclude  $\pdb(\Lambda)=1$  a
priori, but the  next theorem,  Theorem~\ref{th:2}, states
$\pdb(\Lambda)     =1$,   i.e.   distribution    semantics
distributes   the  probability   mass   exclusively   over
$\Lambda$, i.e. possible least Herbrand models.

To  prove   the   theorem, we    need  some  preparations.
Recalling  that  atoms  outside  $head(R)\cup  F$ have  no
chance of being proved from $\db$, we introduce
%
\[
\Lambda' \ddefined
  \{ \omega \in\Omega_{\db} \mid
     \omega \models \neg  D\;
     \mbox{for every ground atom}\; D \not\in head(R)\cup F\}.
\]
%
For a  Herbrand interpretation $\omega  \in \Omega_{\db}$,
$\omega|_F$   $(\in \Omega_F)$  is   the  restriction   of
$\omega$ to those atoms in $F$.

\begin{lemma}
\label{lemma:3}
Let $\omega\in \Omega_{\db}$ be a Herbrand interpretation.\\
%
$\omega = M_{\db}(\nu)\; \mbox{for some}\; \nu\in \Omega_F$
iff
$\omega \in \Lambda'$ and
$\omega   \models  B \leftrightarrow \bigvee \psi_{\db}(B)$
for every $B \in head(R)$.
\end{lemma}
%
(Proof) Only-if part is immediate from the property of the
least Herbrand   model.  For   if-part, suppose   $\omega$
satisfies  the  right hand side.   We  show that $\omega =
M_{\db}(\omega|_F)$.  As $\omega$ and $M_{\db}(\omega|_F)$
coincide  w.r.t.\ atoms not in  $head(R)$, it is enough to
prove that  they also give the same  truth values to atoms
in $head(R)$.
%
Take $B\in head(R)$ and  write  $ \bigvee \psi_{\db}(B)  =
S_1\vee  S_2   \vee  \cdots$   Suppose  $\omega \models  B
\leftrightarrow S_1 \vee S_2  \vee \cdots$ Then if $\omega
\models  B$, we have $\omega   \models S_j$ for some  $j$,
thereby      $\omega|_F   \models     S_j$,    and   hence
$M_{\db}(\omega|_F)  \models  S_j$,        which   implies
$M_{\db}(\omega|_F) \models B$.  Otherwise $\omega \models
\neg B$.  So $\omega \models \neg S_j$  for every $j$.  It
follows that $M_{\db}(\omega|_F)  \models \neg B$.   Since
$B$      is arbitrary,  we    conclude   that $\omega$ and
$M_{\db}(\omega|_F)$ agree on the truth values assigned to
atoms in $head(R)$ as well. \hspace{3em}Q.E.D.\\

\begin{theorem} $\pdb(\Lambda) = 1$.
\label{th:2}
\end{theorem}
(Proof) From Lemma~\ref{lemma:3}, we have
%
\begin{eqnarray*}
\Lambda & = & 
   \{ \omega \in \Omega_{\db} \mid \omega = M_{\db}(\nu)\;
                        \mbox{for some}\; \nu \in \Omega_F \} \\
        & = & \Lambda' \cap \bigcap_{B\in head(R)} \Lambda_B.
\end{eqnarray*}
%
$\pdb(\Lambda_B) = 1$ by  Proposition~\ref{proposition:2}.
To prove $\pdb(\Lambda')  = 1$, let $D_1,D_2,\ldots$ be an
enumeration  of  atoms not belonging   to $head(R)\cup F$.
They are not  provable  from $\db =   F\cup R$, and  hence
false   in  every   least   Herbrand model  $M_{\db}(\nu)$
$(\nu\in\Omega_F)$. So
%
\begin{eqnarray*}
\pdb(\Lambda')
& = &  \lim_{m\rightarrow \infty}
         \pdb(\{ \omega\in\Omega_{\db} \mid
         \omega\models \neg D_1\wedge\cdots\wedge \neg D_m \}) \\
& = &  \lim_{m\rightarrow \infty}
         P_F(\{ \nu\in\Omega_F \mid
         M_{\db}(\nu) \models \neg D_1\wedge\cdots\wedge \neg D_m \}) \\
& = &  P_F(\Omega_F) =  1.
\end{eqnarray*}

Since  a countable   conjunction  of measurable   sets  of
probability measure one has  also probability measure one,
it follows  from   $\pdb(\Lambda_B)=1$  for every  $B  \in
head(R)$           and       $\pdb(\Lambda')=1$       that
$\pdb(\Lambda)=1$. \hspace{4em}Q.E.D. \\


\section{The MAR (missing at random) Condition}
\label{appendix:B}

In the original  formulation    of the EM  algorithm    by
\citeA{Dempster77},  it  is assumed   that there  exists a
many-to-one mapping $y = \chi(x)$ from a complete data $x$
to  an incomplete  (observed)  data $y$.   In  the case of
parsing, $x$ is a parse tree and $y$ is the input sentence
and    $x$ uniquely determines $y$.     In this paper, the
uniqueness  condition   ensures  the existence  of  such a
many-to-one mapping from explanations to observations.  We
however sometimes face a  situation where there is no such
many-to-one mapping  from complete data to incomplete data
but nonetheless we wish to apply the EM algorithm.

This  dilemma can   be solved  by the  introduction   of a
missing-data   mechanism   which makes  a   complete  data
incomplete.    The   missing-data  mechanism,  $m$,  has a
distribution  $g_{\phi}(m \mid x)$ parameterized by $\phi$
and $y$, the observed data, is described as $y=\chi_m(x)$.
It says  $x$  becomes    incomplete  $y$ by   $m$.     The
correspondence between $x$ and $y$, i.e.  $\{ \langle x, y
\rangle  \mid  \exists m   (y  = \chi_m(x))  \}$ naturally
becomes many-to-many.

\citeA{Rubin76}  derived two  conditions on    $g_{\phi}$
(data  are  missing at random  and   data are  observed at
random)   collectively  called  the  {\em  MAR (missing at
random)   condition}, and  showed  that   if  we  assume a
missing-data  mechanism   behind   our  observations  that
satisfies the MAR condition, we may estimate parameters of
the   distribution  over $x$  by    simply applying the EM
algorithm to $y$, the observed data.

We adapt the MAR condition to parameterized logic programs
as follows.  We   keep a generative  model satisfying  the
uniqueness condition that  outputs goals $G$ such as parse
trees.   We  further   extend the model   by  additionally
inserting a missing-data mechanism $m$ between $G$ and our
observation $O$   like   $O=\chi_m(G)$ and   assume    $m$
satisfies the MAR  condition.  Then the extended model has
a  many-to-many correspondence   between  explanations and
observations, and    generates  non-exclusive observations
such that  $P(O\wedge O')>0$ $(O  \neq  O')$, which causes
$\sum_{O}  P(O)\geq 1$ where $P(O)  = \sum_{G: \exists m\;
O=\chi_m(G)}  \pdb(G)$.  Thanks    to  the  MAR  condition
however, we are still allowed to apply the EM algorithm to
such non-exclusive observations.  Put it differently, even
if the uniqueness condition is seemingly destroyed, the EM
algorithm  is applicable just  by (imaginarily) assuming a
missing-data mechanism satisfying the MAR condition.


\bibliography{sato01a}
\bibliographystyle{theapa}

\end{document}
