# A new class of generative classifiers based on staged tree models

Federico Carli<EMAIL_ADDRESS>
Dipartimento di Matematica, Università degli Studi di Genova, Genova, Italy
Manuele Leonelli<EMAIL_ADDRESS>
School of Science and Technology, IE University, Madrid, Spain  Gherardo
Varando<EMAIL_ADDRESS>
Image Processing Laboratory, Universitat de València, València, Spain

###### Abstract

Generative models for classification use the joint probability distribution of
the class variable and the features to construct a decision rule. Among
generative models, Bayesian networks and naive Bayes classifiers are the most
commonly used and provide a clear graphical representation of the relationship
among all variables. However, these have the disadvantage of highly
restricting the type of relationships that could exist, by not allowing for
context-specific independences. Here we introduce a new class of generative
classifiers, called staged tree classifiers, which formally account for
context-specific independence. They are constructed by a partitioning of the
vertices of an event tree from which conditional independence can be formally
read. The naive staged tree classifier is also defined, which extends the
classic naive Bayes classifier whilst retaining the same complexity. An
extensive simulation study shows that the classification accuracy of staged
tree classifiers is competitive with that of state-of-the-art classifiers and
an example showcases their use in practice.

Keywords: Bayesian networks; Model selection; Staged trees; Statistical
classification

##  1 Introduction

The aim of statistical classification is to assign labels to instances
described by a vector of feature variables. The classification task is guided
by a statistical model learnt using data containing labelled instances. The
array of models designed to perform classification is constantly increasing
and includes, among others, random forests (Ho, 1995), recursive partitioning
(Breiman et al., 1984) and probabilistic neural networks (Specht, 1990).

Bayesian network classifiers (BNCs) (Bielza and Larrañaga, 2014; Friedman et
al., 1997) are special types of Bayesian networks (BNs) designed for
classification problems. These have been applied to a wide array of real-world
applications with competitive classification performance against state-of-the-
art classifiers (Flores et al., 2012). There are numerous advantages
associated to BNCs. First, they provide an explicit and intuitive
representation of the relationship among features represented by a graph.
Second, they are a fully coherent probabilistic model thus giving uncertainty
measures about the chosen labels. Third, many of the methods and algorithms
developed for general BNs can be simply adapted and used for BNCs (see e.g.
Benjumeda et al., 2019). Last, they are implemented in various pieces of user-
friendly software (see e.g. the bnclassify R package of Mihaljevic et al.,
2018). More generally, BNCs are generative classifiers which give an estimate
of the joint distribution of both the features and the class.

One of the main limitations of BNs is that they can only explicitly represent
symmetric conditional independences among variables of interest. However, in
many applied domains, conditional independences are context-specific, meaning
that they only hold for specific instantiations of the conditioning variables.
For this reason, numerous extensions of BNs have been proposed that can take
into account asymmetric independences (Boutilier et al., 1996; Cano et al.,
2012; Jaeger et al., 2006; Pensar et al., 2015, 2016; Poole and Zhang, 2003).
With the exception of Jaeger et al. (2006) and Pensar et al. (2015), all these
models somehow lose the intuitiveness of BNs since they cannot represent all
the models’ information into a unique graph.

Despite the efforts in accommodating asymmetries in BNs, the development of
BNCs embedding context-specific information has been limited. Solutions
proposed to address this issue use the idea of Bayesian multinets (Geiger and
Heckerman, 1996), which consist of several networks each associated with a
subset of the domain of one variable, often called distinguished. Bayesian
multinets have been used for classification in Friedman et al. (1997), Gurwicz
and Lerner (2006), Huang et al. (2003) and Hussein and Santos (2004), among
others.

In this paper a novel class of generative classifiers based on staged trees
(Collazo et al., 2018; Smith and Anderson, 2008) embedding asymmetric
conditional independences is considered, and henceforth called _staged tree
classifiers_. Staged trees are defined as probability trees (Shafer, 1996)
embellished with asymmetric conditional independence information, represented
by a partitioning of the vertices of the tree. Staged tree classifiers are
staged trees whose topology is specifically designed for classification
problems. Although staged trees have been used in a variety of applications,
including the modelling of health problems (Barclay et al., 2013; Keeble et
al., 2017) and criminal activities (Collazo and Smith, 2016), their specific
use for classification problems has been limited.

Staged tree classifiers are generative classifiers which, whilst extending the
class of BNCs to deal with asymmetric conditional independences, share the
same advantages of BNCs: first, the relationship between the random variables
is still intuitively depicted in a unique graph; second, they are fully
coherent probabilistic models; third, learning algorithms already defined for
staged trees (e.g. Freeman and Smith, 2011; Silander and Leong, 2013; Leonelli
and Varando, 2022a, b) can simply be adapted for classification purposes;
fourth, the freely-available R package stagedtrees (Carli et al., 2020) gives
an implementation of a variety of learning algorithms and inferential routines
to apply the methods in practice.

Mirroring the theory of BNCs, staged tree classifiers of different complexity
are discussed and their properties investigated. Notably, the naive staged
tree classifier is introduced which is shown to have the same complexity of
the standard naive Bayes classifier, but relaxing the strict conditional
independence assumptions of naive Bayes. Experimental studies demonstrate that
staged tree classifiers have comparable classification rates to state-of-the-
art classifiers and outperform other generative classifiers.

The paper is structured as follows. Section 2 introduces the notation and
reviews BNCs. Staged trees are reviewed in Section 3. Staged tree classifiers
are introduced and studied in Section 4. Section 5 discusses learning
algorithms. Section 6 presents an experimental study and Section 7 discusses a
classification application in details. The paper is concluded with a
discussion.

##  2 Bayesian network classifiers

###  2.1 The Bayesian network model

Let $G=([p],E_{G})$ be a directed acyclic graph (DAG) with vertex set
$[p]=\\{1,\dots,p\\}$ and edge set $E_{G}$. Let $\bm{X}=(X_{i})_{i\in[p]}$ be
categorical random variables with joint mass function $P$ and sample space
$\mathbb{X}=\times_{i\in[p]}\mathbb{X}_{i}$. For $A\subset[p]$, we let
$\bm{X}_{A}=(X_{i})_{i\in A}$ and $\bm{x}_{A}=(x_{i})_{i\in A}$ where
$\bm{x}_{A}\in\mathbb{X}_{A}=\times_{i\in A}\mathbb{X}_{i}$. We say that $P$
is Markov to $G$ if, for $\bm{x}\in\mathbb{X}$,

$P(\bm{x})=\prod_{k\in[p]}P(x_{k}|\bm{x}_{\Pi_{k}}),$ (1)

where $\Pi_{k}$ is the parent set of $k$ in $G$ and
$P(x_{k}|\bm{x}_{\Pi_{k}})$ is a shorthand for
$P(X_{k}=x_{k}|\bm{X}_{\Pi_{k}}=\bm{x}_{\Pi_{k}})$. Henceforth, we assume the
existence of a linear ordering $\sigma$ of $[p]$ for which only pairs $(i,j)$
where $i$ appears before $j$ in the order can be in the edge set.

The ordered Markov condition implies conditional independences of the form

$X_{i}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{X}_{[i-1]}\,|\,\bm{X}_{\Pi_{i}}.$ (2)

###### Definition 1

Let $G$ be a DAG and $P$ Markov to $G$. The _Bayesian network_ model
(associated to $G$) is

$\mathcal{M}_{G}=\\{P\in\Delta_{|\mathbb{X}|-1}\,|\,P\mbox{ is Markov to
}G\\}.$

where $\Delta_{|\mathbb{X}|-1}$ is the ($|\mathbb{X}|-1$)-dimensional
probability simplex.

12345 Figure 1: A simple DAG with vertex set $\\{1,2,3,4,5\\}$ and edge set
$\\{(1,3),(1,4),(2,4),(3,5),(4,5)\\}$.

Figure 1 reports a simple DAG with five vertices. Any Markov distribution to
this DAG must factorize as
$P(x_{5}|x_{3},x_{4})P(x_{4}|x_{1},x_{2})P(x_{3}|x_{1},x_{2})P(x_{2})P(x_{1}).$

Let $\mathcal{G}$ be the set of DAGs with vertex set $[p]$ and ordering
$\sigma$. We define the space of BN models over $\bm{X}$ as
$\mathcal{M}_{\mathcal{G}}=\cup_{G\in\mathcal{G}}\mathcal{M}_{G}$

###  2.2 Bayesian networks for classification

Supppose now that $\bm{X}=(X_{1},\dots,X_{p})$ is a $p$-dimensional vector of
categorical feature variables and $C$ a categorical class variable with sample
space $\mathbb{C}$ and $c\in\mathbb{C}$. Given a training set of labelled
observations $\mathcal{D}=\\{(\bm{x}^{1},c^{1}),\dots,$
$(\bm{x}^{N},c^{N})\\}$, where $\bm{x}^{i}\in\mathbb{X}$ and
$c^{i}\in\mathbb{C}$, the aim of a generative classifier is to learn a joint
probability $p(c,\bm{x})$ and assign a non-labelled instance $\bm{x}$ to the
most probable a posteriori class found as

$\arg\max_{c\in\mathbb{C}}p(c|\bm{x})=\arg\max_{c\in\mathbb{C}}p(c,\bm{x}).$

Such a classifier is referred to as Bayes classifier.

BNCs are Bayes classifiers that factorize $p(c,\bm{x})$ according to a BN over
the variables $X_{1},\dots,X_{p}$ and $C$. Although any BN model could be used
for classification purposes, most often the underlying DAG is restricted so
that the class variable $C$ has no parents. Therefore, in BNCs the class
variable is the root of the DAG and there is a direct link from $C$ to
$X_{i}$, for $i\in[p]$, since otherwise features not connected to the class
would not provide any information for classification. The simplest possible
model is the so-called naive Bayes classifier (Minsky, 1961) which assumes the
features are conditionally independent given the class (Figure 2a). BNCs of
increasing complexity can then be defined by adding dependences between the
feature variables. For instance, the super-parent-one-dependence-estimator
(SPODE) BNC (Keogh and Pazzani, 2002) assumes there is a feature parent of all
others (Figure 2b). Another commonly used classifier is the tree-augmented
naive (TAN) BNC (Friedman et al., 1997) for which each feature has at most two
parents: the class and possibly another feature (Figure 2c).

$C$$X_{1}$$X_{2}$$X_{3}$

(a) Naive BNC

$C$$X_{2}$$X_{1}$$X_{3}$

(b) SPODE BNC

$C$$X_{1}$$X_{2}$$X_{3}$

(c) TAN BNC

Figure 2: Examples of BNCs with three features and one class.

Although BNCs of any complexity can be learnt and used in practice, empirical
evidence demonstrates that model complexity does not necessarily implies
better classification accuracy (Bielza and Larrañaga, 2014). Despite of their
simplicity, naive BNCs have been shown to lead to good accuracy in
classification problems (Bielza and Larrañaga, 2014; Flores et al., 2012).

Alongside these empirical evaluations of BNCs, theoretical studies about the
expressiveness of such models have appeared. Recently, Varando et al. (2015)
and Varando et al. (2016) fully characterized the decision functions induced
by various BNCs and consequently derived bounds for their expressive power.
They built on the work of Ling and Zhang (2002) that demonstrated that any BNC
whose vertices have at most $k$ parents cannot represent any decision function
containing $(k+1)$-XORs (also known as parity functions O’Donnell, 2014). Thus
naive BNCs are not capable of capturing any 2-XORs. On the positive side,
Domingos and Pazzani (1997) demonstrated that naive BNCs are optimal under a
$0-1$ loss even when the assumption of conditional independence among features
does not hold.

##  3 Staged trees

Because of the strict assumptions on the symmetry of independences among
variables in BNs, models accommodating asymmetric conditional independence
statements have been developed. Staged trees are one extension of BNs whose
conditional independences can still directly be read from the associated
graphical representation (Collazo et al., 2018; Smith and Anderson, 2008).
However, differently to BNs, whose graphical representation is a DAG, staged
trees are constructed from probability trees as detailed next.

###  3.1 $X$-compatible staged trees

Consider a $p$-dimensional random vector $\bm{X}$ taking values in the product
sample space $\mathbb{X}$. Let $(V,E)$ be a directed, finite, rooted tree with
vertex set $V$, root node $v_{0}$ and edge set $E$. For each $v\in V$, let
$E(v)=\\{(v,w)\in E\\}$ be the set of edges emanating from $v$ and
$\mathcal{L}$ be a set of labels.

###### Definition 2

An $\bf X$-compatible staged tree is a triple $(V,E,\theta)$, where $(V,E)$ is
a rooted directed tree and:

  1. 1.

$V={v_{0}}\cup\bigcup_{i\in[p]}\mathbb{X}_{[i]}$;

  2. 2.

For all $v,w\in V$, $(v,w)\in E$ if and only if
$w=\bm{x}_{[i]}\in\mathbb{X}_{[i]}$ and $v=\bm{x}_{[i-1]}$, or $v=v_{0}$ and
$w=x_{1}$ for some $x_{1}\in\mathbb{X}_{1}$;

  3. 3.

$\theta:E\rightarrow\mathcal{L}^{*}=\mathcal{L}\times\cup_{i\in[p]}\mathbb{X}_{i}$
is a labelling of the edges such that
$\theta(v,\bm{x}_{[i]})=(\kappa(v),x_{i})$ for some function
$\kappa:V\to\mathcal{L}$. The function $k$ is called the colouring of the
staged tree $T$.

If $\theta(E(v))=\theta(E(w))$ then $v$ and $w$ are said to be in the same
_stage_.

Therefore, the equivalence classes induced by $\theta(E(v))$ form a partition
of the internal vertices of the tree in _stages_.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$
Figure 3: An example of an $\bf X$-compatible staged tree.

Definition 2 first constructs a rooted tree where each root-to-leaf path, or
equivalently each leaf, is associated to an element of the sample space
$\mathbb{X}$. Then a labeling of the edges of such a tree is defined where
labels are pairs with one element from a set $\mathcal{L}$ and the other from
the sample space $\mathbb{X}_{i}$ of the corresponding variable $X_{i}$ in the
tree. By construction, $\bf X$-compatible staged trees are such that two
vertices can be in the same stage if and only if they correspond to the same
sample space. Although staged trees can be more generally defined without
imposing this condition (see e.g. Collazo et al., 2018), henceforth, and as
common in practice, we focus on $\bf{X}$-compatible staged trees only (see
Leonelli, 2019, for an example of a non $\bf{X}$-compatible tree).

Figure 3 reports an $(X_{1},X_{2},X_{3})$-compatible staged tree over three
binary variables. The coloring given by the function $\kappa$ is shown in the
vertices and each edge $(\cdot,(x_{1},\ldots,x_{i}))$ is labeled with
$X_{i}=x_{i}$. The edge labeling $\theta$ can be read from the graph combining
the text label and the color of the emanating vertex. The staging of the
staged tree in Figure 3 is given by the partition $\\{v_{0}\\}$,
$\\{v_{1}\\}$, $\\{v_{2}\\}$, $\\{v_{3},v_{4}\\}$ and $\\{v_{5},v_{6}\\}$.

The parameter space associated to an $\bf X$-compatible staged tree
$T=(V,E,\theta)$ with labeling $\theta:E\rightarrow\mathcal{L}^{*}$ is defined
as

$\Theta_{T}=\Big{\\{}\bm{y}\in\mathbb{R}^{|\theta(E)|}\;|\;\forall~{}e\in
E,y_{\theta(e)}\in(0,1)\textnormal{ and }\sum_{e\in
E(v)}y_{\theta(e)}=1\Big{\\}}$ (3)

Equation (3) defines a class of probability mass functions over the edges
emanating from any internal vertex coinciding with conditional distributions
$P(x_{i}|\bm{x}_{[i-1]})$, $\bm{x}\in\mathbb{X}$ and $i\in[p]$. In the staged
tree in Figure 3 the staging $\\{v_{3},v_{4}\\}$ implies that the conditional
distribution of $X_{3}$ given $X_{1}=0$, and $X_{2}=0$, represented by the
edges emanating from $v_{3}$, is equal to the conditional distribution of
$X_{3}$ given $X_{1}=0$ and $X_{2}=1$. A similar interpretation holds for the
staging $\\{v_{5},v_{6}\\}$. This in turn implies that
$X_{3}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{2}|X_{1}$, thus illustrating that the
staging of a tree is associated to conditional independence statements.

Let $\bm{l}_{T}$ denote the leaves of a staged tree $T$. Given a vertex $v\in
V$, there is a unique path in $T$ from the root $v_{0}$ to $v$, denoted as
$\lambda(v)$. The number of edges in $\lambda(v)$ is called the distance of
$v$, and the set of vertices at distance $k$ is denoted by $V_{k}$. For any
path $\lambda$ in $T$, let $E(\lambda)=\\{e\in E:e\in\lambda\\}$ denote the
set of edges in the path $\lambda$.

###### Definition 3

The _staged tree model_ $\mathcal{M}_{T,\theta}$ associated to the $\bf
X$-compatible staged tree $(V,E,\theta)$ is the image of the map

$\displaystyle\phi_{T}$
$\displaystyle:\Theta_{T}\rightarrow\Delta_{|\bm{l}_{T}|-1}^{\circ};$ (4)
$\displaystyle\phi_{T}$ $\displaystyle:y\mapsto p_{l}=\Big{(}\prod_{e\in
E(\lambda(l))}y_{\theta(e)}\Big{)}_{l\in\bm{l}_{T}}$

Therefore, staged trees models are such that atomic probabilities are equal to
the product of the edge labels in root-to-leaf paths and coincide with the
usual factorization of mass functions via recursive conditioning.

Let $\Theta$ be the set of functions $\theta$ from $E$ to $\mathcal{L}^{*}$,
that is all possible partitions, or staging, of the staged tree. We define
$\mathcal{M}_{T}=\cup_{\theta\in\Theta}\mathcal{M}_{T,\theta}$. So as
$\mathcal{M}_{\mathcal{G}}$ is the union of all possible BN models given a
specific ordering, $\mathcal{M}_{T}$ is the union of all possible staged tree
models, that is of all possible stagings, given a specific ordering of the
variables.

###  3.2 Staged trees and Bayesian networks

Smith and Anderson (2008) demonstrated that any BN can be represented as an
equivalent staged tree, whilst the converse is not true. We next illustrate
how to construct a staged tree equivalent to a BN. Consider a DAG $G$ and an
$\bf X$-compatible staged tree with vertex set $V$, edge set $E$ and labeling
$\theta$ defined via the coloring $\kappa(\bm{x}_{[i]})=\bm{x}_{\Pi_{i}}$ of
the vertices. The staged tree $T_{G}$, with vertex set $V$, edge set $E$ and
labeling $\theta$ so constructed, is called _the staged tree model of $G$_.
Importantly, $\mathcal{M}_{G}=\mathcal{M}_{T_{G},\theta}$, i.e. the two models
are exactly the same, since they entail exactly the same factorization of the
joint probability (Varando et al., 2021). The staging of $T_{G}$ represents
the Markov conditions associated to the graph $G$.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$X_{1}=0$$X_{1}=1$
$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$
Figure 4: An example of a staged tree not associated to any BN $G$.

For instance, the staged tree in Figure 3 can be constructed as the $T_{G}$
from the BN with DAG $X_{2}\leftarrow X_{1}\rightarrow X_{3}$. Conversely,
consider the staged tree in Figure 4. The blue staging implies that the
conditional distribution of $X_{3}$ given $X_{2}=X_{1}=0$ is equal to the
conditional distribution of $X_{3}$ given $X_{2}=X_{1}=1$. Such a constraint
cannot be explicitly represented by the DAG of a BN and therefore there is no
DAG $G$ such that $\mathcal{M}_{G}=\mathcal{M}_{T_{G},\theta}$, i.e. there is
no BN which is equivalent to the staged tree in Figure 4. More generally, it
holds that $\mathcal{M}_{\mathcal{G}}\subset\mathcal{M}_{T}$ (Varando et al.,
2021).

##  4 Staged tree classifiers

The technology of staged trees has been refined over the years and methods to
investigate causal relationships (Genewein et al., 2020; Thwaites et al.,
2010), perform statistical inference (Görgen et al., 2015), check model’s
robustness (Leonelli et al., 2017) and carry out causal discovery (Leonelli
and Varando, 2021) are now available. However, the specific use of staged
trees for classification has not been investigated in the literature. Thus,
just as BNCs have been defined as a specific subclass of BNs whose graph
entertains some properties, the class of staged tree classifiers is defined
here.

As in Section 2, suppose $\bm{X}=(X_{1},\dots,X_{n})$ is a vector of features
and $C$ is the class variable.

###### Definition 4

A staged tree classifier for the class $C$ and features $\bm{X}$ is a
$(C,\bm{X})$-compatible staged tree.

The requirement of $C$ being the root of the tree follows from the idea that
in most BNCs the class has no parents, so to maximise the information provided
by the features for classification.

###  4.1 The relationship between BNCs and staged tree classifiers

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$C=0$$C=1$$X_{1}=0$$X_{1}=1$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$

(a)

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$C=0$$C=1$$X_{1}=0$$X_{1}=1$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$

(b)

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$C=0$$C=1$$X_{1}=0$$X_{1}=1$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$

(c)

Figure 5: Representation of BNCs as staged trees classifiers: (a): naive BNC
in Figure 2a as a staged tree classifier; (b): SPODE BNC in Figure 2b as a
staged tree classifier; (c): TAN BNC in Figure 2c as a staged tree classifier.
All variables are assumed binary.

The BNCs reviewed in Section 2 can now be represented as staged tree
classifiers. Since a BNC is a BN with a DAG $G$, one can construct its
equivalent staged tree $T_{G}$ as in Section 3.2. For instance naive BNCs
(Figure 5a), SPODE BNCs (Figure 5b) and TAN BNCs (Figure 5c) can concisely be
represented as staged tree classifiers. However, the class of staged trees
classifiers is much larger than that of BNCs, as formalized in Proposition 5.
For a class $C$ and features $\bm{X}$, let
$\mathcal{M}_{\mathcal{G}}^{\textnormal{C}}$ be the space of BNCs (where $C$
is the root), and $\mathcal{M}_{T}^{\textnormal{C}}$ the space of
$(C,\bm{X})$-compatible staged trees.

###### Proposition 5

$\mathcal{M}_{\mathcal{G}}^{\textnormal{C}}\subset\mathcal{M}_{T}^{\textnormal{C}}$.

The proof follows from Varando et al. (2021).

In particular naive Bayes, SPODE, and TAN classifiers are all staged tree
classifiers with a specific staging structure as described next for the naive
Bayes. Recall that the set $V_{k}$ includes the nodes of the tree at distance
$k$ from the root and let $T(v)$ be the subtree of $T$ rooted at a vertex $v$.

###### Proposition 6

Let $G$ be the DAG of a naive BNC. Then $T_{G}$ is a $(C,\bm{X})$-compatible
staged tree where, for all $v\in V_{1}$, the subtree $T_{G}(v)$ is a
$\bm{X}$-compatible staged tree where all nodes at the same distance from the
root are in the same stage. In particular, $T_{G}$ has $|\mathcal{C}|$ stages
per each feature.

###  4.2 Conditional independence in staged tree classifiers

For the specific task of classification, it is possible to derive two results
about the dependence between the features and the class in staged tree
classifiers.

###### Proposition 7

If all $v\in V_{k}$ of a staged tree classifier are in the same stage then
$(C,X_{1},X_{2},\ldots,X_{k-1})$ and $X_{k}$ are marginally independent, i.e.
$(C,X_{1},\ldots,X_{k-1})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{k}$.

###### Proposition 8

If for all $v\in V_{1}$, $T(v)$ has the same stage structure over the vertices
at distance $k-2$ from the root, then $X_{k}$ is independent of $C$
conditionally on $X_{1},\dots,X_{k-1}$, i.e. $C\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{k}\,|\,X_{1},\dots,X_{k-1}$.

These two results are illustrated in Figure 6. For instance, consider the
features associated to the last random variable in Figure 6b. The vertices in
the upper half are framed as the vertices in the bottom half thus implying
that the class variable is conditionally independent of the last feature given
all others.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$C=0$$C=1$$X_{1}=0$$X_{1}=1$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$

(a)

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$C=0$$C=1$$X_{1}=0$$X_{1}=1$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$

(b)

Figure 6: Staged trees classifiers embedding conditional independence
statements between features and the class. (a):
$X_{3}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C$; (b): $X_{3}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C\,|\,X_{1},X_{2}$.

###  4.3 Naive staged tree classifiers

The class of staged tree classifiers is extremely rich and for any
classification task the number of candidate models that could explain the
relationship between class and features increases exponentially. One first
common assumption that we make here is to consider only
$(C,\bm{X})$-compatible staged trees, ones where only vertices at the same
distance from the root can be in a same stage. However, even with this
assumption the model class of staged tree classifiers is still much richer
than that of BNCs. Therefore, just like BNCs whose DAGs have restricted
topologies (as SPODE and TAN classifiers) have been studied, next we introduce
a class of simpler staged tree classifiers.

We have discussed that in many practical applications the naive BNC has very
good classification performance despite of its simplicity. A naive BNC has a
total of $\sum_{j=1}^{n}|\mathbb{C}|\cdot(|\mathbb{X}_{j}|-1)+|\mathbb{C}|-1$
free parameters that need to be learnt, whilst its DAG is always fixed.
Similarly, we introduce a class of staged tree classifiers which has the
constraint of having the same number of free parameters as naive BNCs, and
therefore has the same complexity, whilst being a much richer class of models
then the naive BNC.

###### Definition 9

A $(C,\bm{X})$-compatible staged tree classifier such that for every $k\leq
p$, the set $V_{k}$ is partitioned into $|\mathbb{C}|$ stages is called
_naive_.

It staightforwardly follows from the definition that naive staged tree
classifiers have the same number of free parameters as naive Bayes
classifiers.

Despite of the strict constraint on the number of parameters, the class of
naive staged tree classifiers is still rich and extends naive BNCs in a non-
trivial way. More formally, let
$\mathcal{M}_{\mathcal{G}}^{\textnormal{naive}}$ and
$\mathcal{M}_{T}^{\textnormal{naive}}$ be the space of naive BNCs and of naive
staged tree classifiers, respectively, for a class $C$ and features $\bm{X}$.

###### Proposition 10

$\mathcal{M}_{\mathcal{G}}^{\textnormal{naive}}\subset\mathcal{M}_{T}^{\textnormal{naive}}$.

Differently to naive BNCs, it is not sufficient to simply learn the
probabilities of the naive staged tree classifier, but also the staging
structure has to be discovered. However, because of the strict restriction on
the number of parameters, fast algorithms can be devised to efficiently
explore the model space. Notice that in a binary classification problem the
set $V_{k}$ must be partitioned into two subsets.

Critically and differently to naive Bayes classifiers, naive staged trees are
capable of representing complex decision rules. For instance, consider the
simplest scenario of a binary class with two binary features. The naive staged
tree classifier in Figure 4 (assuming the first variable is the class
variable), which, as already noticed, does not have a naive BNC
representation, is capturing the only 2-XOR present.

To investigate further the capabilities of naive staged tree classifiers in
expressing complex decision rules, we simulate $N_{train}=200$ observations
from $n=10$ binary variables $X_{1},\ldots,X_{10}$ taking values in
$\mathbb{X}=\\{-1,+1\\}^{n}$. We define the class variable as the parity (or
XOR) function $C=\prod_{i=1}^{n}X_{i}$ and compare naive Bayes, random forest
and naive staged tree classifiers over $N_{test}=10000$ test instances,
obtaining the results in Table 1. See Section 4.3 for details on the learning
of naive staged tree classifiers. As expected, the Naive Bayes classifier is
unable to represent the parity function (Varando et al., 2015) and wrongly
classifies the class in more than $40\%$ of the test data. Similar
performances are obtained by random forests (implemented with the randomForest
R package using $500$ trees), even if theoretically they have much larger
expressive power. Conversely, the naive staged tree correctly learns the
parity function with an accuracy of 90%.

| ST_Naive |  | RFor |  | NB
---|---|---|---|---|---
| -1 | 1 |  | -1 | 1 |  | -1 | 1
-1 | 0.4623 | 0.0713 |  | 0.3972 | 0.3123 |  | 0.4314 | 0.3982
1 | 0.0347 | 0.4317 |  | 0.0998 | 0.1907 |  | 0.0656 | 0.1048
Table 1: Proportion of predicted instances in the simulated XOR example for
the naive staged tree classifier (ST_Naive), random forest (RFor) and naive
Bayes (NB).

##  5 Learning staged tree classifiers

The learning of the structure of a staged tree from data is challenging due to
the exponential increase of the size of the tree with the number of random
variables. The first learning algorithm used for this purpose was the
agglomerative hierarchical clustering proposed by Freeman and Smith (2011).
Other learning algorithms were then introduced by Silander and Leong (2013)
and Collazo and Smith (2016), among others. Recently, a staged trees
implementation in R was made available in the stagedtrees package (Carli et
al., 2020) with various searching algorithms to estimate stage structures from
data.

In this article we present different methods for learning staged tree
classifiers. All methods follow three main steps: (i) an optimal order of the
variables is identified; (ii) the full probability tree based on all random
variables is pruned to speed-up computations; (iii) model search heuristics
are used to identify high-scoring structures. We next give details about each
of these phases. Furthermore, in Section 5.4 we introduce algorithms
specifically designed for naive staged tree classifiers.

###  5.1 Choosing a features’ order

It is well known that the ordering of the variables affect the quality of a
statistical classifier (e.g. Hruschka Jr and Ebecken, 2007). Thus, choosing
the order of the features in a staged tree classifier is expected to be
important in order to obtain good performances.

To confirm this, we perform an empirical study where we compute classification
accuracies for different staged tree classifiers over all possible features’
orders. Given the combinatorial explosion of the number of orders we limit
this study to two illustrative datasets with a small amount of features,
namely puffin and monks3 (see Table 2 for details). From Figure 7 we can
observe that, as expected, the order of the features is highly relevant with
respect to classification performance. For instance, the naive staged tree
learnt with hierarchical clustering, ST_Naive (as described in Section 5.4),
and the staged tree classifier learnt with the fast backward-hill climbing
algorithm, ST_FBHC (maximizing BIC score, see Section 5.3), exhibit a variety
of accuracies (from $0.5$ to $1$) for the same dataset.

Figure 7: Distribution of accuracy for different orders of features.
Accuracies obtained with the CMI ordering are shown with a cross. Results
using two learning algorithms for datasets puffin and monks3.

It is therefore critical to couple any algorithm to learn a staged tree with
an appropriate method to select a good ordering of the variables. For staged
tree classifiers, we tested various ordering heuristics and we propose here
the use of the conditional mutual information (CMI) criterion. The order of
the features given by the CMI criterion is obtained by iteratively selecting
the feature that maximise the conditional mutual information with respect to
the class variable, given the features previously selected. The performance of
such ordering is shown with a cross symbol in Figure 7. We can observe that in
both puffin and monks3 problems, the CMI order leads to a staged tree which
perfrom better than the majority of the possible orders.

###  5.2 Dealing with unobserved instances

Once the order of the variables is selected, it is possible that some of the
nodes in the tree are not observed in the training data. For such nodes,
estimating the associated probabilities is not possible without further
assumptions. As implemented in the R package stagedtrees, for each feature,
unobserved situations are joined in a common stage, called unobserved stage,
and a uniform probability is imposed for their distribution. Furthermore, such
unobserved stages are excluded from the stage structure search. This has the
effect of drastically reducing the number of nodes for which a stage structure
must be learnt, thus critically reducing training time.

From a statistical point of view, this does not affect the number of free
parameters of the corresponding statistical model, since probability
distributions of these unobserved stages have not to be estimated. This
procedure can be seen as the usual pruning step of learning algorithms for
classification trees. From the naive staged tree classifier perspective, it
follows that the number of stages in the tree is unaffected, since its
complexity, i.e. number of free parameters, does not change.

###  5.3 Learning algorithms

Different heuristic algorithms can be used to search the space of possible
stage structures. We rely on the available implementations in the stagedtrees
package, in particular we use greedy approaches that maximise the BIC score
(Görgen et al., 2020) or iterative methods that join nodes into the same
stages based on the distance between their associated probabilities. More
details of the algorithms are described in Carli et al. (2020) and in the
stagedtrees documentation.

###  5.4 Learning naive staged tree classifiers

As described in Section 4.3, a naive staged tree classifier is a staged tree
where the class is the first variable in the tree and nodes at the same
distance from the root are assigned to a fixed number of stages equal to the
number of possible values of the class (e.g. two for a binary classification
problem).

The two algorithms based on clustering (hierarchical and k-means) available in
the stagedtrees package can be used to learn naive staged tree classifiers.
They cluster probabilities into a user-selected number of stages which, for
the purposes of classification, can be fixed to number of possible values of
the class variable.

##  6 Experimental study

Dataset | # observations | # variables | # atomic events
---|---|---|---
asym | 100 | 4 | 16
breastCancer | 683 | 10 | 1024
chestSim500 | 500 | 8 | 256
energy1 | 768 | 9 | 1728
energy2 | 768 | 9 | 1728
fallEld | 5000 | 4 | 64
fertility | 100 | 10 | 15552
monks1 | 432 | 7 | 864
monks2 | 432 | 7 | 864
monks3 | 432 | 7 | 864
puffin | 69 | 6 | 768
ticTacToe | 958 | 10 | 39366
titanic | 2201 | 4 | 32
voting | 435 | 17 | 131072
Table 2: Details about the 14 datasets included in the experimental study.

The classification accuracy for binary classification of staged tree
classifiers is investigated in a comprehensive simulation study involving 14
datasets, whose details are given in Table 2. Each dataset is randomly divided
ten times in train set ($80\%$ of the data) to learn the classifiers and test
set (remaining $20\%$) to predict the response. The reported performance
measures, area under the curve (AUC) and balanced accuracy, are computed as
the mean over the ten replications.

First, 9 model search algorithms to learn staged tree classifiers are
compared, namely: ST_BHC (backward hill-climbing); ST_BJ_01 (backward joining
of vertices that have Kullback-Leibler divergence between their floret
probability distributions less than 0.01); ST_BJ_20 (as ST_BJ_01 but with
threshold at 0.20); ST_FBHC (a fast backward hill-climbing where two vertices
are joined whenever the score is increased); ST_Full (each vertex is in its
own stage); ST_HC_Full (hill-climbing algorithm starting from ST_Full);
ST_HC_Indep (hill-climbing algorithm starting from a tree where all vertices
associated to the same variable are in the same stage); ST_Naive_HC (naive
staged tree learnt with hierarchical clustering); ST_Naive_KM (naive staged
tree learnt with k-means). Further details about these algorithms can be found
in Carli et al. (2020). Due to computational restrictions, for ST_HC_Full the
model search is restricted to the first five features according to the
variable ordering chosen through CMI, whilst for ST_BHC and ST_HC_Indep only
the first seven are considered. The vertices corresponding to the remaining
variables are still used for classification but left as in the starting tree
of the model search.

Figure 8: AUC, balanced accuracy and logarithm of time spent for structure
learning for nine staged tree classifiers algorithms over fourteen datasets.

The results of the experiment are reported in Figure 8, which suggests the
following conclusions:

  * •

The ST_Full model (in yellow), which does not require any model search and has
the largest number of parameters, has in general lower AUC and balanced
accuracy than other staged trees. This highlights the need of a model-based
search of simpler models;

  * •

Models based on hill-climbing (in blue) overall perform better than others (in
particular ST_HC_Full). This is expected since these are the most refined
learning algorithms and, as a consequence, they are also the slowest.

  * •

Models based on backward joining (in green) have a satisfactory performance,
often comparable to that of hill-climbing models, whilst being much quicker to
learn.

  * •

Naive staged trees (in red) can be learnt extremely quickly and whilst often
they have a lower performance, there are cases where they are comparable to
the one of much more complex trees (see e.g. the balanced accuracy for the
titanic dataset)

Figure 9: AUC, balanced accuracy and logarithm of time spent for structure
learning for three staged tree classifiers (in red) and three BNCs (in blue)
over fourteen datasets.

Next, we compare staged trees classifiers with their competitor generative
classifier, namely BNCs. For ease of exposition, three representative staged
trees are selected (ST_BJ_01, ST_HC_Full and ST_Naive_KM) and three BNCs are
fitted: (i) the TAN BNC (BNC_TAN); (ii) the 3-dependence BNC (BNC_KDB); (iii)
the naive Bayes (BNC_NB). The results are reported in Figure 9. We can see
that for most datasets there is one staged tree classifier (in red) that
outperforms BNCs (in blue). Due to the complexity of the models, staged trees
are in general slower to learn, but the ST_Naive_KM, due to its simplicity,
has learning times comparable to those of generic BNCs.

Figure 10: AUC, balanced accuracy and logarithm of time spent for structure
learning for three staged tree classifiers (in red), two BNCs (in blue) and
other discriminative models (in green) over fourteen datasets.

Figure 10 reports the results of the simulation experiments for three staged
tree classifiers (ST_BJ_01, ST_HC_Full and ST_Naive_KM) as well as other
state-of-the-art generative and discriminative classifiers, namely: (i) the
naive Bayes (BNC_NB); (ii) the TAN BNC (BNC_TAN); (iii) Classification trees
(CTree) (iv) Logistic regression (Logistic); (v) Neural Networks with 20
hidden layers and 0.01 weight decay (NNet); (vi) Random Forests combining 100
classification trees (RFor). Although in some cases discriminative classifiers
(in green) outperform staged trees (in red), in many others they have
comparable AUC and balanced accuracy. However, as shown in the next section,
staged tree classifiers have the capability of producing an understanding of
the relationship between the class and the features, since they are generative
models. As already noticed, staged trees have an advantage over BNCs (in
blue). Although the learning time for generic staged trees is larger, the
learning time for naive staged tree classifiers is comparable to that of
state-of-the-art classifiers.

Figure 11: AUC, balanced accuracy and logarithm of time spent for structure
learning for two naive staged tree classifiers (in red) and two
implementations of naive Bayes classifiers (in blue) over fourteen datasets.

Last, we compare the performance of naive Bayes classifiers with the one of
naive staged tree classifiers in Figure 11. The overall conclusion is that in
most cases naive staged trees (in red) outperform naive Bayes in terms of AUC
and balanced accuracy. Furthermore, although naive staged trees require more
learning time since their structure has to be discovered, these can be learnt
very quickly and most times in less than one second.

##  7 An example of a staged tree classifier

To illustrate the capabilities of staged trees classifiers we next develop an
example classification analysis over the freely available titanic dataset,
which provides information on the fate of the Titanic passengers. It has three
binary variables (Survived, Sex and Age) and a categorical variable Class
taking four levels (1st, 2nd, 3rd and Crew). The aim is to correctly classify
whether the Titanic passengers survived or not based on their gender, age and
travelling class.

From Figure 8 we can see that one of the best staged tree classifiers is the
ST_BJ_01 learnt using a backward joining of the vertices based on the
Kullback-Leibler divergence and a threshold of 0.1. In Figure 12 we report the
staged tree classifier ST_BJ_01 learnt over the full Titanic dataset using the
R package stagedtrees. By investigating the staging structure we can deduce
conditional independence statements relating to the classification variable
(Survived) and the features. From stages associated to Class we can deduce
that
$P(\textnormal{Class}\,|\,\textnormal{Sex}=\textnormal{Male},\textnormal{Survived})=P(\textnormal{Class}\,|\,\textnormal{Sex}=\textnormal{Male})$
since the second and the fourth vertices (starting from the top) are in the
same stage. This implies the asymmetric conditional independence

$\textnormal{Class}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\textnormal{Survived}\,|\,\textnormal{Sex}=\textnormal{Male}$

The complex staging structure over the Age variable also implies asymmetric
conditional independences. We can notice that all paths going through an edge
labelled Crew are in the same stage for the variable Age. This implies that

$\textnormal{Age}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\textnormal{Survived}\,|\,\textnormal{Class}=\textnormal{Crew}$

The same conclusion can also be drawn for Class $=$ 3rd.

Figure 12: Staged tree classifier ST_BJ_01 learnt over the full Titanic
dataset.

As an additional illustration in Figure 13 is reported the naive staged tree
classifier learnt over the full Titanic dataset using the k-means hierarchical
clustering algorithm. The staging structure over the variables Sex and Class
implies that Sex and Survived are not independent and that Class is
conditionally independent of Survived given Sex. The staging structure over
the Age variable is a lot more complex describing highly asymmetric
constraints on the associated probabilities. Whilst imposing much more
flexible dependence structures, the naive staged tree classifier has the same
complexity of the naive Bayes classifier, meaning they have the same number of
independent parameters.

Figure 13: Naive staged tree classifier learnt over the full Titanic dataset
using the stages_kmeans algorithm.

BNCs of different complexity are also learnt over the full titanic dataset
using the bnclassify R package. Irrespective of the complexity chosen, the
model selection search always returns the simple naive Bayes classifier. Given
that staged tree classifiers outperfom BNCs in classification measures for the
Titanic dataset (see Figure 9), as well as for other datasets, this
observations highlights the need of context-specific generative classifiers
that can more flexibly model the dependence structure between the
classification variable and the features.

##  8 Discussion

Staged trees classifiers are a highly-expressive new class of generative
classifiers with classification performance comparable to that of state-of-
the-art classifiers. They embed context-specific conditional independence
statements which can be easily read by the staging of the vertices of the
tree. These are implemented in the freely available stagedtrees R package.

A special staged tree classifier is the naive staged tree classifier which,
whilst having the same complexity as the naive Bayes classifier, can flexibly
represent complex classification rule. Naive staged trees not only relax the
assumption of conditional indepedence of the features as in naive Bayes
classifiers but also have better performances in classification, as
highlighted by the simulation study.

Naive staged trees are learnt from data using a clustering algorithm of the
probability distributions over the non-leaf vertices of the tree. Such
algorithms divide the vertices at the same distance from the root in
$|\mathbb{C}|$ stages. More generally, we could devise clustering algorithms
were, for each variable, the number of stages is automatically selected
according to an optimality criterion. The development of these algorithms is
the focus of ongoing research.

Furthermore, many model search algorithms implemented in stagedtrees are based
on the maximization of a model score, by the default the negative BIC. We are
currently investigating algorithms based on the minimization of the
classification error, as commonly implemented for BNCs.

## Acknowledgments

Gherardo Varando’s work was partly funded by the European Research Council
(ERC) Synergy Grant “Understanding and Modelling the Earth System with Machine
Learning (USMILE)” under Grant Agreement No 855187.

## References

  * Barclay et al. (2013) L. Barclay, J. Hutton, and J. Smith.  Refining a Bayesian network using a chain event graph.  _International Journal of Approximate Reasoning_ , 54:1300–1309, 2013.
  * Benjumeda et al. (2019) M. Benjumeda, S. Luengo-Sanchez, P. Larranaga, and C. Bielza.  Tractable learning of Bayesian networks from partially observed data.  _Pattern Recognition_ , 91:190–199, 2019.
  * Bielza and Larrañaga (2014) A. Bielza and P. Larrañaga.  Discrete Bayesian network classifiers: a survey.  _ACM Computing Surveys_ , 47(1):5:1–5:43, 2014\.
  * Boutilier et al. (1996) C. Boutilier, N. Friedman, M. Goldszmidt, and D. Koller.  Context-specific independence in Bayesian networks.  In _Proceedings of the 12th Conference on Uncertainty in Artificial Intelligence_ , pages 115–123, 1996.
  * Breiman et al. (1984) L. Breiman, J. Friedman, R. Olshen, and C. Stone.  _Classification and regression trees_.  Wadsworth, 1984.
  * Cano et al. (2012) A. Cano, M. Gómez-Olmedo, S. Moral, C. Pérez-Ariza, and A. Salmerón.  Learning recursive probability trees from probabilistic potentials.  _International Journal of Approximate Reasoning_ , 53(9):1367–1387, 2012.
  * Carli et al. (2020) F. Carli, M. Leonelli, E. Riccomagno, and G. Varando.  The R package stagedtrees for structural learning of stratified staged trees.  _Journal of Statistical Software_ , 102(6):1–30, 2020.
  * Collazo and Smith (2016) R. Collazo and J. Smith.  A new family of non-local priors for chain event graph model selection.  _Bayesian Analysis_ , 11(4):1165–1201, 2016.
  * Collazo et al. (2018) R. Collazo, C. Görgen, and J. Smith.  _Chain event graphs_.  Chapmann & Hall, 2018.
  * Domingos and Pazzani (1997) P. Domingos and M. Pazzani.  On the optimality of the simple Bayesian network classifiers.  _Machine Learning_ , 29(2-3):103–130, 1997.
  * Flores et al. (2012) M. Flores, J. Gámez, and A. Martínez.  Supervised classification with Bayesian networks: a review on models and applications.  In _Intelligent data analysis for real-life applications: theory and practice_ , pages 72–102. 2012.
  * Freeman and Smith (2011) G. Freeman and J. Smith.  Bayesian MAP model selection of chain event graphs.  _Journal of Multivariate Analysis_ , 102(7):1152–1165, 2011.
  * Friedman et al. (1997) N. Friedman, D. Geiger, and M. Goldszmidt.  Bayesian network classifiers.  _Machine Learning_ , 29(2-3):131–163, 1997.
  * Geiger and Heckerman (1996) D. Geiger and D. Heckerman.  Knowledge representation and inference in similarity networks and Bayesian multinets.  _Artificial Intelligence_ , 82:45–74, 1996.
  * Genewein et al. (2020) T. Genewein, T. McGrath, G. Delétang, V. Mikulik, M. Martic, S. Legg, and P. Ortega.  Algorithms for causal reasoning in probability trees.  _arXiv:2010.12237_ , 2020.
  * Görgen et al. (2015) C. Görgen, M. Leonelli, and J. Smith.  A differential approach for staged trees.  In _European Conference on Symbolic and Quantitative Approaches to Reasoning and Uncertainty_ , pages 346–355. Springer, 2015.
  * Gurwicz and Lerner (2006) Y. Gurwicz and B. Lerner.  Bayesian class-matched multinet classifier.  In _Proceedings of the Joint IAPR International Conference on Structural, Syntactic, and Statistical Pattern Recognition_ , pages 145–153, 2006\.
  * Görgen et al. (2020) C. Görgen, M. Leonelli, and O. Marigliano.  The curved exponential family of a staged tree.  _Electronic Journal of Statistics_ , 16(1):2607–2620, 2020.
  * Ho (1995) T. Ho.  Random decision forests.  In _Proceedings of the 3rd International Conference on Document Analysis and Recognition_ , volume 1, pages 278–282, 1995.
  * Hruschka Jr and Ebecken (2007) E. R. Hruschka Jr and N. F. Ebecken.  Towards efficient variables ordering for Bayesian networks classifier.  _Data & Knowledge Engineering_, 63(2):258–269, 2007.
  * Huang et al. (2003) K. Huang, I. King, and M. Lyu.  Discriminative training of Bayesian Chow-Liu multinet classifiers.  In _Proceedings of the International Joint Conference on Neural Networks_ , pages 484–488, 2003.
  * Hussein and Santos (2004) A. Hussein and E. Santos.  Exploring case-based Bayesian networks and Bayesian multi-nets for classification.  In _Proceedings of the 17th Conference of the Canadian Society of Computational Studies of Intelligence_ , pages 485–492, 2004.
  * Jaeger et al. (2006) M. Jaeger, J. Nielsen, and T. Silander.  Learning probabilistic decision graphs.  _International Journal of Approximate Reasoning_ , 42(1-2):84–100, 2006.
  * Keeble et al. (2017) C. Keeble, P. Thwaites, P. Baxter, S. Barber, R. Parslow, and G. Law.  Learning through chain event graphs: The role of maternal factors in childhood type 1 diabetes.  _American Journal of Epidemiology_ , 186(10):1204–1208, 2017.
  * Keogh and Pazzani (2002) E. Keogh and M. Pazzani.  Learning the structure of augmented Bayesian classifiers.  _International Journal on Artificial Intelligence Tools_ , 11(4):587–601, 2002.
  * Leonelli (2019) M. Leonelli.  Sensitivity analysis beyond linearity.  _International Journal of Approximate Reasoning_ , 113:106–118, 2019.
  * Leonelli and Varando (2021) M. Leonelli and G. Varando.  Context-specific causal discovery for categorical data using staged trees.  _arXiv:2106.04416_ , 2021.
  * Leonelli and Varando (2022a) M. Leonelli and G. Varando.  Highly efficient structural learning of sparse staged trees.  _arXiv:2206.06970_ , 2022a.
  * Leonelli and Varando (2022b) M. Leonelli and G. Varando.  Structural learning of simple staged trees.  _arXiv:2203.04390_ , 2022b.
  * Leonelli et al. (2017) M. Leonelli, C. Görgen, and J. Smith.  Sensitivity analysis in multilinear probabilistic models.  _Information Sciences_ , 411:84–97, 2017.
  * Ling and Zhang (2002) C. X. Ling and H. Zhang.  The representational power of discrete Bayesian networks.  _Journal of Machine Learning Research_ , 3:709–721, 2002\.
  * Mihaljevic et al. (2018) B. Mihaljevic, C. Bielza, and P. Larrañaga.  _bnclassify: learning discrete Bayesian network classifiers from data_ , 2018.  R package version 0.4.0.
  * Minsky (1961) M. Minsky.  Steps towards artificial intelligence.  In _Computers and Thought_ , pages 406–450, 1961.
  * O’Donnell (2014) R. O’Donnell.  _Analysis of Boolean functions_.  Cambridge University Press, 2014.
  * Pensar et al. (2015) J. Pensar, H. Nyman, T. Koski, and J. Corander.  Labeled directed acyclic graphs: a generalization of context-specific independence in directed graphical models.  _Data Mining and Knowledge Discovery_ , 29(2):503–533, 2015.
  * Pensar et al. (2016) J. Pensar, H. Nyman, J. Lintusaari, and J. Corander.  The role of local partial independence in learning of Bayesian networks.  _International Journal of Approximate Reasoning_ , 69:91–105, 2016.
  * Poole and Zhang (2003) D. Poole and N. Zhang.  Exploiting contextual independence in probabilistic inference.  _Journal of Artificial Intelligence Research_ , 18:263–313, 2003.
  * Shafer (1996) G. Shafer.  _The art of causal conjecture_.  MIT Press, 1996.
  * Silander and Leong (2013) T. Silander and T.-Y. Leong.  A dynamic programming algorithm for learning chain event graphs.  In _Proceedings of the International Conference on Discovery Science_ , pages 201–216, 2013.
  * Smith and Anderson (2008) J. Smith and P. Anderson.  Conditional independence and chain event graphs.  _Artificial Intelligence_ , 172(1):42 – 68, 2008\.
  * Specht (1990) D. Specht.  Probabilistic neural networks.  _Neural Networks_ , 3(1):109–118, 1990.
  * Thwaites et al. (2010) P. Thwaites, J. Smith, and E. Riccomagno.  Causal analysis with chain event graphs.  _Artificial Intelligence_ , 174(12-13):889–909, 2010.
  * Varando et al. (2015) G. Varando, C. Bielza, and P. Larrañaga.  Decision boundary for discrete Bayesian network classifiers.  _Journal of Machine Learning Research_ , 16:2725–2749, 2015\.
  * Varando et al. (2016) G. Varando, C. Bielza, and P. Larrañaga.  Decision functions for chain classifiers based on Bayesian networks for multi-label classification.  _International Journal of Approximate Reasoning_ , 68:164–178, 2016.
  * Varando et al. (2021) G. Varando, F. Carli, and M. Leonelli.  Staged trees and asymmetry-labeled dags.  _arXiv:2108.01994_ , 2021.


# Time Adaptive Gaussian Model

Federico Ciech${}^{1},*$ and Veronica Tozzo2,3,∗

(1 Department of Bioengineering, Robotics, Informatics and System Engineering,
Università degli Studi di Genova
2 Center for System Biology and Department of Pathology, Massachusetts General
Hospital
3 Department of Systems Biology, Harvard Medical School
∗ These authors equally contributed to this paper. )

###### Abstract

Multivariate time series analysis is becoming an integral part of data
analysis pipelines. Understanding the individual time point connections
between covariates as well as how these connections change in time is non-
trivial. To this aim, we propose a novel method that leverages on Hidden
Markov Models and Gaussian Graphical Models — Time Adaptive Gaussian Model
(TAGM). Our model is a generalization of state-of-the-art methods for the
inference of temporal graphical models, its formulation leverages on both
aspects of these models providing better results than current methods. In
particular,it performs pattern recognition by clustering data points in time;
and, it finds probabilistic (and possibly causal) relationships among the
observed variables. Compared to current methods for temporal network
inference, it reduces the basic assumptions while still showing good inference
performances.

##  1 Introduction

The inference of temporal networks has started to become a common topic in the
last few years [10, 6, 12, 13, 14, 4, 23, 24]. Current method approach the
problem by taking a multi-variate time series and dividing it in chunks [8,
13, 24]. Each chunk is assumed to be a short enough period of time that all
its points are identically and independently sampled from a unique
distribution. Such distribution, when the variables are continuous, can be
represented as a Gaussian Graphical Model (GGM), where the conditional
independency patterns of variables are encoded as edges of graphs that evolve
in time. These approaches show good performances, but the assumption that the
time points in each chunk are i.i.d. is most often is not true.

In this paper, we propose Time Adaptive Gaussian Model (TAGM), a combination
of GGMs with Hidden Markov Models (HHMs) [1] that allows to easily consider
the sequence of single time points and relax the chunk assumption. It also
allows us to obtain clusters of time points as well as repeated evolving
patterns of graphs that may be impossible to obtain with current state-of-the-
art methods.

A schematic representation of the model is presented in Figure 1. Here, we are
considering an HMM at 2 states, which by looking at the left panel, presents
as the sequence 1, 2, 2, 1. Given the latent states, observations $x_{1}$ and
$x_{4}$ belong to an underlying distribution while $x_{2}$ and $x_{3}$ belong
to another one. Given the Markov chain that connects the latent state we can
assume that all these observations are independent and, thus, use them to
infer two GGMs that models the probability distribution of cluster 1 and
cluster 2 (right panel Figure 1). The inferred GGMs provide us more
information on how, within each latent state, the variables are dependent to
each other.

Figure 1: TAGM schematic representation. To each temporal observation is
associated a state. Each state is characterized by an underlying distribution
which is represented by a graphical model.

Our approach, given a multi-variate time-series, is able to

  1. 1.

cluster temporal data points considering sequentiality. Note that, TAGM does
not intend to group time series based on their morphology. Rather, it
_clusters_ the observations within the time series based on their similarities
and the order they appear. Such approach, in literature, is either formulated
as a standard clustering problem on the time points or as a longitudinal
clustering on the time-series [19, 7, 21].

  2. 2.

Infer temporal-dependent conditional dependencies among variables: TAGM
relaxed the chunk assumption common to the inference methods available in
literature [13, 24] thus inferring a time-varying network that adapts at each
observation. Note that TAGM assumes a sequentiality of the states, while in
[18, 7] the authors combined GGMs with Gaussian Mixture Models. This approach
is more prone in clustering the points in classes which can be seen as an
unsupervised extension of the Joint Graphical Lasso [6].

We show on synthetic data that the model performs better than current state-
of-the-art methods on both tasks (1) and (2). We argue that the general
characteristics of TAGM make it a trustful model that may be applied in a
variety of applicative domains.

##### Related work

For time series clustering, to our knowledge, HMM [1] is the only clustering
method which considers also sequentiality, as state-of-the-art methods
typically cluster points based on the feature differences [19, 7, 21]. The
inference of time-varying network has being tackled recently in literature [8,
13, 24] by dividing the time-series in chunks with the strong assumption that
all the time points within a chunk are i.i.d. TAGM overcomes this assumption
proposing a more elegant way of inferring evolving networks as well as their
pattern of evolution. By inferring $K$ different graphs we reach a deeper
level of understanding on the states that allows us to gain insights on the
system under analysis [14]. Lastly, state-of-the-art prediction methods on
time-series [22, 15, 17, 25] commonly assume the relations among the past
values of the variables and the present values of each variable to be constant
in time. An idea similar to TAGM was proposed with Gaussian Mixture Models
(GMMs) [7] where they combined GMM with GGMs [18]. The use of GMMs though
would not allow to explicitly consider sequentiality and it is therefore not
suited for the analysis of time-series. In literature, we found two examples
that explicitly consider non-stationarity and sequentiality in a setting
similar to ours.

##  2 Preliminaries

Consider a complex non-stationary system where, with non-stationarity, we
intend a change in time of the underlying distribution of observations. Of
this system we observe $N$ temporal observations, each of these observation
$n=1,\dots,N$ is a vector $\mathbf{x}_{n}\in\mathbb{R}^{d}$ of $d$ variables
sampled from an unknown distribution. In the following we will denote vectors
with bold letters $\mathbf{x}$, matrices with capital letters $X$ and
sequences of vectors or matrices with bold capital letters $\mathbf{X}$.

###  2.1 Hidden Markov Model

HMMs are statistical models widely applied sequential data structures. HMMs
assume that the series of observations is generated by a given number $K$ of
(hidden) internal states which follow a Markov process (see left panel of
Figure 1) [1]. Consider the $N$ sequential (temporal) observations, we pair
each of them with a hidden (latent) state
$\mathbf{z}_{n}=\sum_{i=1}^{K}\mathds{1}_{\\{i=k\\}}\mathbf{e}_{i}$, where
$\mathbf{e}_{i}$ is the $K$-dimensional natural basis which has a non-zero
component only at position $i$ and $k$ is the _cluster_ label of observation
$n$. We use the notation $z_{n,k}$ to indicate the $k$-th positional value of
the vector $\mathbf{z}_{n}$. Through their hidden states, the observations
$\mathbf{x}_{n}$ and $\mathbf{x}_{n+1}$ become independent given their states
(see Figure 1 left panel). The hidden states, on the other hand, follow a
Markov chain process which satisfies the conditional independence property
$\mathbf{z}_{n+1}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathbf{z}_{n-1}|\mathbf{z}_{n}$.
The HMM joint distribution on both the observations
$\mathbf{X}=(\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{N})$ and the
latent variables of
$\mathbf{Z}=(\mathbf{z}_{1},\mathbf{z}_{2},\dots,\mathbf{z}_{N})$ is then
given by

$p(\mathbf{X},\mathbf{Z}|\mathbf{\pi},A,\mathbf{\phi})=p(\mathbf{z}_{1}|\pi)\Big{[}\prod_{n=2}^{N}p(\mathbf{z}_{n}|\mathbf{z}_{n-1},A)\Big{]}\prod_{n=1}^{N}p(\mathbf{x}_{n}|\mathbf{z}_{n},\phi).$
(1)

The probability
$p(\mathbf{z}_{1}|\mathbf{\pi})=\prod_{k=1}^{K}\pi_{k}^{z_{1,k}}\quad\text{with}\quad\sum_{k}\pi_{k}=1$
is the initial latent node $\mathbf{z}_{1}$ probability, which differs from
the other states as there is no parent node. Thus, its marginal distribution
is embodied by a vector of probabilities $\mathbf{\pi}$ whose elements
$\pi_{k}\equiv p(z_{1,k}=1)$ represents the probability of the first
observation to belong to the $k$-th state.
The probability
$p(\mathbf{z}_{n}|\mathbf{z}_{n-1,A})=\prod_{k=1}^{K}\prod_{j=1}^{K}A_{j,k}^{z_{n-1,j}z_{n,k}}$
is the transition probability of moving from one state to the other where
$A\in[0,1]^{K\times K}$ is the transition matrix that we assume to be constant
in time.It is defined as $A_{j,k}=p(z_{n,k}=1|z_{n-1,j}=1)$ with $0\leq
A_{j,k}\leq 1$ and $\sum_{k}A_{j,k}=1$.
Lastly,
$p(\mathbf{x}_{n}|\mathbf{z}_{n},\mathbf{\phi})=\prod_{k=1}^{K}p(\mathbf{x}_{n}|\phi_{k})^{z_{n,k}}$,
are the emission probabilities where
$\mathbf{\phi}=\\{\phi_{1},\dots,\phi_{K}\\}$ is a set of $K$ different
parameters governing the distributions, one for each of the possible $K$
states.

###  2.2 Gaussian Graphical Models

GGMs are typically employed in the analysis of multivariate problems where one
seeks to understand the relationships among variables. A GGM is a probability
distribution which factorizes according to an undirected graph whose set of
edges univocally determines a multivariate normal distribution
$\mathcal{N}(\mu,\Sigma)$. Indeed, the _precision matrix_ ,
$\Theta=\Sigma^{-1}$ encodes the conditional independence between pairs of
variables, and $\Theta(i,j)=0$ implies the absence of an edge in the graph.
Thus, $\Theta$ is the weighted adjacency matrix of the graph (see right panel
Figure 1) [16].
Given $M$ observations of $d$ variables $\mathbf{X}\in\mathbb{R}^{M\times d}$
we aim at inferring the underlying graph corresponding to the precision matrix
$\Theta$. In order to perform such inference we need to assume that the
underlying graph is sparse. This is due to the combinatorial nature of the
problem that requires to be constrained to have identifiability guarantees
[9]. The most common way of inferring such graph is through the _Graphical
Lasso_ (GL) [9], a penalized Maximum Likelihood Estimation (MLE) method that
solves the following problem

$\underset{\Theta\succ 0}{\text{argmin}}~{}~{}\text{tr}(\Theta S)-\text{log
det}(\Theta)+\lambda\|\Theta\|_{1,od}$ (2)

where $S$ is the empirical covariance matrix of the input data defined as
$S=\frac{1}{n}X^{\top}X$, $\text{tr}(\Theta S)~{}-~{}\text{log det}(\Theta)$
is the negative log likelihood of the multivariate normal distribution and
$\|\cdot\|_{1,dot}$ is the off-diagonal $\ell_{1}$-norm that imposes sparsity
on the precision matrix $\Theta$ without considering the diagonal elements.

##  3 Time Adaptive Gaussian Model

TAGM presents an elegant way of inferring evolving networks as well as their
pattern of evolution. Consider the system presented in previous section, we
have $N$ temporal observations that are connected through $K$ latent states.
Each of this state is connected to an unknown distribution that we want to
model as a graphical model. By inferring $K$ different graphs we reach a
deeper level of understanding on the intra-dependency patterns of the states
that allows us to gain insights on the system under analysis [14].

TAGM is a combination of HHMs with GGMs, which can be straightforwardly
derived by assuming that all the observations belonging to a cluster $k$ are
i.i.d. and multivariate normally distributed. Thus, we assume that the
emission probabilities are

$p(\mathbf{x}_{n}|\mathbf{z}_{n},\mathbf{\phi})=\prod_{k=1}^{K}\mathcal{N}(\mathbf{x}_{n}|\mu_{k},\Theta_{k}^{-1})^{z_{n,k}}$

in such a way to have an explicit correspondence between the distribution of
each state of the HMM and a graph, modeled through the precision matrix
$\Theta$. Thus, the parameters $\phi$ are equal to the sequences
$\mu=(\mu_{1},\dots,\mu_{K})$ and
$\mathbf{\Theta}=(\Theta_{1},\dots,\Theta_{K})$.
In order to obtain the Graphical Lasso form we add a sparsity constraint by
multiplying the joint probability distribution in Equation (1) with a
Laplacian prior on the precision matrix. The Laplacian prior would also
provide more interpretable results as we keep only the strongest edges of the
final graphs. The posterior is defined as

$\displaystyle
p(\mathbf{\pi},A,\mathbf{\mu},\mathbf{\Theta}|\mathbf{X},\mathbf{Z})\propto$
$\displaystyle
p(\mathbf{z}_{1}|\pi)\Big{[}\prod_{n=2}^{N}p(\mathbf{z}_{n}|\mathbf{z}_{n-1},A)\Big{]}$
$\displaystyle\prod_{n=1}^{N}\prod_{k=1}^{K}\mathcal{N}(\mathbf{x}_{n}|\mu_{k},\Theta_{k}^{-1})^{z_{n,k}}e^{-\frac{\lambda}{2}||\Theta_{k}||_{1,od}}.$

where with the notation $\mathbf{z}_{n}$ and $\mathbf{z}_{n-1}$ we hid a
product on $z_{n,k}$ and $z_{n-1,j}$ for $k,j=1,\dots,K$. In order to fit such
distribution on data we need to detect the best set of parameters
$\theta=\\{\pi,A,\mathbf{\mu},\mathbf{\Theta}\\}$ through a Maximum A
Posteriori approach. To marginalize over the latent variables
$p(\mathbf{X}|\mathbf{\theta})=\sum_{\mathbf{Z}}p(\mathbf{X},\mathbf{Z}|\mathbf{\theta})$
we employ the _expectation maximization_ (EM) algorithm, in its particular
version known as Baum’s algorithm [2]. The _expectation maximization_ (EM)
algorithm [2] is an iterative algorithm which alternates two steps: the _E-
step_ and the _M-step_. It starts with some initial selection for the model
parameters, which we denote by $\theta_{\text{old}}$. Then from the posterior
distribution over the latent variables
$p(\mathbf{X},\mathbf{Z}|\mathbf{\theta}_{\text{old}})$ it evaluates the
expectation of the complete-data log-likelihood as a function of the
parameters $\theta=\\{\mathbf{\pi},A,\mathbf{\mu},\mathbf{\Theta}\\}$:

$\displaystyle Q(\theta,\theta_{\text{old}})=$
$\displaystyle\gamma(\mathbf{z}_{n})\ln(\mathbf{\pi})+\sum_{n=2}^{N}\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})\ln
A+$
$\displaystyle\sum_{n=1}^{N}\sum_{k=1}^{K}\left[\gamma(z_{n,k})\ln\mathcal{N}(\mathbf{x}_{n}|\mu_{k},\Theta_{k}^{-1})-\frac{1}{2}\lambda||\Theta_{k}||_{1,od}\right]$

where $\gamma(z_{n,k})$ and $\xi(z_{n-1,j},z_{n,k})$ come from the definition

$\displaystyle\gamma(\mathbf{z}_{n})$
$\displaystyle=p(\mathbf{z}_{n}|\mathbf{X},\mathbf{\theta}_{\text{old}})$
$\displaystyle\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})$
$\displaystyle=p(\mathbf{z}_{n-1},\mathbf{z}_{n}|\mathbf{X},\mathbf{\theta}_{\text{old}}).$

Denoting $\gamma(z_{n,k})$ the conditional probability of $z_{n,k}=1$, with a
similar use of notation for $\xi(z_{n-1,j},z_{n,k})$. Because the expectation
of a binary random variable is just the probability that it takes the value 1,
we have

$\displaystyle\gamma(z_{n,k})$
$\displaystyle=\mathbb{E}[z_{n,k}]=\sum_{\mathbf{z}}\gamma(\mathbf{z})z_{n,k}$
$\displaystyle\xi(z_{n-1,j},z_{n,k})$
$\displaystyle=\mathbb{E}[z_{n-1,j}z_{n,k}]=\sum_{\mathbf{z}}\gamma(\mathbf{z})z_{n-1,k}z_{n,k}.$

The goal of the _E-step_ is to evaluate the quantities
$\gamma(\mathbf{z}_{n})$ and $\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})$. While the
_M step_ maximizes $Q(\theta,\theta_{\text{old}})$ with respect to the
parameters $\mathbf{\theta}$ in which we treat $\gamma(\mathbf{z}_{n})$ and
$\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})$ as constants.

###  3.1 E step

An efficient algorithm to evaluate the quantities $\gamma(\mathbf{z}_{n})$ and
$\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})$ is the forward-backward algorithm
described in [2]. Here we want to emphasize the most relevant formulas which
characterizes the algorithm. If the reader is interested in more details and
how to derive them we refer to [3].

Using Bayes’ theorem, we have

$\gamma(\mathbf{z}_{n})=\frac{p(\mathbf{X}|\mathbf{z}_{n})p(\mathbf{z}_{n})}{p(\mathbf{X})}=\frac{p(\mathbf{x}_{1},\dots,\mathbf{x}_{n},\mathbf{z}_{n})p(\mathbf{x}_{n+1},\dots,\mathbf{x}_{N}|\mathbf{z}_{n})}{p(\mathbf{X})}=\frac{\alpha(\mathbf{z}_{n})\beta(\mathbf{z}_{n})}{p(\mathbf{X})}.$
(3)

where we have defined

$\displaystyle\alpha(\mathbf{z}_{n})$
$\displaystyle=p(\mathbf{x}_{1},\dots,\mathbf{x}_{n},\mathbf{z}_{n})$
$\displaystyle\beta(\mathbf{z}_{n})$
$\displaystyle=p(\mathbf{x}_{n+1},\dots,\mathbf{x}_{N}|\mathbf{z}_{n}).$

$\alpha(\mathbf{z}_{n})$ is also called _forward process_ while
$\beta(\mathbf{z}_{n})$ _backward process_. It is possible to prove that by
the conditional independence assumption
$\mathbf{z}_{n+1}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathbf{z}_{n-1}|\mathbf{z}_{n}$ the forward
and backward processes are characterized by two recursive equations:

$\displaystyle\alpha(\mathbf{z}_{n})$
$\displaystyle=p(x_{n}|\mathbf{z}_{n})\sum_{\mathbf{z}_{n-1}}\alpha(\mathbf{z}_{n-1})p(\mathbf{z}_{n}|\mathbf{z}_{n-1}),$
(4) $\displaystyle\beta(\mathbf{z}_{n})$
$\displaystyle=\sum_{\mathbf{z}_{n+1}}\beta(\mathbf{z}_{n+1})p(\mathbf{x}_{n+1}|\mathbf{z}_{n+1})p(\mathbf{z}_{n+1}|\mathbf{z}_{n}),$
(5)

with the initial conditions:

$\displaystyle\alpha(\mathbf{z}_{1})$
$\displaystyle=\prod_{k=1}^{K}\\{\pi_{k}\mathcal{N}(\mathbf{x}_{1}|\mu_{k},\Theta^{-1}_{k})\\}^{z_{1,k}}$
$\displaystyle\beta(\mathbf{z}_{N})$ $\displaystyle=(1,\dots,1).$

If we sum both sides of (3) over $\mathbf{z}_{n}$, and use the fact that the
left-hand side is a normalized distribution, we obtain

$p(\mathbf{X})=\sum_{\mathbf{z}_{n}}\alpha(\mathbf{z}_{n})\beta(\mathbf{z}_{n})=\sum_{\mathbf{z}_{N}}\alpha(\mathbf{z}_{N}).$

Using similar arguments we have

$\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})=\frac{\alpha(\mathbf{z}_{n-1})p(\mathbf{x}_{n}|\mathbf{z}_{n})p(\mathbf{z}_{n}|\mathbf{z}_{n-1})\beta(\mathbf{z}_{n})}{p(\mathbf{X})}.$

For moderate lengths of chain the calculation of $\alpha(\mathbf{z})$ can go
to zero exponentially quickly. We therefore work with re-scaled versions of
$\alpha(\mathbf{z})$ and $\beta(\mathbf{z})$ whose values remain of order
unity. The corresponding scaling factors cancel out when we use there re-
scaled quantities in the EM algorithm.

We define a normalised version of $\alpha(\mathbf{z})$ as

$\hat{\alpha}(\mathbf{z}_{n})=p(\mathbf{z}_{n}|\mathbf{x}_{1},\dots,\mathbf{x}_{n})=\frac{\alpha(\mathbf{z}_{n})}{p(\mathbf{x}_{1},\dots,\mathbf{x}_{n})}$

which we expect to be well behaved numerically because it is a probability
distribution over $K$ variables for any value of $n$. In order to relate to
the original $\alpha(\mathbf{z})$ variables we introduce scaling factors

$c_{n}=p(\mathbf{x}_{n}|\mathbf{x}_{1},\dots,\mathbf{x}_{n-1})$

and therefore

$p(\mathbf{x}_{1},\dots,\mathbf{x}_{n})=\prod_{m=1}^{n}c_{m}.$

From the $\alpha$ and $\beta$ recursive equations (4) and (5) its scaled
formula $\hat{\alpha}(\mathbf{z})$ becomes

$\hat{\alpha}(\mathbf{z}_{n})=\frac{p(\mathbf{x}_{n}|\mathbf{z}_{n})\sum_{\mathbf{z}_{n-1}}\hat{\alpha}(\mathbf{z}_{n-1})p(\mathbf{z}_{n}|\mathbf{z}_{n-1})}{c_{n}}$

and, similarly

$\hat{\beta}(\mathbf{z}_{n})=\frac{\sum_{\mathbf{z}_{n+1}}\hat{\beta}(\mathbf{z}_{n+1})p(\mathbf{x}_{n+1}|\mathbf{z}_{n+1})p(\mathbf{z}_{n+1}|\mathbf{z}_{n})}{c_{n+1}}.$

Note that for the computation of $\beta$s we can recur to the scaling factors
we computed in the $\alpha$ phase. We can also notice that the probability
distribution of $\mathbf{X}$ becomes

$p(\mathbf{X})=\prod_{n=1}^{N}c_{n}$

and that

$\displaystyle\gamma(\mathbf{z}_{n})$
$\displaystyle=\hat{\alpha}(\mathbf{z}_{n})\hat{\beta}(\mathbf{z}_{n})$
$\displaystyle\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})$
$\displaystyle=c^{-1}_{n}\hat{\alpha}(\mathbf{z}_{n})p(\mathbf{x}_{n}|\mathbf{z}_{n})p(\mathbf{z}_{n}|\mathbf{z}_{n-1})\hat{\beta}(\mathbf{z}_{n}).$

###  3.2 M step:

Given $\gamma(\mathbf{z}_{n})$ and $\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})$
computed in the E step, the M step finds the optimal parameters
$\mathbf{\theta}$. To maximize respect to $\pi$ and $\mathbf{A}$ we keep the
addends of $Q(\theta,\theta_{\text{old}})$ which are directly interested in,
where we use the appropriate Lagrangian to take into account of the
constraints.

$\displaystyle\frac{\partial Q(\pi_{k})}{\partial\pi_{k}}$
$\displaystyle=\frac{\partial\Big{(}\sum_{k=1}^{K}\gamma(z_{1,k})\ln(\pi_{k})+\lambda(1-\sum_{k=1}^{K}\pi_{k})\Big{)}}{\partial\pi_{k}}=\frac{\gamma(z_{1,k})}{\pi_{k}}-\lambda=0$

We multiply both sides by $\pi_{k}$ and summing over all $k\in\\{1,\dots,K\\}$
and obtain $\lambda=\sum_{j=1}^{K}\gamma(z_{1,k})$ using that
$\sum_{j=1}^{K}\pi_{j}=1$. Therefore

$\pi_{k}=\frac{\gamma(z_{1,k})}{\sum_{j=1}^{K}\gamma(z_{1,k})}.$

Analogously for $\mathbf{A}$

$\displaystyle\frac{\partial Q(A_{j,k})}{\partial A_{j,k}}$
$\displaystyle=\frac{\partial\Big{(}\sum_{n=2}^{N}\sum_{j=1}^{K}\sum_{k=1}^{K}\xi(z_{n-1,j},z_{n,k})\ln
A_{j,k}+\lambda(1-\sum_{k=1}^{K}A_{j,k})\Big{)}}{\partial A_{j,k}}$
$\displaystyle=\frac{\sum_{n=2}^{N}\xi(z_{n-1,j},z_{n,k})}{A_{j,k}}-\lambda=0$

Similarly as before we multiply both sides by $A_{j,l}$ and summing over all
$l\in\\{1,\dots,K\\}$ and obtain
$\lambda=\sum_{l=1}^{K}\sum_{n=2}^{N}\xi(z_{n-1,j},z_{n,l})$ using that
$\sum_{l=1}^{K}A_{j,l}=1$. Therefore

$A_{j,k}=\frac{\sum_{n=2}^{N}\xi(z_{n-1,j},z_{n,k})}{\sum_{l=1}^{K}\sum_{n=2}^{N}\xi(z_{n-1,j},z_{n,l})}.$

When we differentiate by $\mu_{k}$ we have the same equation as in the
Gaussian mixture model, therefore we give directly the result

$\mu_{k}=\frac{\sum_{n=1}^{N}\gamma(z_{n,k})\mathbf{x}_{n}}{\sum_{n=1}^{N}\gamma(z_{n,k})}.$

The maximization of $\Theta_{k}$ is performed in the main paper in Section 3.

For a detailed explanation of all the properties of this algorithm we refer to
[20]. The E step corresponds to evaluating the expected value of the log-
likelihood at $\mathbf{\theta}$. Given the values $\theta_{\text{old}}$, i.e.,
the parameter values at the previous iteration, the expectation is defined by
the function

$\displaystyle\begin{split}Q(\theta,\theta_{\text{old}})&=\gamma(\mathbf{z}_{n})\ln(\mathbf{\pi})+\sum_{n=2}^{N}\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})\ln
A\\\
&+\sum_{n=1}^{N}\sum_{k=1}^{K}\left[\gamma(z_{n,k})\ln\mathcal{N}(\mathbf{x}_{n}|\mu_{k},\Theta_{k}^{-1})-\frac{1}{2}\lambda||\Theta_{k}||_{1,od}\right]\end{split}$
(6)

In the functional,
$\gamma(\mathbf{z}_{n})=p(\mathbf{z}_{n}|\mathbf{X},\mathbf{\theta}_{\text{old}})$
and
$\xi(\mathbf{z}_{n-1},\mathbf{z}_{n})=p(\mathbf{z}_{n-1},\mathbf{z}_{n}|\mathbf{X},\mathbf{\theta}_{\text{old}})$
are the expectations on the latent variables and they are typically computed
by a _forward-backward_ algorithm [3].
Once the computation of the expectation is performed, the M step consists in
finding the $\mathbf{\theta}$ values which maximize the function $Q$. The
derivation of the new parameters $\pi,A$ and $\mu$ is straightforward from
literature [3]. The maximization of Equation (6) w.r.t. $\mathbf{\Theta}$,
instead requires further attention and it can be shown that, given the
imposition of the Laplacian prior, reduces to the Graphical Lasso (Equation
(2)) with few algebraic manipulations. Furthermore, the maximization of
$\mathbf{\Theta}$ can be performed separately for each $\Theta_{k}$ as, given
$k$, all the maximizations are independent between each other. Thus, if we
indicate as $\bar{\mathbf{x}}_{i}=\mathbf{x}_{i}-\mathbf{\mu}_{k}$ the
centered observations belonging to cluster $k$ the joint log-likelihood at
fixed $k$ writes out as

$\displaystyle\sum_{n=1}^{N}\gamma(z_{n,k})\ln\mathcal{N}(\mathbf{x}_{n}|\mu_{k},\Theta_{k}^{-1})-\frac{1}{2}\lambda||\Theta_{k}||_{1,od}$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\frac{\gamma(z_{n,k})}{2}\ln\text{det}\Theta_{k}-\frac{1}{2}\text{tr}\Big{(}\sum_{n=1}^{N}\gamma(z_{n,k})(\bar{\mathbf{x}}_{n}\bar{\mathbf{x}}_{n}^{\prime})\Theta_{k}\Big{)}-\frac{1}{2}\lambda||\Theta_{k}||_{1,od}$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\frac{\gamma(z_{n,k})}{2}\Big{[}\ln\text{det}\Theta_{k}-\text{tr}\Big{(}\tilde{S}_{k}\Theta_{k}\Big{)}-\tilde{\lambda}_{k}||\Theta_{k}||_{1,od}\Big{]}$
(7)

where
$\tilde{S}_{k}=\frac{\sum_{n=1}^{N}\gamma(z_{n,k})(\bar{\mathbf{x}}_{n},\bar{\mathbf{x}}_{n}^{\prime})}{\sum_{n=1}^{N}\gamma(z_{n,k})}$
is the weighted empirical covariance matrix,
$\tilde{\lambda}_{k}=\frac{\lambda}{\sum_{n=1}^{N}\gamma(z_{n,k})}$ and

$\Theta_{k}=\arg\max_{\Theta_{k}}\Big{\\{}\ln\text{det}\Theta_{k}-\text{tr}\Big{(}\tilde{S}_{k}\Theta_{k}\Big{)}-\tilde{\lambda}_{k}||\Theta_{k}||_{1,od}\Big{\\}}$
(8)

which is equivalent to a Graphical Lasso [9]. Equation (8) is a convex
functional having guarantees of reaching a global optimum. Differently,
Equation (6) is non-convex and depending on the initialization different local
optima may be reached.

##  4 TAGM extensions

The proposed approach can benefit from two types of extensions. The first one
consists in higher-order Markov relationships among latent states. This can be
achieved following [11], with the main drawback of a much higher computational
time for learning as the transition matrix and the corresponding initial state
dimensions increase. The second extension consists in an online learning
version that will allow to employ the model in a more applicative setting
where we may deal with high-frequency data [5]. Indeed, in a situation where
new observations arrive at a high rate (every second or even millisecond) we
want to be able to fine tune the model online in order to consider such
observations instantaneously. This will allow to promptly gain insights on
data and possibly predict the next time point. The weakness of this extension
is that it requires an approximation which makes the updated parameters less
accurate respect to the batch (original) version. We discuss in further
details these two extensions in Sections 8 and 9 of the Supplementary
material. In particular, we present two experiments where we compare them with
TAGM that show the potential weaknesses.

###  4.1 Higher order extension: Memory Time Adaptive Gaussian model (MemTAGM)

Sometimes real world applications have events which rely on their past
realizations. Therefore we can exploit more information from data if we
consider a higher-order Markov process whose $\mathbf{z}_{n}$ state
probability does not depend only on $\mathbf{z}_{n-1}$ but also on the other
$r$ past states according to the choice of $r$. TAGM can be extended to higher
order sequential relationships. We consider a homogeneous Markov process of
order $r\in\mathbb{Z}^{+}$ over a finite state set $\\{1,\dots,K\\}$ with
hidden sequence $\\{\mathbf{z}\\}_{n=1}^{N}$. This stochastic process
satisfies

$p(\mathbf{z}_{n}|\\{\mathbf{z}_{\ell}\\}_{\ell<n})=p(\mathbf{z}_{n}|\\{\mathbf{z}_{\ell}\\}_{\ell=n-r}^{n-1})$

or in other words $\mathbf{z}_{n}$ can depend on a different number of hidden
past states, and we assume that the process is homogeneous i.e., the
transition probability is independent of $n$. To be as more general as
possible we allow that the emission probability of $\mathbf{x}_{n}$ can depend
not only on $\mathbf{z}_{n}$ but also from the previous $m\in\mathbb{Z}^{+}$
sequence of states

$p(\mathbf{x}_{n}|\\{\mathbf{x}_{\ell}\\}_{\ell<n},\\{\mathbf{z}_{\ell}\\}_{\ell\leq
n})=p(\mathbf{x}_{n}|\\{\mathbf{z}_{\ell}\\}_{\ell=n-(m-1)}^{n}).$

Each observation is conditionally independent of the previous ones and of the
state sequence history, given the current and the preceding $m-1$ states.

The idea is transform the High order hidden Markov model (HHMM) to a first
order hidden Markov model (HMM). It can be done by considering the following
two propositions where we omit the prove but it can be found in [11]

###### Proposition 1

Let
$\mathbf{Z}_{n}=[\mathbf{z}_{n},\mathbf{z}_{n-1},\dots,\mathbf{z}_{n-(\nu-1)}]^{\top}$.
The process $\\{\mathbf{Z}_{n}\\}$ is a first order homogeneous Markov process
for any $\nu\geq r$, taking values in $\mathcal{S}^{\nu}$.

###### Proposition 2

Let $\nu=\max\\{r,m\\}$. The state sequence $\\{\mathbf{Z}_{n}\\}$ and the
observation sequence $\\{\mathbf{x}_{n}\\}$ satisfy

$p(\mathbf{x}_{n}|\\{\mathbf{x}_{\ell}\\}_{\ell<n},\\{\mathbf{z}_{\ell}\\}_{\ell\leq
n})=p(\mathbf{x}_{n}|\mathbf{Z}_{n})$

and thus constitute a first order HMM.

We can thus reformulate HHMM as a first order HMM with $K^{\nu}$ states, where
$\nu=\max\\{r,m\\}$. Note that the last $\nu-1$ entries of $\mathbf{Z}_{n}$
are equal to the first $\nu-1$ entries of $\mathbf{Z}_{n+1}$, one concludes
that a transition from $\mathbf{z}_{n}$ to $\mathbf{z}_{n+1}$ is possible only
if
$\lfloor\mathbf{z}_{n}/K\rfloor=\mathbf{z}_{n+1}-\lfloor\mathbf{z}_{n+1}/K^{\nu-1}\rfloor
K^{\nu-1}$, and thus

$A_{i,j}=0\quad\text{if
}\Big{\lfloor}\frac{i}{K}\Big{\rfloor}\not=j-\Big{\lfloor}\frac{1}{K^{\nu-1}}\Big{\rfloor}K^{\nu-1}.$

Therefore we can use the EM algorithm to find the optimal parameters changing
the number of states to $K^{\nu}$. The $z_{n,j}$s which contribute to the M
step for a quantity of state $i$ are given by the set

$\mathcal{I}_{m}(i)=\Bigg{\\{}\Big{\lfloor}\frac{i}{K^{\nu-m}}\Big{\rfloor}K^{\nu-m},\Big{\lfloor}\frac{i}{K^{\nu-m}}\Big{\rfloor}K^{\nu-m}+1,\dots,\Big{(}\Big{\lfloor}\frac{i}{K^{\nu-m}}\Big{\rfloor}+1\Big{)}K^{\nu-m}-1\Bigg{\\}}.$

Therefore the means become

$\mu_{i}=\frac{\sum_{n=1}^{N}\mathbf{x}_{n}\sum_{j\in\mathcal{I}_{m}(i)}\gamma(z_{n,j})}{\sum_{n=1}^{N}\sum_{j\in\mathcal{I}_{m}(i)}\gamma(z_{n,j})}$

The empirical covariances to substitute in the graphical lasso equation is

$\tilde{S}_{i}=\frac{\sum_{n=1}^{N}(\mathbf{x}_{n}-\mathbf{\mu}_{i})(\mathbf{x}_{n}-\mathbf{\mu}_{i})^{\top}\sum_{j\in\mathcal{I}_{m}(i)}\gamma(z_{n,j})}{\sum_{n=1}^{N}\sum_{j\in\mathcal{I}_{m}(i)}\gamma(z_{n,j})},$

with the hyper-parameter
$\tilde{\lambda}_{k}=\frac{\lambda}{\sum_{n=1}^{N}\sum_{j\in\mathcal{I}_{m}(i)}\gamma(z_{n,j})}$.
The transition probability matrix becomes

$A_{i,j}=\begin{cases}&\frac{\sum_{n=2}^{N}\sum_{k\in\mathcal{I}_{r}(i)}\xi(z_{n-1,k},z_{n,\lfloor\frac{k}{K}\rfloor+\lfloor\frac{j}{K^{\nu-1}}\rfloor
K^{\nu-1}})}{\sum_{l=1}^{K^{\nu}}\sum_{n=2}^{N}\sum_{k\in\mathcal{I}_{r}(i)}\xi(z_{n-1,k},z_{n,l})}\quad\text{if
}\lfloor\frac{i}{K}\rfloor=j-\lfloor\frac{j}{K^{\nu-1}}\rfloor K^{\nu-1},\\\
&0\quad\text{otherwise}.\end{cases}$

Finally, the initial state probabilities are given by

$\pi_{i}=\gamma(z_{1,i}).$

Therefore if we increase the number of states to $K^{\nu}$ and modify the TAGM
M step formulas with the one just found we obtain the MemTAGM. We test MemTAGM
performance in subsection 9.2 where we compare it with TAGM.

###  4.2 On-line learning: Incremental Time Adaptive Gaussian model (IncTAGM)

In many applications it is important to update the TAGM parameters almost
instantaneously every time a new observation comes up. Since TAGM does not
allow to have such a quick response, it can be extended to an incremental
version. We call this model IncTAGM and it initially starts as a standard TAGM
reading a set of observations and updating its current parameters
$\pi,A,\mu,\mathbf{\Theta}$ according to new incoming data. Therefore, after
the standard TAGM has finished training on its observation set, it calculates
the revised $\alpha,\beta,\xi$ and $\gamma$ variables based on the new set of
observations. To update the model incrementally, we need recursive equations
for $\alpha$ and $\beta$ which depend on past values. Notice that the $\alpha$
recursive equation is already of this form. While the $\beta$ recursive
equation needs an approximation to become of that form. In fact if we assume
that $\beta(z_{T,i})\simeq\beta(z_{T,j})$ for every $i\not=j$, the $\beta$
recursive equation becomes

$\beta(\mathbf{z}_{T+1})=\frac{\beta(\mathbf{z}_{T})}{\sum_{\mathbf{z}_{T+1}}p(\mathbf{x}_{T+1}|\mathbf{z}_{T+1})p(\mathbf{z}_{T+1}|\mathbf{z}_{T})}.$
(9)

The M step optimal parameters are updated in the following way:

the initial state $\pi$

$\pi_{k}^{\prime}=\gamma(z_{1,k}),$ (10)

the transition matrix $\mathbf{A}$

$\displaystyle A^{T+1}_{j,k}$
$\displaystyle=\frac{\sum_{n=2}^{T}\xi(z_{n-1,j},z_{n,k})+\xi(z_{T,j},z_{T+1,k})}{\sum_{l=1}^{K}\sum_{n=2}^{T}\xi(z_{n-1,j},z_{n,l})+\sum_{l=1}^{K}\xi(z_{T,j},z_{T+1,l})}$
$\displaystyle=\frac{\sum_{n=2}^{T}\xi(z_{n-1,j},z_{n,k})}{\sum_{n=2}^{T+1}\gamma(z_{n-1,j})}+\frac{\xi(z_{T,j},z_{T+1,k})}{\sum_{n=2}^{T+1}\gamma(z_{n-1,j})}$
$\displaystyle=\frac{\sum_{n=2}^{T}\gamma(z_{n-1,j})}{\sum_{n=2}^{T+1}\gamma(z_{n-1,j})}\frac{\sum_{n=2}^{T}\xi(z_{n-1,j},z_{n,k})}{\sum_{n=2}^{T}\gamma(z_{n-1,j})}+\frac{\xi(z_{T,j},z_{T+1,k})}{\sum_{n=2}^{T+1}\gamma(z_{n-1,j})}$
$\displaystyle=\frac{\sum_{n=2}^{T}\gamma(z_{n-1,j})}{\sum_{n=2}^{T+1}\gamma(z_{n-1,j})}A^{T}_{j,k}+\frac{\xi(z_{T,j},z_{T+1,k})}{\sum_{n=2}^{T+1}\gamma(z_{n-1,j})}.$
(11)

Note that if we sum by $k\in\\{1,\dots,K\\}$ the row normalization holds.
Similarly we obtain the formula for the means

$\mu_{k}^{T+1}=\frac{\sum_{n=1}^{T}\gamma(z_{n,k})}{\sum_{n=1}^{T+1}\gamma(z_{n,k})}\mu_{k}^{T}+\frac{\gamma(z_{T+1,k})\mathbf{x}_{T+1}}{\sum_{n=1}^{T+1}\gamma(z_{n,k})}$
(12)

and the empirical covariances

$\tilde{S}^{T+1}_{k}=\frac{\sum_{n=1}^{T}\gamma(z_{n,k})}{\sum_{n=1}^{T+1}\gamma(z_{n,k})}\tilde{S}^{T}_{k}+\frac{\gamma(z_{T+1,k})(\mathbf{x}_{T+1}-\mathbf{\mu}^{T+1})(\mathbf{x}_{T+1}-\mathbf{\mu}^{T+1})^{\top}}{\sum_{n=1}^{T+1}\gamma(z_{n,k})},$
(13)

with the hyper-parameter
$\tilde{\lambda}_{k}=\frac{\lambda}{\sum_{n=1}^{T+1}\gamma(z_{n,k})}$.

####  4.2.1 Slide Incremental Time Adaptive Gaussian model (S-IncTAGM)

The training of new data points results in the accumulation of an increasingly
large observation set. As a result, if the time sequence considered is large
and the first point is far away in the past respect to the last point it is
possible that the initial trained observation points become outdated after
many updates and therefore they do not carry any useful information to analyze
the more recent points. Therefore, the new addition to the IncTAGM is a fixed
sliding window to effectively analyze discrete data (appropriately discarding
the outdated observations) whilst updating its model parameters.

The estimation of the $\alpha,\beta,\xi$ and $\gamma$ variables remains the
same as in the IncTAGM algorithm. What changes are the $A$, $\mathbf{\mu}$ and
$\mathbf{\Theta}$ updates. Using the simple moving average (SMA) definition

$sma=\frac{x_{1}+x_{2}\dots+x_{n}+x_{n+1}-x_{1}}{n}=ave+\frac{x_{n+1}}{n}-\frac{x_{1}}{n}$
(14)

where $ave=\frac{x_{1}+x_{2}\dots+x_{n}}{n}$, we update $A$, $\mathbf{\mu}$
and $\mathbf{\Theta}$ in the following way

$\displaystyle A^{T+1}_{j,k}$
$\displaystyle=\frac{\sum_{n=3}^{T+1}\xi(z_{n-1,j},z_{n,k})}{\sum_{n=3}^{T+1}\gamma(z_{n-1,j})},$
(15) $\displaystyle\mu_{k}^{T+1}$
$\displaystyle=\frac{\sum_{n=2}^{T+1}\gamma(z_{n,k})\mathbf{x}_{n}}{\sum_{n=2}^{T+1}\gamma(z_{n,k})},$
(16) $\displaystyle\tilde{S}^{T+1}_{k}$
$\displaystyle=\frac{\sum_{n=2}^{T+1}\gamma(z_{n,k})(\mathbf{x}_{n}-\mathbf{\mu}^{T+1})(\mathbf{x}_{n}-\mathbf{\mu}^{T+1})^{\top}}{\sum_{n=2}^{T+1}\gamma(z_{n,k})}$
(17)

with $\tilde{\lambda}_{k}=\frac{\lambda}{\sum_{n=2}^{T+1}\gamma(z_{n,k})}$.

##  5 Conclusions

We presented a novel methodology to perform data-mining and forecasting on
multi-variate time-series. Our method, namely TAGM, combines HMMs and GGMs,
providing a way to simultaneously cluster non-stationary time-series into
stationary sub-groups and for each cluster detecting probability relationships
among variables through graphical model inference. This simultaneous inference
is suitable to be transformed into a time-varying regression model that allows
to make predictions on non-stationary time-series. The coupling we performed
allows to generalize many state-of-the-art methods and provide a wide range of
analysis type to be performed on the time series.

## References

  * [1] Leonard E Baum and Ted Petrie.  Statistical inference for probabilistic functions of finite state markov chains.  The annals of mathematical statistics, 37(6):1554–1563, 1966.
  * [2] Leonard E Baum, Ted Petrie, George Soules, and Norman Weiss.  A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains.  The annals of mathematical statistics, 41(1):164–171, 1970.
  * [3] Christopher M Bishop.  Pattern recognition and machine learning.  springer, 2006.
  * [4] Andersen Chang, Tianyi Yao, and Genevera I Allen.  Graphical models and dynamic latent factors for modeling functional brain connectivity.  In 2019 IEEE Data Science Workshop (DSW), pages 57–63. IEEE, 2019\.
  * [5] Tiberiu Chis and Peter G Harrison.  Adapting hidden markov models for online learning.  Electronic Notes in Theoretical Computer Science, 318:109–127, 2015\.
  * [6] Patrick Danaher, Pei Wang, and Daniela M Witten.  The joint graphical lasso for inverse covariance estimation across multiple classes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(2):373–397, 2014.
  * [7] Brian S Everitt.  Finite mixture distributions.  Wiley StatsRef: Statistics Reference Online, 2014.
  * [8] Nicholas J Foti, Rahul Nadkarni, AK Lee, and Emily B Fox.  Sparse plus low-rank graphical models of time series for functional connectivity in meg.  In 2nd KDD Workshop on Mining and Learning from Time Series, 2016\.
  * [9] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9(3):432–441, 2008.
  * [10] Jian Guo, Elizaveta Levina, George Michailidis, and Ji Zhu.  Joint estimation of multiple graphical models.  Biometrika, 98(1):1–15, 2011.
  * [11] Uri Hadar et al.  High-order hidden markov models-estimation and implementation.  In 2009 IEEE/SP 15th Workshop on Statistical Signal Processing, pages 249–252. IEEE, 2009.
  * [12] David Hallac, Jure Leskovec, and Stephen Boyd.  Network lasso: Clustering and optimization in large graphs.  In Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, pages 387–396. ACM, 2015.
  * [13] David Hallac, Youngsuk Park, Stephen Boyd, and Jure Leskovec.  Network inference via the time-varying graphical lasso.  In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 205–213. ACM, 2017.
  * [14] David Hallac, Sagar Vare, Stephen Boyd, and Jure Leskovec.  Toeplitz inverse covariance-based clustering of multivariate time series data.  In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 215–223. ACM, 2017.
  * [15] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [16] Steffen L Lauritzen.  Graphical models, volume 17.  Clarendon Press, 1996.
  * [17] F. Li, L. Zhang, B. Chen, D. Gao, Y. Cheng, X. Zhang, Y. Yang, K. Gao, Z. Huang, and J. Peng.  A Light Gradient Boosting Machine for Remainning Useful Life Estimation of Aircraft Engines.  In 2018 21st International Conference on Intelligent Transportation Systems (ITSC), pages 3562–3567, 2018.
  * [18] Anani Lotsi and Ernst Wit.  High dimensional sparse gaussian graphical mixture model.  arXiv preprint arXiv:1308.3381, 2013.
  * [19] James MacQueen et al.  Some methods for classification and analysis of multivariate observations.  In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1, pages 281–297. Oakland, CA, USA, 1967\.
  * [20] Ronald C Neath et al.  On convergence properties of the monte carlo em algorithm.  In Advances in modern statistical theory and applications: A festschrift in honor of Morris L. Eaton, pages 43–62. Institute of Mathematical Statistics, 2013.
  * [21] Andrew Y Ng, Michael I Jordan, and Yair Weiss.  On spectral clustering: Analysis and an algorithm.  In Advances in neural information processing systems, pages 849–856, 2002.
  * [22] Christopher A Sims.  Macroeconomics and reality.  Econometrica: journal of the Econometric Society, pages 1–48, 1980\.
  * [23] Federico Tomasi, Veronica Tozzo, and Annalisa Barla.  Temporal patterns detection in time-varying graphical models.  In Submitted for review.
  * [24] Federico Tomasi, Veronica Tozzo, Saverio Salzo, and Alessandro Verri.  Latent variable time-varying network inference.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 2338–2346. ACM, 2018.
  * [25] Jean-Philippe Vert, Koji Tsuda, and Bernhard Schölkopf.  A primer on kernel methods.  Kernel methods in computational biology, 47:35–70, 2004.