%\VignetteIndexEntry{grain: Introduction}
%\VignetteDepends{gRbase}
%\VignetteKeyword{Bayesian networks}
%\VignetteKeyword{Graphical Models}
%\VignetteKeyword{Probabilistic networks}
%\VignetteEngine{knitr::knitr} 

\documentclass[10pt]{article}
\usepackage{boxedminipage,color,a4wide,url}
\usepackage[utf8]{inputenc}

\usepackage[authoryear,round]{natbib}
\bibliographystyle{plainnat}

<<include=FALSE,echo=FALSE,warning=FALSE>>=
library(knitr)
dir.create("figures")
opts_chunk$set(fig.height=2.5,
               fig.path='figures/grain-',
               warning=FALSE, message=FALSE
)
options("prompt"="> ","width"=85)
@

\usepackage{etoolbox} 
\makeatletter 
\preto{\@verbatim}{\topsep=0pt \partopsep=-25pt } 
\makeatother

\usepackage{alltt}
\AtBeginEnvironment{alltt}{\setlength{\topsep}{-25pt}}



\def\pkg#1{{\bf #1}}
\def\grbn{{\bf gRain}}
\def\grain{\texttt{grain}}
\def\code#1{{\texttt{#1}}}
\def\R{\texttt{R}}
\def\bn{Bayesian network}

<<echo=FALSE>>=
require(gRain)
prettyVersion <- packageDescription("gRain")$Version
prettyDate <- format(Sys.Date())
@

\author{S{\o}ren H{\o}jsgaard\\Aalborg University, Denmark}
\title{Bayesian networks in R with the \pkg{gRain} package}
\date{\pkg{gRain} version \Sexpr{prettyVersion} as of \Sexpr{prettyDate}}

\begin{document}
\maketitle
\tableofcontents
\parindent0pt\parskip5pt

\section{Introduction}

The \grbn\ package implements Bayesian Networks (hereafter often
abbreviated BNs). The name \grbn\ is an acronym for [gra]phical
[i]ndependence [n]etworks. The main reference for \grbn\  is
\cite{hoj:12}, see also `citation("gRain")`.


Moreover, \cite{hoj:edw:lau:12} gives a broad treatment of graphical
models (including Bayesian networks) More information about the
package, other graphical modelling packages and development versions
is available from

\begin{quote}
\url{http://people.math.aau.dk/~sorenh/software/gR}
\end{quote}

%% \section{Two worked examples}
%% \label{sec:two-worked-examples}

%% \subsection{Wet grass}
%% \label{sec:wet-grass}

%% The `wet grass` example is motivated by the following narrative (taken from 
%% \url{https://en.wikipedia.org/wiki/Bayesian_network})
%% \begin{quote}
%%   Two events can cause grass to be wet: an active sprinkler or
%%   rain. Rain has a direct effect on the use of the sprinkler (namely
%%   that when it rains, the sprinkler usually is not active).
%% \end{quote}

%% <<echo=F, results='hide'>>=
%% p.R    <- cptable(~R, values=c(2, 8), levels=yn)
%% p.S_R  <- cptable(~S:R, values=c(1, 99, 4, 6), levels=yn)
%% p.G_SR <- cptable(~G:S:R, values=c(99, 1, 8, 2, 9, 1, 0, 1), levels=yn)
%% grass_bn <- compileCPT(p.R, p.S_R, p.G_SR)  %>% grain
%% @ 

%% <<chest-grass, fig.height=2, echo=F, fig.cap="Wet graph example; taken from Wikipedia.">>=
%% plot(grass_bn)
%% @ %def


\section{Example: Chest clinic}
\label{sec:chest-clinic}
\label{sec:chest}

% <<echo=F, results='hide'>>=
% yn <- c("yes","no")
% a    <- cptable(~asia, values=c(1,99),levels=yn)
% t.a  <- cptable(~tub|asia, values=c(5,95,1,99),levels=yn)
% s    <- cptable(~smoke, values=c(5,5), levels=yn)
% l.s  <- cptable(~lung|smoke, values=c(1,9,1,99), levels=yn)
% b.s  <- cptable(~bronc|smoke, values=c(6,4,3,7), levels=yn)
% e.lt <- cptable(~either|lung:tub,values=c(1,0,1,0,1,0,0,1),levels=yn)
% x.e  <- cptable(~xray|either, values=c(98,2,5,95), levels=yn)
% d.be <- cptable(~dysp|bronc:either, values=c(9,1,7,3,8,2,1,9), levels=yn)
% plist <- compileCPT(list(a, t.a, s, l.s, b.s, e.lt, x.e, d.be))
% plist
% chest_bn <- grain(plist)
% chest_bn
% @ %def

<<echo=F, results='hide'>>=
yn <- c("yes","no")
a    <- cpt(~asia, values=c(1,99),levels=yn)
t.a  <- cpt(~tub|asia, values=c(5,95,1,99),levels=yn)
s    <- cpt(~smoke, values=c(5,5), levels=yn)
l.s  <- cpt(~lung|smoke, values=c(1,9,1,99), levels=yn)
b.s  <- cpt(~bronc|smoke, values=c(6,4,3,7), levels=yn)
e.lt <- cpt(~either|lung:tub,values=c(1,0,1,0,1,0,0,1),levels=yn)
x.e  <- cpt(~xray|either, values=c(98,2,5,95), levels=yn)
d.be <- cpt(~dysp|bronc:either, values=c(9,1,7,3,8,2,1,9), levels=yn)
plist <- compileCPT(list(a, t.a, s, l.s, b.s, e.lt, x.e, d.be))
plist
chest_bn <- grain(plist, compile=FALSE)
chest_bn
@ %def










This section reviews the chest clinic example of \cite{lau/spieg:88}
(illustrated in Figure~\ref{fig:chest-LS}) and shows one way of
specifying the model in \grbn{}.  \cite{lau/spieg:88} motivate the
chest clinic example with the following narrative:

\begin{quote}
  ``Shortness--of--breath (dyspnoea) may be due to tuberculosis, lung
  cancer or bronchitis, or none of them, or more than one of them. A
  recent visit to Asia increases the chances of tuberculosis, while
  smoking is known to be a risk factor for both lung cancer and
  bronchitis. The results of a single chest X--ray do not discriminate
  between lung cancer and tuberculosis, as neither does the presence or
  absence of dyspnoea.''
\end{quote}


<<chest-LS, echo=F, fig.cap="Chest clinic example from Lauritzen and Spiegelhalter (1988).">>=
plot(chest_bn)
@ %def

\subsection{Building a network}

The description above involves the following binary variables:
$\alpha=\mbox{asia}$,
$\sigma=\mbox{smoker}$,
$\tau=\mbox{tuberculosis}$,
$\lambda=\mbox{lung cancer}$,
$\beta=\mbox{bronchitis}$,
$\epsilon=\mbox{either tuberculosis or lung cancer}$,
$\delta=\mbox{dyspnoea}$ and
$\xi=\mbox{xray}$. 
Each variable is binary and can take the values ``yes'' and ``no'':
Note that $\epsilon$ is a logical variable which is
true (yes) if either $\tau$ or $\lambda$ are true (yes) and false (no) otherwise.
The connection between the variables is displayed by the DAG (directed acyclic graph) in
Figure~\ref{fig:chest-LS}.

A joint probability density factorizing according to a DAG with nodes
$V$ can be constructed as follows: Each node $v\in V$ has a set $pa(v)$ of parents and each node
$v\in V$ has a finite set of states. A joint distribution
over the variables $V$ can be given as
\begin{equation}
  \label{eq:dagfact1}
  p(V) = \prod_{v\in V} p(v|pa(v))
\end{equation}
where $p(v|pa(v))$ is a function defined on $(v,pa(v))$. This function
satisfies that $\sum_{v^*} p(v=v^*|pa(v))=1$, i.e.\ that
for each configuration of the parents $pa(v)$, the sum
over the levels of $v$ equals one. Hence $p(v|pa(v))$ becomes the
conditional distribution of $v$ given $pa(v)$.
In practice $p(v|pa(v))$ is specified as a table called a conditional
probability table or a CPT for short.
Thus, a Bayesian network can be regarded as a complex stochastic model built up by
putting together simple components (conditional probability
distributions).
A joint probability density for all eight variables in
Figure~\ref{fig:chest-LS}
can be constructed as 
\begin{equation}
  \label{eq:chestfact1}
  p(V) =
  p(\alpha)p(\sigma)p(\tau|\alpha)p(\lambda|\sigma)p(\beta|\sigma)p(\epsilon|\tau,\lambda)
  p(\delta|\epsilon, \beta)p(\xi|\epsilon).
\end{equation}



\subsection{Queries to networks}
\label{sec:xxx}

Suppose we are given the evidence (sometimes also called ``finding'')
that a set of variables $E\subset V$
have a specific value $e^*$.
With this evidence, we are often interested in the conditional
distribution $p(v|E=e^*)$
for some of the variables $v \in V \setminus E$
or in $p(U|E=e^*)$
for a set $U\subset V \setminus E$. Interest might also be in
calculating the probability of a specific event, e.g.\ the probability
of seeing a specific evidence, i.e.\ $p(E=e^*)$.
Other types of evidence (called soft evidence, virtual evidence or likelihood evidence) are discussed in
Section~\ref{sec:hard-virt-likel}.

For
example that a person has recently visited Asia and suffers from
dyspnoea, i.e.\ $\alpha=\mbox{yes}$ and $\delta=\mbox{yes}$.
In the chest clinic example, interest might be in $p(\lambda|e^*)$, $p(\tau|e^*)$
and  $p(\beta|e^*)$, or possibly in the joint (conditional) distribution
$p(\lambda,\tau,\beta|e^*)$.


\section{A one--minute version of  \grbn{}}
\label{sec:oneminute}

\subsection{Specifying a network}
\label{sec:specifying-network}

A simple way of  specifying the model for the chest clinic
example is as follows.

\begin{enumerate}
\item Specify conditional probability tables (with values as given in
  \cite{lau/spieg:88}) (there are other ways of specifying conditional
  probability tables, see the package documentation):

% <<echo=FALSE, results='hide'>>=
% yn <- c("yes", "no")
% a    <- cptable(~asia, values=c(1, 99), levels=yn)
% t.a  <- cptable(~tub|asia, values=c(5, 95, 1, 99), levels=yn)
% s    <- cptable(~smoke, values=c(5, 5), levels=yn)
% l.s  <- cptable(~lung|smoke, values=c(1, 9, 1, 99), levels=yn)
% b.s  <- cptable(~bronc|smoke, values=c(6, 4, 3, 7), levels=yn)
% e.lt <- cptable(~either|lung:tub, values=c(1, 0, 1, 0, 1, 0, 0, 1), levels=yn)
% x.e  <- cptable(~xray|either, values=c(98, 2, 5, 95), levels=yn)
% d.be <- cptable(~dysp|bronc:either, values=c(9, 1, 7, 3, 8, 2, 1, 9), levels=yn)
% @ %def




<<>>=
yn <- c("yes","no")
a    <- cpt(~asia, values=c(1, 99),levels=yn)
t.a  <- cpt(~tub|asia, values=c(5, 95, 1, 99),levels=yn)
s    <- cpt(~smoke, values=c(5, 5), levels=yn)
l.s  <- cpt(~lung|smoke, values=c(1, 9, 1, 99), levels=yn)
b.s  <- cpt(~bronc|smoke, values=c(6, 4, 3, 7), levels=yn)
e.lt <- cpt(~either|lung:tub,values=c(1, 0, 1, 0, 1, 0, 0, 1),levels=yn)
x.e  <- cpt(~xray|either, values=c(98, 2, 5, 95), levels=yn)
d.be <- cpt(~dysp|bronc:either, values=c(9, 1, 7, 3, 8, 2, 1, 9), levels=yn)
@ 

\item Compile list of conditional probability tables.

<<>>=
chest_cpt <- compileCPT(a, t.a, s, l.s, b.s, e.lt, x.e, d.be)
chest_cpt
@ %def

The components are arrays, but coercion into dataframes sometimes makes it easier to digest the components.
<<>>=
chest_cpt$tub
chest_cpt$tub  |> as.data.frame.table()
@ 

Notice: \code{either} is a logical node:
<<>>=
chest_cpt$either  |> as.data.frame.table()
@ 


\item Create the network:
<< >>= 
chest_bn <- grain(chest_cpt)
chest_bn
@

Default is that the network is compiled at creation time, but if one
chooses not to do so, compilation can be done with:

<< >>= 
chest_bn <- compile(chest_bn)
@

\end{enumerate}

\subsection{Setting evidence and querying a network}
\label{sec:querying-network}

\begin{enumerate}
\item A network can be queried to give marginal probabilities for
  each of a set of nodes (the default) or the joint probability for a
  set of nodes.\footnote{A third type of output exists, see package
    documentation for details.}  Notice that \code{querygrain()} can
  be abbreviated \code{qgrain()}.

<<>>=
querygrain(chest_bn, nodes=c("lung", "bronc"), type="marginal")
@ %def


\item Likewise, a joint distribution can be obtained (represented as a
  multi dimensional array):

<<>>=
querygrain(chest_bn, nodes=c("lung", "bronc"), type="joint")
@ %def

\item Evidence can be entered in two different ways:\footnote{Alternative forms
    exist; see package documentation for details.}

<<>>=
chest_ev  <- setEvidence(chest_bn,
                          evidence=list(asia="yes", dysp="yes"))
chest_ev  <- setEvidence(chest_bn,
                          nodes=c("asia", "dysp"), states=c("yes", "yes"))

## Also: modify object with
## evidence(chest_bn) <- list(asia="yes", dysp="yes")
@ %def

\item The evidence is a list and can conveniently be displayed as a dataframe:

<<>>=
getEvidence(chest_ev) |> as.data.frame()
@ 

\item The network can be queried again:
<<>>=
querygrain(chest_ev, nodes=c("lung", "bronc"))
@ %def


\item The probability of observing this evidence under the model is
<<>>=
pEvidence(chest_ev)
@ %def

\item The probability of an evidence can be found with only
  propagation towards the root of a junction tree. This saves about
  half the computational effort of propagation. However, notice that the
  network is not changed by this operation, so if the network is
  subsequently queried, there is no gain in computing time.

<<>>=
pEvidence(chest_bn, evidence=list(asia="yes", dysp="yes"))
@ 


\end{enumerate}


\section{Hints and shortcuts}
\label{sec:small-shortcuts}

\begin{enumerate}
\item An alternative way of specifying a network is by first defining CPTs and then entering values afterwards. Programmatically, this can be done as:

<<>>=
yn <- c("yes","no")
flist <- c(
    ~asia, ~tub|asia, ~smoke, ~lung|smoke, ~bronc|smoke, ~either|lung:tub,
    ~xray|either, ~dysp|bronc:either
)
## or
flist <- list("asia", c("tub", "asia"), "smoke", c("lung", "smoke"),
              c("bronc", "smoke"), c("either", "tub", "lung"),
              c("xray", "either"), c("dysp", "bronc", "either"))

chest_cpt2 <- lapply(flist, function(f){
    cpt(f, levels=yn)
})

bn2 <- compileCPT(chest_cpt2) |> grain()

lst2 <- list(asia=c(1, 99),
            tub=c(5, 95, 1, 99),
            smoke=c(5, 5),
            lung=c(1, 9, 1, 99),            
            bronc=c(6, 4, 3, 7),
            either=c(1, 0, 1, 0, 1, 0, 0, 1),
            xray=c(98, 2, 5, 95),
            dysp=c(9, 1, 7, 3, 8, 2, 1, 9))

bn2 <- replaceCPT(bn2, lst2)
@ 

\item Consider querying a network where focus is on 
  marginal distributions (the default). If all variables have the same
  levels (as the case is here), the output can be coerced to a
  dataframe:

<<>>=
querygrain(chest_bn, nodes=c("lung", "bronc"), simplify = TRUE)
@ 

In the more general case the output can be coerced to a list of dataframes
<<>>=
querygrain(chest_bn, nodes=c("lung", "bronc"), result="data.frame")
@ 

\item A typical use of \bn\ involves setting evidence and then
  querying the network afterwards. This can all be done in one call of
  \code{querygrain} (notice that this does not alter the network object):

<<>>= 
querygrain(chest_bn,
           evidence=list(asia="yes", dysp="yes"),
           nodes=c("lung", "bronc"), simplify = TRUE)
@

\item Evidence can be also be given as a vector of weights. 

<<>>= 
querygrain(chest_bn,
           evidence=list(asia=c(1,0), dysp=c(1,0)),
           nodes=c("lung", "bronc"), simplify = TRUE)
@

The weights must be non-negative but need not sum to one. This is
important in connection with soft evidence (also called likelihood
evidence), see Section~\ref{sec:hard-soft}. Above, the weights could
also have been set as \code{c(.1, 0)}. The important part is that the
zero excludes certain states as being impossible.


\item Nodes on which evidence is given are not reported unless \code{exclude=FALSE}

<<>>= 
querygrain(chest_bn,
           evidence=list(asia=c(1, 0), dysp=c(1, 0)),
           nodes=c("lung", "bronc", "asia", "dysp"),
           exclude=FALSE, simplify = TRUE)
@


\item If \code{nodes} are not specified, queries for all nodes without evidence are reported.
<<>>= 
querygrain(chest_bn,
           evidence=list(asia="yes", dysp="yes"),
           simplify = TRUE)
@

If \code{nodes} are not specified and \code{exclude=FALSE}, then queries for all nodes are reported.

<<>>= 
querygrain(chest_bn,
           evidence=list(asia="yes", dysp="yes"),
           exclude = FALSE, simplify = TRUE)
@


\end{enumerate}



\section{Conditioning on evidence with zero probability}
\label{sec:zero-probabilities}

Consider setting the evidence
<<>>=
chest_bn3 <- setEvidence(chest_bn, evidence=list(either="no", tub="yes"))
@ %def

Under the model, this specific evidence has zero probability:
\verb|either| is true if \verb|tub| is true or \verb|lung| is true (or
both). Hence the specific evidence is impossible and therefore, all
conditional probabilities are (under the model) undefined:

<<>>=
pEvidence(chest_bn3)
querygrain(chest_bn3, nodes=c("lung", "bronc"), type="joint")
@ %def




\section{Hard  and virtual (likelihood) evidence}
\label{sec:hard-virt-likel}
\label{sec:hard-soft}

Below we describe  how to work with virtual evidence (also known
as soft evidence or likelihood evidence) in \grbn. This is done via the function
\code{setEvidence()}.

The clique potential representation in a Bayesian network gives
\begin{displaymath}
  p(x) \propto \psi(x) = \prod_{C} \psi_C(x_C).
\end{displaymath}

Here we recall that the whole idea in computations with Bayesian
networks is to avoid calculation the product on the right hand
side above. Instead computations are based on propagation (multiplying,
dividing and summing clique potentials $\psi_C$ in an appropriate
order, and such an appropriate order comes from a junction tree).
The normalizing constant, say $c=\sum_x \psi(x)$, comes out of
propagation as a ``by product''.

Suppose a set of nodes $E$ are known to have a specific value,
i.e. $x_E=x^*_E$. This is called hard evidence. The probability of
the event $x_E=x^*_E$ is
\begin{displaymath}
  p(x_E=x^*_E)=E_p\{I(x_E=x^*_E)\} = \sum_x I(x_E=x^*_E) p(x)
  = \frac{1}{c} \sum_x I(x_E=x^*_E) \psi(x)
\end{displaymath}

The computations are based on modifying the clique potentials $\psi_C$
by giving value zero to states in $\psi_C$ which are not consistent
with $x_E=x^*_E$. This can be achieved with an indicator function, say
$L_C(x_C)$ such that we obtain a set of new potentials $\tilde \psi_C(x)
= L_C(x_C) \psi_C(x_C)$. Propagation with these new potentials gives,
as a by product, $\tilde c=\sum \tilde \psi(x)$ where
$\tilde\psi(x)= \prod_C \tilde\psi_C(x_C)$. Consequently, we have
$p(x_E=x^*_E)=\tilde c / c$.

In a more general setting we may have non--negative weights $L(x)$ for
each value of $x$. We may calculate
\begin{displaymath}
  E_p\{L(X)\} = \sum_x L(x)p(x)
\end{displaymath}
If $L(X)$ factorizes as $L(X)= \prod_C L_C(X_C)$ then the computations are
carried out as outlined above, i.e.\ by the message passing scheme.


\subsection{An excerpt of the chest clinic network}
\label{sec:an-excerpt-chest}


Consider the following excerpt of the chest clinic network.

<<>>=
yn <- c("yes", "no")
a    <- cpt(~asia, values=c(1, 99), levels=yn)
t.a  <- cpt(~tub|asia, values=c(5, 95, 1, 99), levels=yn)

plist1 <- compileCPT(list(a, t.a))
chest1 <- grain(plist1)
querygrain(chest1, simplify = TRUE) 
@ %def

\subsection{Specifying hard evidence} 
\label{sec:hard-evidence}

Suppose we want to make a diagnosis about tuberculosis given the evidence that a person has recently been to Asia.
The function \code{setEvidence()} can  be used for this purpose. The following forms are equivalent

<<>>=
setEvidence(chest1, evidence=list(asia="yes"))
@ %def

We call such evidence hard evidence because the state of the variables are known with certainty.


% <<>>=
% querygrain(chest1, nodes="tub", simplify = TRUE)
% querygrain(chest1, evidence=list(asia="yes"), simplify = TRUE)
% @ %def

\subsection{Virtual evidence (also called soft or likelihood evidence)}
\label{sec:virt-evid-likel}

Suppose we do not know with certainty whether a patient has
recently been to Asia (perhaps the patient is too ill to
tell). However the patient (if he/she is Caucasian) may be unusually
tanned and this lends support to the hypothesis of a recent visit to
Asia.

To accommodate this setting we create an extended network with an extra
node for which we enter evidence.


However, it is NOT necessary to do
so in practice, because we may equivalently enter the virtual evidence
in the original network.

We can then introduce a new variable
\code{guess\_asia} with \code{asia} as its only parent.

<<>>=
g.a <- cpt(~guess_asia+asia, levels=yn,
              values=c(.8, .2, .1, .9))
@ %def

This reflects the assumption that for patients who have recently been
to Asia we would (correctly) guess so in 80\% of the cases, whereas for patients who have
not recently been to A we would (erroneously) guess that they have
recently been to Asia in 10\% of the cases.

<<>>=
plist2 <- compileCPT(list(a, t.a, g.a ))
chest2 <- grain(plist2)
querygrain(chest2, simplify = TRUE)
@ %def


Now specify the guess or judgment, that the person has recently been
to Asia:

<<>>=
querygrain(chest2, evidence=list(guess_asia="yes"),
           simplify=TRUE, exclude = FALSE)
@ %def

\subsection{Specifying virtual evidence}
\label{sec:spec-virt-evid}

The same guess or judgment can be specified as virtual evidence
(also called likelihood evidence) for the original network:

<<>>=
chest1_ve <- chest1 |> setEvidence(evidence=list(asia=c(.8, .1)))
chest1_ve |> querygrain(simplify = TRUE)
getEvidence(chest1_ve, short=FALSE)
@ %def

This also means that specifying that specifying \code{asia='yes'} can
be done as
<<>>=
querygrain(chest1, evidence=list(asia=c(1, 0)), simplify=T)
@ %def

\subsection{Extending networks to include other types of variables}
\label{sec:ixture}

\grbn\ only handles discrete variables with a finite state space, but
using likelihood evidence it is possible to work with networks with
both discrete and continuous variables (or other types of variables).
This is possible only when he networks have a specific structure. This
is possible when no discrete variable has non--discrete parents.

Take a simple example: Form a \bn\ for variables $x=(x_1, x_2)$. 
Conceptually augment this network with additional variables $y=(y_1, y_2)$ where
$y_1|x_1=k \sim N(\mu_k, v)$ and
$y_2|x_2=k \sim Poi(\lambda_k)$ for $k=1,2$. Also we make the assumption
that $y_1$ and $y_2$ are independent given $x=(x_1, x_2)$. This gives the DAG
below:

<<>>=
plot(dag(~y1:x1 + x2:x1 + y2:x2))
@ 


A \bn\ for $x$ can be constructed as:

<<>>=
u <- list(x1=yn, x2=yn)
x1 <- cpt(~x1, values=c(1, 3), levels=yn)
x2 <- cpt(~x2|x1, values=c(1, 3, 3, 1), levels=yn)
bn <- grain(compileCPT(x1, x2))
querygrain(bn, simplify=TRUE)
@

The augmentation of $y|x$ can go along these lines: The parameters describing $y|x$ are set to be:
<<>>=
v <- 2
mu <- c(mu1=2, mu2=5)
lambda <- c(lambda1=0, lambda2=7)
@ 

Suppose we observe $y_1 = y_1^*$. Then
\begin{displaymath}
  p(x|y_1= y_1^*)\propto p(x_1)p(x_2|x_1) p(y_1=y_1^*|x_1) =  p(x_1)p(x_2|x_1) L_1(x_1)
\end{displaymath}
where $L_1(x_1)$ denotes the likelihood. In a \bn\ setting this
corresponds to changing $p(x_1)$ as
\begin{displaymath}
  p(x_1) \leftarrow p(x_1)L_1(x_1)
\end{displaymath}
and then carry on with propagation. This can be achieved in different ways. One is by setting the likelihood as evidence:

<<>>=
y1 <- 1 # Observed value for y1
lik1 <- dnorm(y1, mean=mu, sd=sqrt(v))
querygrain(bn, exclude = FALSE,
           evidence=list(x1=lik1), simplify = TRUE) 
@ 

An alternative is to explicitly modify the CPT which specifies $p(x_1)$:
<<>>=
x1_upd <- getgrain(bn, "cptlist")$x1 * lik1
bn2 <- replaceCPT(bn, list(x1=x1_upd))
querygrain(bn2) 
@ 

A final remark: The conditional distribution of
$y_1$ is normal, but the unconditional distribution is a mixture of
normals. Likewise, the conditional distribution of
$y_2$ is poisson, but the unconditional distribution is a mixture of
poissons. Evidence on, say $y_1$ changes the belief in $x_1$ and
$x_2$ and this in turn changes the distribution of
$y_2$ (evidence changes the mixture weights.)  


<<>>=
set.seed(2022)
nsim <- 1000
xsim1 <- simulate(bn, nsim)
head(xsim1)
xsim2 <- simulate(bn2, nsim)
head(xsim2)

par(mfrow=c(1,2))
y2sim <- rpois(n=nsim, lambda=lambda[xsim1$x2])
y22sim <- rpois(n=nsim, lambda=lambda[xsim2$x2])
y2sim |> hist(prob=T, ylim=c(0, .4), breaks=10)
y22sim |> hist(prob=T, ylim=c(0, .4), breaks=10)
@ 








The joint distribution is
\begin{displaymath}
  p(x,y_1, y_2) = p(x)p(y_1|x)p(y_2|x)
\end{displaymath}

Suppose the interest is in the distribution of $x$ given
$y_1=y_1^*$ and $y_2=y_2^*$. We then have
\begin{displaymath}
  p(x|y_1^*, y_2^*) \propto p(x) p(y_1^*|x)p(y_2^*|x) =
  p(x) L_1(x) L_2(x)
\end{displaymath}



\appendix

\section{Building networks from data}
\label{sec:using-textttsm-argum}

The following two graphs specify the same model:
<<>>=
dG  <- dag(~A:B + B:C)
uG  <- ug(~A:B + B:C)
par(mfrow=c(1,2)); plot( dG ); plot( uG )
@ %def

Suppose data is
<<>>=
dat <- tabNew(c("A", "B", "C"), levels=c("lev1", "lev2"), #levels=c(2,2,2),
              values=c(0, 0, 2, 3, 1, 2, 1, 4))
class(dat)
@ %def

A network can be built from data using:

<<>>=
gr.dG <- compile(grain( dG, data=dat ))
gr.uG <- compile(grain( uG, data=dat ))
@ %def

However, when there are zeros in the table, care must be taken.

\subsection{Extracting information from tables}
\label{sec:extr-inform-from}

In the process of creating networks, conditional probability tables
are extracted when the graph is a dag and clique potentials are
extracted when the graph is a chordal (i.e.\ triangulated) undirected
graph. This takes place as follows (internally):

<<>>=
extractCPT(dat, dG) |> c() ## FIXME: Printing problem
extractPOT(dat, uG) |> c() ## FIXME: Printing problem
@ %def

The conditional probability table $P(A|B)$ contains \code{NaN}s
because
\begin{displaymath}
  P(A|B=B1)=\frac{n(A,B=B1)}{\sum_A n(A,B=B1)} = \frac{0}{0} = \mbox{NaN}
\end{displaymath}

For this reason the network \code{gr.dG} above will fail to compile
whereas \code{gr.uG} will work, but it may not give the expected results.

\subsection{Using smooth}
\label{sec:using-smooth}

To illustrate what goes on, we can extract the distributions from data
as follows:

<<>>=
p.A_B <- tabDiv(dat, tabMarg(dat, "B")) ## p(A|B)
p.B   <- tabMarg(dat, "B") / sum(dat)   ## p(B)
p.AB2 <- tabMult(p.A_B, p.B)            ## P(AB)
@ %def

However, the result is slightly misleading because \code{tabDiv}
sets $0/0=0$.
In \grain\ there is a \code{smooth} argument that will add a small
number to the cell entries before extracting tables, i.e.
\begin{displaymath}
  P(A|B=B1)=\frac{n(A,B=B1)+\epsilon}{\sum_A ( n(A,B=B1) + \epsilon) }
  = \frac{\epsilon}{2\epsilon} = 0.5
\end{displaymath}
and
\begin{displaymath}
  P(B)= \frac{\sum_A (n(A,B)+\epsilon)}{\sum_{AB} (n(A,B)+\epsilon)}
\end{displaymath}

We can mimic this as follows:
<<>>=
e <- 1e-2
(dat.e <- dat + e)
@ %def

<<>>=
pe.A_B <- tabDiv(dat.e, tabMarg(dat.e, "B"))
pe.B   <- tabMarg(dat.e, "B") / sum(dat.e)
pe.AB  <- tabMult(pe.A_B, pe.B)

@ %def

However this resulting joint distribution is different from what is
obtained from the adjusted table itself
<<>>=
dat2.e <- dat.e / sum(dat.e)
@ %def

<<>>=
(dat2.e - pe.AB)  |> ftable()
@ 
This difference appears in the \grbn\ networks.

\subsection{Extracting tables}
\label{sec:extracting-tables}

One can do
<<>>=

gr.dG <- grain(dG, data=dat, smooth=e)
@ %def

which (internally) corresponds to
<<>>=
extractCPT(dat, dG, smooth=e)  |> c()
@ %def

We get
<<>>=
querygrain(gr.dG, exclude=FALSE, simplify=TRUE) 
querygrain(gr.uG, exclude=FALSE, simplify=TRUE) 
@ %def

However, if we condition on \code{B=lev1} we get:
<<>>=
querygrain(gr.dG, evidence=list(B="lev1"), exclude=FALSE, simplify=TRUE)
querygrain(gr.uG, evidence=list(B="lev1"), exclude=FALSE, simplify=TRUE)
@ %def

so the ``problem'' with zero entries shows up in a different
place. However, the answer is not necessarily wrong; the answer simply
states that $P(A|B=lev1)$ is undefined.
To ``remedy'' we can use the \code{smooth} argument:
<<>>=
gr.uG <- grain(uG, data=dat, smooth=e)
@ %def
which (internally) corresponds to
<<>>=
extractPOT(dat, uG, smooth=e) |> c()
@ %def

Notice that the results are not exactly identical:

<<>>=
querygrain(gr.dG, exclude=FALSE, simplify=TRUE) 
querygrain(gr.uG, exclude=FALSE, simplify=TRUE) 
@ %def


<<>>=
querygrain(gr.dG, evidence=list(B="lev1"), exclude=FALSE, simplify=TRUE) 
querygrain(gr.uG, evidence=list(B="lev1"), exclude=FALSE, simplify=TRUE) 
@ %def




\section{Brute force computations and why they fail}
\label{sec:brute-force-comp}


The \grbn\ package makes computations as those outlined above in a
very efficient way; please see the references.  However, it is in this
small example also possible to make the computations directly: We can
construct the joint distribution (an array with $2^8=256$ entries) directly as:
<< >>= 
joint <- tabListMult(chest_cpt)
dim(joint)
joint  |> as.data.frame.table() |> head()
@

This will clearly fail even moderate size problems: For example, a
model with $80$
nodes each with $10$
levels will give a joint state space with $10^{80}$
states; that is about the number of atoms in the universe. Similarly,
$265$
binary variables will result in a joint state space of about the same
size. Yet, \grbn\ has been used successfully
in models with tens of
thousand variables.  The ``trick'' in \grbn\ is to make all
computations without ever forming the joint distribution. 

However, we
can do all the computations by brute force methods as we will
illustrate here:

Marginal distributions are
<< >>= 
tabMarg(joint, "lung")
tabMarg(joint, "bronc")
@

Conditioning on evidence can be done in different ways: The
conditional density is a $6$--way slice of the original $8$--way joint
distribution:

<< >>= 
ev <- list(asia="yes", dysp="yes")
cond1 <- tabSlice(joint, slice=ev)
cond1 <- cond1 / sum(cond1)
dim(cond1)
tabMarg(cond1, "lung")
tabMarg(cond1, "bronc")
@

Alternatively, multiply all entries not consistent by zero and all other entries by one and then marginalize:
<< >>= 
cond2 <- tabSliceMult(joint, slice=ev)
cond2 <- cond2 / sum(cond2)
dim(cond2)
tabMarg(cond2, "lung")
tabMarg(cond2, "bronc")
@


\bibliography{grain}

\end{document}



% << >>= 
% f1 <- function(){
%     cond1 <- ar_slice(joint, slice=ev)
%     cond1 <- cond1 / sum(cond1)
%     dim(cond1)
%     tabMarg(cond1, "lung")    
% }
% f2 <- function(){
%     cond2 <- ar_slice_mult(joint, slice=ev)
%     cond2 <- cond2 / sum(cond2)
%     dim(cond2)
%     tabMarg(cond2, "lung")    
% }
% f1()
% f2()
% querygrain(chest_bn, nodes="lung", evidence=ev)
% chest_bn1 <- propagate(chest_bn)
% if (require(microbenchmark)){
%     microbenchmark(f1(), f2(),
%                    f3=querygrain(chest_bn, nodes="lung", evidence=ev),
%                    f4=querygrain(chest_bn, nodes="lung", evidence=ev)                   
%                    )                   
% }
% @


%%\SweaveInput{Rmarkup.STY}
%% ------------------------
%% \definecolor{darkred}{rgb}{.7,0,0}
%% \definecolor{midnightblue}{rgb}{0.098,0.098,0.439}
%% 
%% \DefineVerbatimEnvironment{Sinput}{Verbatim}{
%%   fontfamily=tt,
%%   %%fontseries=b,
%%   %% xleftmargin=2em,
%%   formatcom={\color{midnightblue}}
%% }
%% \DefineVerbatimEnvironment{Soutput}{Verbatim}{
%%   fontfamily=tt,
%%   %%fontseries=b,
%%   %% xleftmargin=2em,
%%   formatcom={\color{darkred}}
%% }
%% \DefineVerbatimEnvironment{Scode}{Verbatim}{
%%   fontfamily=tt,
%%   %%fontseries=b,
%%   %% xleftmargin=2em,
%%   formatcom={\color{blue}}
%% }
%% 
%% \fvset{listparameters={\setlength{\topsep}{-2pt}}}
%% \renewenvironment{Schunk}{\linespread{.90}}{}
%% %% ------------------------
%% 




% We can look closer into this zero--probability issue. Because the node
% \code{either} is logical, half of the configurations will have zero probability:

% <<>>=
% tt <- querygrain(chest_bn, type="joint")
% sum(tt == 0) / length(tt)
% @ %def

% In particular the configuration above has zero probability
% <<>>=
% sum(ar_slice(tt, list(either="no", tub="yes")))
% @ %def


