# Approximate learning of parsimonious
Bayesian context trees

Daniyar Ghani  Department of Mathematics, Imperial College London
180 Queen’s Gate, London SW7 2AZ, United Kingdom
<EMAIL_ADDRESS>0000-0001-8611-9966 Nicholas A. Heard
Department of Mathematics, Imperial College London
180 Queen’s Gate, London SW7 2AZ, United Kingdom
<EMAIL_ADDRESS>0000-0002-8767-0810 Francesco Sanna Passino
Department of Mathematics, Imperial College London
180 Queen’s Gate, London SW7 2AZ, United Kingdom
<EMAIL_ADDRESS>0000-0002-4571-6681

###### Abstract

Models for categorical sequences typically assume exchangeable or first-order
dependent sequence elements. These are common assumptions, for example, in
models of computer malware traces and protein sequences. Although such
simplifying assumptions lead to computational tractability, these models fail
to capture long-range, complex dependence structures that may be harnessed for
greater predictive power. To this end, a Bayesian modelling framework is
proposed to parsimoniously capture rich dependence structures in categorical
sequences, with memory efficiency suitable for real-time processing of data
streams. Parsimonious Bayesian context trees are introduced as a form of
variable-order Markov model with conjugate prior distributions. The novel
framework requires fewer parameters than fixed-order Markov models by dropping
redundant dependencies and clustering sequential contexts. Approximate
inference on the context tree structure is performed via a computationally
efficient model-based agglomerative clustering procedure. The proposed
framework is tested on synthetic and real-world data examples, and it
outperforms existing sequence models when fitted to real protein sequences and
honeypot computer terminal sessions.

Keywords — categorical sequences, context trees, Markov models, model-based
clustering.

##  1 Introduction

Sequences of categorical data are ubiquitous in application areas such as
language processing, bioinformatics and cyber-security. Statistical models for
categorical sequences typically involve Markov assumptions (see, for example,
Lewis, , 2001; Mächler and Bühlmann, , 2004), which yield easily interpretable
and powerful methods for sequential prediction, changepoint detection and
classification. Fixed order-$D$ Markov models are the natural extension of
order-1 Markov chains, where the prediction of a new sequence element depends
on the values of the previous $D$ elements, known as the context. However,
while fixed Markov models can theoretically learn rich dependence structures,
they are limited by the exponential growth of parameter space. Given a state
space, or vocabulary, of $V$ elements, a fixed order-$D$ Markov model requires
$V^{D}(V-1)$ parameters. The computational cost and storage requirements
become excessively large as vocabulary size $V$ or Markov order $D$ grow. On
the contrary, models using short contexts are computationally tractable but
can have inferior predictive performance.

To retain the ability to learn complex, long-range dependence structures while
reducing the size of the parameter space, many parsimonious alternatives to
fixed higher-order Markov models have been proposed in the literature.
Rissanen,  (1983) introduce variable-order Markov models (VOMM) for the
purpose of data compression, where the Markov order depends on the context of
the current element. Unlike fixed order-$D$ Markov models, the context lengths
used for prediction in VOMMs can vary depending on the data. Inference
algorithms for these models are explored in Begleiter et al.,  (2004).

Bourguignon and Robelin,  (2004) extend variable-order Markov models in the
form of a parsimonious context tree (PCT) structure, where nodes are fused
together so that leaves may correspond to multiple contexts sharing the same
predictive distributions. Developments in PCT learning algorithms are centred
around dynamic programming (Eggeling et al., , 2013, 2019), and while the
efficiency of PCT learning has improved, the combinatorially large context
trees required at initialisation limit the feasible vocabulary to $V\approx
10$.

Sparse Markov models (Jääskinen et al., , 2014; Xiong et al., , 2016; Bennett
et al., , 2023) and minimal Markov models (García and González-López, , 2011,
2017) are further generalisations of VOMMs with no hierarchical context tree
structure; contexts of varying lengths are grouped into equivalence classes
with the same predictive distributions. Such flexible modelling can uncover
interesting dependence structures, yet inference remains challenging.

Bayesian context trees (BCT, Kontoyiannis et al., , 2022; Papageorgiou and
Kontoyiannis, , 2024) are a recent area of research focusing on exact Bayesian
inference for VOMMs. BCT inference algorithms can efficiently recover long-
range dependencies, and enable computation of quantities such as the prior
predictive likelihood and model posterior probabilities, useful for downstream
analysis. However, BCTs are still limited to small vocabularies.

In this article, parsimonious Bayesian context trees (PBCT) are introduced as
a new class of variable-order Markov model. The main differentiator of the
proposed model with existing Bayesian context trees is the clustering of
contexts, as the vocabulary is repeatedly partitioned in the tree to yield a
variable-order Markov model over clusters. The advantages of context
clustering are two-fold: a significant reduction in dimensionality and the
ability to borrow statistical strength across similar contexts. To generate
the structure of a PBCT, a natural recursive partitioning procedure is
developed using the Chinese restaurant process (CRP, Aldous, , 1985).
Dirichlet priors are placed on the unknown predictive distributions associated
with each leaf of the tree to allow simple evaluation of the marginal
likelihood of a sequence $\bm{x}$ given a context tree structure.

A novel algorithm is developed for efficient approximate inference of
parsimonious Bayesian context trees, using model-based agglomerative
clustering. The method proposed for inference enables the analysis of data
with larger vocabulary sizes compared to existing inference schemes for
similar context tree models. Results are obtained using synthetic sequences as
well as real-world examples in cyber-security and bioinformatics,
demonstrating the ability of the PBCT model to perfectly recover simulated
model structures and outperform existing models for sequential prediction.

In Section 2, the parsimonious Bayesian context tree model is described,
including the assumed generative process. Section 3 outlines the model-based
clustering inference scheme for learning PBCTs. Results of a simulation study
are presented in Section 4, and real data examples are investigated in Section
5.

##  2 Parsimonious Bayesian context trees

###  2.1 Markov models

Consider a discrete vocabulary $\mathcal{V}=\\{1,2,\dots,V\\}$, $V\geq 1$, and
a potentially infinite sequence $x_{1},x_{2},\dots$, where each
$x_{i}\in\mathcal{V}$. Suppose the first $N$ sequence elements have been
observed as $\bm{x}=(x_{1},\dots,x_{N})\in\mathcal{V}^{N}$. For $D\geq 0$, an
order-$D$ Markov model predicts the next element of the sequence $x_{N+1}$
conditional on the previous $D$ elements:

$p(x_{N+1}\,|\,x_{N},\dots,x_{1})=\begin{cases}p(x_{N+1}),&\text{if }D=0,\\\
p(x_{N+1}\,|\,x_{N},\dots,x_{N-D+1}),&\text{if }D\geq 1.\end{cases}$ (2.1)

Equivalently, the next element depends on its length-$D$ context
$(x_{N},\dots,x_{N-D+1})\in\mathcal{V}^{D}$.

Markov models can be represented naturally using context trees, whose formal
definition is given in Section 2.2. The root node of the tree represents the
unobserved next element $x_{N+1}$ of a sequence, and the non-root nodes
represent possible values of observed elements. The depth of a node
corresponds to context length. In context tree representations of standard
Markov models, nodes contain single elements of the vocabulary $\mathcal{V}$.
In this article, parsimonious context trees are introduced as representations
of parsimonious Markov models (Bourguignon and Robelin, , 2004), where nodes
contain subsets of elements of $\mathcal{V}$. Each path from the root to a
leaf node represents a collection of contexts which share the same predictive
distribution for the next sequential element.

###  2.2 Context trees

A tree $\bm{\mathcal{T}}=(\bm{\mathcal{C}},\,\bm{\mathcal{E}})$ is an
undirected graph with no cycles, where $\bm{\mathcal{C}}$ denotes the nodes of
the graph and $\bm{\mathcal{E}}\subset\bm{\mathcal{C}}\times\bm{\mathcal{C}}$
is a set of edges. In predictive models for sequences over a vocabulary
$\mathcal{V}$, the nodes $\bm{\mathcal{C}}$ of a parsimonious context tree
$\bm{\mathcal{T}}$ each contain subsets of $\mathcal{V}$, and are indexed
according to their position within the tree; specifically, each node is
associated with an index in $\mathbb{N}^{d}$ for some $0\leq d\leq D$, whose
length corresponds to its depth in the tree hierarchy, and $D\geq 0$ is the
maximum depth of the tree.

Let $\mathbb{E}_{D}=\bigcup_{d=0}^{D}\mathbb{N}^{d}$ denote the set of all
possible indices for nodes in $\bm{\mathcal{T}}$. Each node in
$\bm{\mathcal{C}}$ will be uniquely referred to as $C_{\bm{e}}$, for an index
$\bm{e}\in\mathbb{E}_{D}$. The root node $C\equiv C_{\emptyset}=\mathcal{V}$
is indexed by the empty tuple. For any node pair
$(C_{\bm{e}},\,C_{\bm{e}^{\prime}})\in\bm{\mathcal{E}}$, where
$\bm{e},\bm{e}^{\prime}\in\mathbb{E}_{D}$, the node $C_{\bm{e}^{\prime}}$ is
said to be a child of $C_{\bm{e}}$, and $\bm{e}^{\prime}$ must be the
concatenation $(\bm{e},k)$ for some $k\in\mathbb{N}$. Therefore, if
$C_{\bm{e}}$ has $K\geq 1$ children, these are indexed by sequences in
$\\{\bm{e}\\}\times\\{1,\dots,K\\}$, written
$\mathcal{C}_{\bm{e}}=\text{children}(C_{\bm{e}})=\\{C_{\bm{e}1},C_{\bm{e}2},\dots,C_{\bm{e}K}\\}$.
Crucially, every set $\mathcal{C}_{\bm{e}}$ of child nodes must form a
partition of $\mathcal{V}$, so that $\bigcup_{k=1}^{K}C_{\bm{e}k}=\mathcal{V}$
and $C_{\bm{e}k}\cap C_{\bm{e}k^{\prime}}=\emptyset$ for any $k\neq
k^{\prime}$. A leaf node has no children, and each leaf node is associated
with a unique conditional probability distribution (2.1) for the next element
in a sequence, given a context matching the corresponding path through the
tree.

A CBA BCAB CABACCBACBACBFull order-2 treeParsimonious context tree Figure 1:
Comparison of tree structures for fixed order-2 Markov and parsimonious
context tree models over the vocabulary
$\mathcal{V}=\\{\texttt{A},\texttt{B},\texttt{C}\\}$.

Figure 1 compares the tree structures representing a fixed order-2 Markov
model ($D=2$) and a parsimonious context tree with $D=3$. To understand how
the tree structures determine sequential prediction, consider a sequence
$\bm{x}\in\mathcal{V}^{N}$ over a 3-term vocabulary
$\mathcal{V}=\\{\texttt{A},\texttt{B},\texttt{C}\\}$. Suppose the last three
elements of the sequence $\bm{x}$ are observed as
$(x_{N},x_{N-1},x_{N-2})=(\texttt{A},\texttt{C},\texttt{B})$. In the order-2
model, the distribution of the next element is given by
$p(x_{N+1}\,|\,\bm{x})=p(x_{N+1}\,|\,x_{N}=\texttt{A},\,x_{N-1}=\texttt{C})$,
corresponding to the path: root $\to$ A $\to$ C. Instead, in the parsimonious
model in Figure 1, the predictive distribution corresponds to the path: root
$\to$ {A, C} $\to$ {C} $\to$ {B, C}, and is given by
$p(x_{N+1}\,|\,x_{N}\in\\{\texttt{A},\texttt{C}\\},\,x_{N-1}\in\\{\texttt{C}\\},\,x_{N-2}\in\\{\texttt{B},\texttt{C}\\})$.
The comparison illustrates the desired dimensionality reduction from fixed to
parsimonious Markov models: a fixed order-3 model implies a tree with 27
leaves, each associated with a different predictive distribution, whereas the
parsimonious model in Figure 1 defines only 4 predictive distributions.

###  2.3 Model specification

####  2.3.1 Conditional distributions

Given the parsimonious context tree $\bm{\mathcal{T}}$ with nodes
$\bm{\mathcal{C}}$, let $E\subseteq\mathbb{E}_{D}$ denote the set of leaf node
indices. Assign to each leaf node $C_{\bm{e}}$, $\bm{e}\in E$, a probability
mass function $\bm{\phi}_{\bm{e}}$ over the vocabulary $\mathcal{V}$ with
Dirichlet prior distributions:

$\bm{\phi}_{\bm{e}}\sim\text{Dirichlet}(\bm{\eta}),$ (2.2)

where $\bm{\eta}\in\mathbb{R}^{V}_{+}$. The tree $\bm{\mathcal{T}}$ and the
implied conditional distributions $\\{\bm{\phi}_{\bm{e}}\\}_{\bm{e}\in E}$
define a parsimonious Bayesian context tree (PBCT) model.

For each leaf index $\bm{e}=(\epsilon_{1},\epsilon_{2},\dots,\epsilon_{d})\in
E$, define

$\mathcal{P}_{\bm{e}}^{N}=\\{\bm{x}\in\mathcal{V}^{N}:x_{N}\in
C_{\epsilon_{1}},\,x_{N-1}\in
C_{\epsilon_{1}\epsilon_{2}},\dots,\,x_{N-d+1}\in C_{\bm{e}}\\}$ (2.3)

to be the set of $\bm{x}\in\mathcal{V}^{N}$ such that the length-$d$ context
$(x_{N},x_{N-1},\dots,x_{N-d+1})$ matches the path from the root of
$\bm{\mathcal{T}}$ to the leaf node indexed by $\bm{e}$. Suppose the leaf
corresponding to sequence $\bm{x}$ has index $\bm{e}^{\bm{x}}\in E$, such that
$\bm{x}\in\mathcal{P}_{\bm{e}^{\bm{x}}}^{N}$. Then the next element $x_{N+1}$
is sampled from the conditional distribution $\bm{\phi}_{\bm{e}^{\bm{x}}}$:

$x_{N+1}\,|\,\bm{x},\bm{\mathcal{T}},\\{\bm{\phi}_{\bm{e}}\\}_{\bm{e}\in
E}\sim\text{Categorical}(\bm{\phi}_{\bm{e}^{\bm{x}}}),$ (2.4)

yielding the Markov property
$p(x_{N+1}\,|\,\bm{x})=\bm{\phi}_{\bm{e}^{\bm{x}}}$ for sequential prediction.
Figure 2 gives an illustrative example of next-element prediction given a PBCT
$\bm{\mathcal{T}}$ and sequence $\bm{x}$.

$\text{children}(C_{12})\equiv\mathcal{C}_{12}=\\{C_{121},\,C_{122}\\}$ A CBA
BCAB CABACCBACBACBFull order-2 treeParsimonious context tree$C_{\emptyset}$A
C$C_{1}$B$C_{2}$A B$C_{11}$C$C_{12}$A$C_{121}$B
C$C_{122}$$\bm{\phi}_{11}$$\bm{\phi}_{121}$$\bm{\phi}_{122}$$\bm{\phi}_{2}$A$x_{1}$C$x_{2}$CAB$\cdots$ABCA$x_{N}$Sequence
$\bm{x}$Predict:
$x_{N+1}=\textnormal{{\color[rgb]{0.75,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.75,0,0}{{B}}}}$$\bm{x}\in\mathcal{P}_{122}^{N}=C_{1}\times
C_{12}\times C_{122}\times V^{N-3}$$x_{N+1}\sim\bm{\phi}_{122}$ {B, C}{C}{A,
C}ABCTree $\bm{\mathcal{T}}$ Figure 2:  Example workings of the proposed
Markov model for predicting the next element given a length-$N$ sequence
$\bm{x}$ and tree $\bm{\mathcal{T}}$ with 4 leaves. The vocabulary is
$\mathcal{V}=\\{\texttt{A},\texttt{B},\texttt{C}\\}$. Read the sequence from
right to left and traverse the tree starting from the root (in red). The
sequence corresponds to a path from root to a leaf via nodes $C_{1}$, $C_{12}$
and $C_{122}$. This path admits the predictive distribution $\bm{\phi}_{122}$
from which the next element $x_{N+1}$ is sampled.

####  2.3.2 Marginal likelihood

Consider a PBCT $\bm{\mathcal{T}}$ with leaf indices $E$. For
$\bm{x}\in\mathcal{V}^{N}$, define
$\bm{X}_{\bm{e}}=(X_{\bm{e},1},\dots,X_{\bm{e},V})$, where

$X_{\bm{e},v}=\sum_{n=1}^{N}\mathds{1}_{\mathcal{P}^{n-1}_{\bm{e}}}(\bm{x}_{n-1})\mathds{1}_{\\{v\\}}(x_{n}).$
(2.5)

The count $X_{\bm{e},v}$ is the number of times in $\bm{x}$ that
$v\in\mathcal{V}$ follows a context indexed by $\bm{e}$, and the subsequence
$\bm{x}_{n}=(x_{1},\dots,x_{n})$ denotes the first $n$ elements of $\bm{x}$.
Due to the conjugacy of (2.2) with (2.4), the marginal likelihood of $\bm{x}$
under $\bm{\mathcal{T}}$ is obtained from the Dirichlet–Categorical
distribution:

$p(\bm{x}\,|\,\bm{\mathcal{T}})=\prod_{e\in
E}\frac{B(\bm{X}_{\bm{e}}+\bm{\eta})}{B(\bm{\eta})},$ (2.6)

where $B(\bm{\alpha})=\prod_{i}\Gamma(\alpha_{i})/\Gamma(\sum_{i}\alpha_{i})$
is the multivariate beta function.

In practice, for the simulation study in Section 4 and real-data applications
in Section 5, the counts $\bm{X}_{\bm{e}}$ are calculated ignoring the first
$D$ elements of a sequence, where $D$ is the specified maximum tree depth. In
this way, the sum of counts are computed over the values $D+1\leq n\leq$,
avoiding any cases where the current subsequence $\bm{x}_{n}$ does not have
enough sequential history to map to a leaf in the context tree.

Another consideration for practical model inference is that discrete sequence
data commonly occur as collections of many sequences over the same vocabulary.
The computations of counts (2.5) and marginal likelihood (2.6) in this case
are modified as follows: given a collection of observed sequences, the counts
$X_{\bm{e},v}$ for each sequence are calculated separately and aggregated, and
the marginal likelihood is computed using the aggregated counts.

####  2.3.3 Generative process

This section describes a method for generating a parsimonious Bayesian context
tree graph. For each parent node in a PBCT, it is assumed that the child nodes
correspond to a partition of the vocabulary $\mathcal{V}$, as described in
Section 2.2. A natural approach to randomly partitioning a set of integers is
the Chinese restaurant process (CRP, Aldous, , 1985), a common representation
of the Dirichlet process (Ferguson, , 1973). Suppose customers are to be
seated at a restaurant with an infinite number of tables. Let
$z_{m}\in\mathbb{N}$ denote the table allocation of customer $m$. The first
customer sits at the first table, so $z_{1}=1$. For each successive customer
$m\geq 2$, if $K$ tables have have already been occupied, the customer sits at
an occupied table $k\in\\{1,\dots,K\\}$ or a new table $K+1$ with
probabilities:

$\displaystyle\mathbb{P}(z_{m}=k)=\frac{m_{k}}{\alpha+m},\quad\mathbb{P}(z_{m}=K+1)$
$\displaystyle=\frac{\alpha}{\alpha+m},$

where $\alpha>0$ controls the rate at which new tables are formed and $m_{k}$
denotes the number of customers seated at table $k$. By taking a finite
realisation of the CRP after $V$ customers are seated, each occupied table
will be associated with a subset of integers which together form an
exchangeable partition of $\mathcal{V}=\\{1,\dots,V\\}$. The proposed
generative process assumes a CRP prior distribution with parameter $\alpha$
over the set of partitions of $\mathcal{V}$, denoted $\text{CRP}_{V}(\alpha)$.

The following processes describe the generation of a PBCT
$\bm{\mathcal{T}}=(\bm{\mathcal{C}},\,\bm{\mathcal{E}})$. Let
$\mathcal{A}=\\{A_{1},\dots,A_{K}\\}$ be a random partition of the vocabulary
$\mathcal{V}$, where each $A_{k}\neq\emptyset$. For an infinite sequence of
sets $(C_{1},C_{2},\dots)$, suppose $(C_{1},C_{2},\dots)\sim H_{V}(\alpha)$
if:

$\displaystyle\mathcal{A}\sim\text{CRP}_{V}(\alpha),\quad
C_{k}=\begin{cases}A_{k},&\text{if }|\mathcal{A}|>1\text{ and
}k\leq|\mathcal{A}|,\\\ \emptyset,&\text{otherwise},\end{cases}\quad\text{for
}k\in\mathbb{N}.$

The process $H_{V}(\alpha)$ describes the construction of a sequence of sets
starting with the components of a partition of $\mathcal{V}$ generated by the
CRP, followed by a countably infinite number of empty sets. For $D\geq 1$,
define $G_{V}(\alpha,D)$ as the process:

$\displaystyle C\equiv C_{\emptyset}$ $\displaystyle=\mathcal{V},$ (2.7)
$\displaystyle(C_{\bm{e}1},C_{\bm{e}2},\ldots)$
$\displaystyle\sim\begin{cases}H_{V}(\alpha),&\text{if
}C_{\bm{e}}\neq\emptyset,\\\
\delta_{(\emptyset,\emptyset,\ldots)},&\text{otherwise},\end{cases}\quad\text{for
}\bm{e}\in\mathbb{E}_{D-1},$ (2.8) $\displaystyle\bm{\mathcal{C}}$
$\displaystyle=\\{C_{\bm{e}}\subseteq\mathcal{V}:C_{\bm{e}}\neq\emptyset,\,\bm{e}\in\mathbb{E}_{D}\\},$
(2.9) $\displaystyle\bm{\mathcal{E}}$
$\displaystyle=\\{(C_{\bm{e}},C_{\bm{e}k})\in\bm{\mathcal{C}}\times\bm{\mathcal{C}}:\bm{e}\in\mathbb{E}_{D-1},\,k\in\mathbb{N}\\},$

where $\delta_{(\emptyset,\emptyset,\ldots)}$ denotes a Dirac measure placing
probability one on the sequence $(\emptyset,\emptyset,\ldots)$. Then
$\bm{\mathcal{T}}\sim G_{V}(\alpha,D)$ is the generative process of the tree
$\bm{\mathcal{T}}$ of maximum depth $D$ with nodes $\bm{\mathcal{C}}$ and
edges $\bm{\mathcal{E}}$. The children of each non-leaf node of the tree form
a partition of $\mathcal{V}$.

There are two ways to control the complexity of a generated tree: (i) vary the
maximum depth of the tree, $D$, or (ii) vary the CRP rate parameter $\alpha$.
Smaller values of $\alpha$ lead to partitions of $\mathcal{V}$ with fewer
components, or clusters, constraining the complexity of the tree.
Additionally, to limit tree size, $\alpha$ may be chosen to decay with depth,
causing the expected number of generated child nodes to decrease with depth.

##  3 Inference

Section 2.3.3 described the assumed generative process of a parsimonious
Bayesian context tree. In this section, an approximate Bayesian inference
scheme is described to learn PBCT structure given data. The introduced method,
based on agglomerative clustering, is applied in Section 4 to infer synthetic
tree structures, and in Section 5 to estimate the dependence structures in
real protein sequences and command-line data.

###  3.1 Agglomerative clustering

Agglomerative hierarchical clustering (Duda and Hart, , 1973) is a general
method for partitioning a discrete set of elements using a similarity metric.
The procedure starts with each element in a singleton cluster, and
successively merges clusters according to their similarities until all
elements are grouped together in the same cluster. A nested sequence of
clusterings is created, and the optimal clustering is chosen to maximise an
objective function. Bayesian model-based agglomerative schemes (Heller and
Ghahramani, , 2005; Heard et al., , 2006), use marginal posterior
probabilities to define similarities between clusters.

Model-based agglomerative clustering is applied to infer the parsimonious
Bayesian context tree $\bm{\mathcal{T}}$ of maximum depth $D$ given a sequence
$\bm{x}$. The clustering process is initialised at the root node and repeated
recursively to obtain the optimal children of each node: each recursion
terminates if the optimal child configuration is $\\{\mathcal{V}\\}$, or the
current node is at the maximum depth $D$, at which point the current node is
set as a leaf of the tree.

Starting from the root node $C\equiv C_{\emptyset}$, let
$\bm{e}\in\mathbb{E}_{d}$ be an index at tree depth $d\leq D$, and let
$\mathcal{C}_{\bm{e}}=\\{C_{\bm{e}1},\ldots,\,C_{\bm{e}K}\\}$ denote the
children of node $C_{\bm{e}}$ which form a partition of the vocabulary
$\mathcal{V}$. Denote by $\bm{\mathcal{T}}_{d}$ the tree cut at the current
depth $d$. The cluster similarity at each step is the multiplicative increase
in probability obtained by merging $C_{\bm{e}i}$ and $C_{\bm{e}j}$:

$s_{i,j}=\frac{p(\bm{x}\,|\,\bm{\mathcal{T}}_{d},\,\mathcal{C}_{\bm{e}}^{(i,j)})\,p(\mathcal{C}_{\bm{e}}^{(i,j)})}{p(\bm{x}\,|\,\bm{\mathcal{T}}_{d},\,\mathcal{C}_{\bm{e}})\,p(\mathcal{C}_{\bm{e}})},$
(3.1)

where $\mathcal{C}_{\bm{e}}^{(i,j)}$ denotes the cluster configuration in
which clusters $C_{\bm{e}i}$ and $C_{\bm{e}j}$ are merged, and
$p(\bm{x}\,|\,\bm{\mathcal{T}}_{d},\,\mathcal{C}_{\bm{e}})$ is the marginal
likelihood of $\bm{x}$ under the current tree $\bm{\mathcal{T}}_{d}$ with the
proposed cluster configuration $\mathcal{C}_{\bm{e}}$ at depth $d+1$. Assuming
the generative process in Section 2.3.3, the prior probability
$p(\mathcal{C}_{\bm{e}})$ is available from the Chinese restaurant process.
The agglomerative clustering procedure for the children of each node
$C_{\bm{e}}$ at depth $d$ starts with each element of $\mathcal{V}$ belonging
to a singleton cluster, and is sequentially iterated until all elements are
merged in the same group; the clustering structure $\mathcal{C}_{\bm{e}}$
which locally maximises the marginal posterior probability
$\pi(\mathcal{C}_{\bm{e}})=p(\bm{x}\,|\,\bm{\mathcal{T}}_{d},\,\mathcal{C}_{\bm{e}})\,p(\mathcal{C}_{\bm{e}})$
is chosen as the optimal child configuration for the corresponding node in the
tree. This scheme is denoted _recursive agglomerative clustering_ (RAC), and
it is outlined in Algorithm 1.

Algorithm 1 Recursive agglomerative clustering (RAC)

1: Initialise: Set current node to root node: $C_{\bm{e}}\equiv
C_{\emptyset}=\mathcal{V}$.

2: while current node $C_{\bm{e}}$ is not a leaf do

3: Initialise clustering
$\mathcal{C}_{\bm{e}}=\\{C_{\bm{e}1},\ldots,C_{\bm{e}V}\\},$ where each
$C_{\bm{e}i}=\\{i\\}$.

4: Set number of clusters $K=V$ and calculate
$\pi_{K}=\pi(\mathcal{C}_{\bm{e}})$.

5: while $K>1$ do

6: Merge clusters $C_{\bm{e}i}$ and $C_{\bm{e}j}$ with maximal similarity
$s_{i,j}$.

7: Set $K=K-1$ and relabel clusters $i,j$ accordingly.

8: Calculate $\pi_{K}=\pi(\mathcal{C}_{\bm{e}})$ for new clustering.

9: Select the optimal clustering $\mathcal{C}_{\bm{e}}$ which maximises
$\pi(\mathcal{C}_{\bm{e}})$.

10: if $\mathcal{C}_{\bm{e}}=\\{\mathcal{V}\\}$ then

11: Set current node $C_{\bm{e}}$ as a leaf.

12: else

13: for each child node $C_{\bm{e}^{\prime}}\in\mathcal{C}_{\bm{e}}$ do

14: Set current node to child node: $C_{\bm{e}}=C_{\bm{e}^{\prime}}$.

15: if $C_{\bm{e}}$ is at maximum depth then

16: Set current node $C_{\bm{e}}$ as a leaf.

17: else

18: Repeat the procedure from line 2 to get optimal children of current node
$C_{\bm{e}}$.

The RAC algorithm provides a fast and scalable procedure for estimating the
tree which generated an observed sequence. In Section 4, the RAC algorithm is
shown to recover the true underlying context trees from simulated data in a
variety of settings. An advantage of the RAC algorithm over existing methods
is the ability to handle significantly larger vocabulary sizes than existing
methods for context tree learning. For instance, the algorithms presented in
Eggeling et al.,  (2019) require initialisation with extended trees containing
all possible combinations of child node configurations up to the specified
maximum tree depth. This limits feasible vocabulary size to $V\approx 10$,
whereas RAC can learn PBCTs in reasonable computational time for data with
$V\approx 100$, as shown in Section 5.

###  3.2 Model evaluation

Suppose the PBCT $\bm{\mathcal{T}}$ is inferred from the training sequence
$\bm{x}^{\text{train}}$. If $\bm{\mathcal{T}}$ has leaf indices $E$, the
predictive marginal likelihood of an unseen length-$N$ test sequence
$\bm{x}^{\text{test}}\in\mathcal{V}^{N}$ under the trained model is

$\displaystyle
p(\bm{x}^{\text{test}}\,|\,\bm{x}^{\text{train}},\,\bm{\mathcal{T}}\,)=\prod_{\bm{e}\in
E}\frac{B(\bm{X}_{\bm{e}}^{\text{test}}+\bm{X}_{\bm{e}}^{\text{train}}+\bm{\eta})}{B(\bm{X}_{\bm{e}}^{\text{train}}+\bm{\eta})}$
(3.2)

where $\bm{X}_{\bm{e}}^{\cdot}$ are the counts (2.5) calculated separately
over the training and test sequences. To evaluate the predictive performance
of a trained model, define the marginal log-loss $\ell(\bm{x}^{\text{test}})$
to be the mean negative predictive log-likelihood of the test sequence,

$\displaystyle\ell(\bm{x}^{\text{test}}\,|\,\bm{x}^{\text{train}},\,\bm{\mathcal{T}}\,)=-\frac{1}{N}\log
p(\bm{x}^{\text{test}}\,|\,\bm{x}^{\text{train}},\,\bm{\mathcal{T}}\,).$ (3.3)

Lower marginal log-loss indicates superior predictive ability.

To obtain the marginal log-loss (3.3) for a given tree, the predictive
distributions for each leaf are marginalised. For simulation studies, it is
also possible to define a true log-loss for a model for which the true
conditional distributions, denoted $\\{\bm{\phi}_{\bm{e}}^{*}\\}$, are assumed
known:

$\displaystyle\ell(\bm{x}^{\text{test}}\,|\,\bm{x}^{\text{train}},\,\bm{\mathcal{T}},\,\\{\bm{\phi}_{\bm{e}}^{*}\\}_{\bm{e}\in
E})=-\frac{1}{N}\sum_{i=1}^{N}\log\phi^{*}_{\bm{e}^{\bm{x}_{i}},\,x_{i}}.$
(3.4)

Estimates of the predictive distributions can be useful for efficient
downstream analysis: for each leaf $\bm{e}\in E$, denote by
$\hat{\bm{\phi}}_{\bm{e}}=(\hat{\phi}_{\bm{e},1},\dots,\hat{\phi}_{\bm{e},V})$
the posterior mean estimates of the corresponding Dirichlet–Categorical
distribution,

$\displaystyle\hat{\phi}_{\bm{e},v}=\frac{X_{\bm{e},v}^{\text{train}}+\eta_{v}}{\sum_{u=1}^{V}X_{\bm{e},u}^{\text{train}}+\eta_{u}},\quad
v\in\mathcal{V}.$ (3.5)

A similarity metric based on the adjusted Rand index (ARI, Hubert and Arabie,
, 1985) is defined to compare tree structures; for example, to compare a
fitted tree with a known, simulated tree. For each clustering at depth $d\leq
D$ in the fitted tree, the metric obtains the ARI of the most similar
clustering at depth $d$ in the simulated tree; these maximal ARIs are then
weighted by the frequencies of contexts in the training sequence and averaged.
In this way, more importance is placed on the parts of the trees which are
commonly represented in the data, and errors are penalised less for rare
contexts.

The tree similarity metric is defined as follows. Let $\bm{\mathcal{T}}_{1}$
and $\bm{\mathcal{T}}_{2}$ be trees of maximum depths $D_{1}$ and $D_{2}$,
such that $D_{1}\geq D_{2}$. For $m\in\\{1,2\\}$, let
$I_{m}\subseteq\mathbb{E}_{D_{m}}$ denote the indices of all nodes in
$\bm{\mathcal{T}}_{m}$ and $E_{m}\subseteq I_{m}$ denote the indices of leaf
nodes in $\bm{\mathcal{T}}_{m}$. For any node $\bm{e}_{m}\in I_{m}$ at depth
$d_{m}$, let $\mathcal{C}_{\bm{e}_{m}}$ denote the partition of $\mathcal{V}$
comprising the children of node $C_{\bm{e}_{m}}$. If $\bm{e}_{m}\in E_{m}$ is
a leaf, then let $\mathcal{C}_{\bm{e}_{m}}=\\{\mathcal{V}\\}$. Define the
weight $w_{\bm{e}_{1}}$ as the proportion of contexts in $\bm{x}$ associated
with node $C_{\bm{e}_{1}}$ over all nodes at the same depth in tree
$\bm{\mathcal{T}}_{1}$, and let $\text{ARI}_{\bm{e}_{1},\bm{e}_{2}}$ denote
the adjusted Rand index between clusterings $\mathcal{C}_{\bm{e}_{1}}$ in
$\bm{\mathcal{T}}_{1}$ and $\mathcal{C}_{\bm{e}_{2}}$ in
$\bm{\mathcal{T}}_{2}$. Then, for $0\leq d\leq D_{1}-1$, define the tree
similarity at depth $d+1$ as

$\sum_{\bm{e}_{1}\in\mathbb{N}^{d}\,\cap\,I_{1}}w_{\bm{e}_{1}}\max_{\bm{e}_{2}\in\mathbb{N}^{d}\,\cap\,I_{2}}\left\\{\text{ARI}_{\bm{e}_{1},\bm{e}_{2}}\right\\}.$
(3.6)

##  4 Simulation study

A simulation study enables assessment of the performance of the inference
procedure for learning a parsimonious Bayesian context tree (PBCT). This
section explores two experiments: (i) a hyperparameter study, and (ii) a model
comparison as training lengths vary.

###  4.1 Setup

In both experiments, PBCT models are randomly generated over the vocabulary
$\mathcal{V}$ as described in Section 2.3.3, and sequences are simulated using
these trees. Each simulated sequence is split into a training sequence
$\bm{x}^{\text{train}}$ and test sequence $\bm{x}^{\text{test}}$. Given a
simulated sequence, the task is to recover the true generated tree. Trees are
learned from the training sequence via the recursive agglomerative clustering
(RAC) procedure described in Section 3.1, Algorithm 1. The evaluation metrics
described in Section 3.2 – marginal log-loss (3.3), true log-loss (3.4) and
tree similarity (3.6) – are used to evaluate model performance.

The first experiment investigates the sensitivity of RAC to the simulation
parameter $\bm{\eta}\in\mathbb{R}^{V}_{+}$. Additionally, the simulated
predictive distributions are modified to capture the influence of a background
level of stochastic variation in a sequence. Let $\lambda\in[0,1]$ be a
parameter to control the strength of the stochastic variation. Then, the
conditional distribution $\bm{\phi}_{\bm{e}}^{\text{true}}$ for each leaf
$\bm{e}$ of a simulated tree is considered to be a mixture of the
$\text{Dirichlet}(\bm{\eta})$ distribution and the uniform distribution over
$\mathcal{V}$: the simulated distribution is given by
$\bm{\phi}_{\bm{e}}^{\text{true}}=(1-\lambda)\bm{\phi}_{\bm{e}}+\lambda\mathbf{1}_{V}/V,\,\,\bm{\phi}_{\bm{e}}\sim\text{Dirichlet}(\bm{\eta})$,
where $\mathbf{1}_{V}$ is a length-$V$ vector of ones. In this “spike-and-
slab” representation of $\bm{\phi}_{\bm{e}}^{\text{true}}$, $\lambda=0$ yields
the standard Dirichlet distribution, and $\lambda=1$ yields the uniform
distribution.

Synthetic parsimonious Bayesian context trees are generated over a vocabulary
of size $V=10$, and each partition of the vocabulary is drawn from
$\text{CRP}_{V}(\alpha)$ with rate parameter $\alpha=1$. The maximum depth of
the PBCTs is specified as $D=3$. Training sequences of length $N=10,000$ and
test sequences of length $N=1,000$ are simulated from the generated trees. The
Dirichlet hyperparameters are specified as $\bm{\eta}=\eta\mathbf{1}_{V}$, and
the simulation parameters $\eta$ and $\lambda$ are varied. Results are
obtained for 15 different model simulations. Model performance is evaluated by
computing the ARI-based similarity metric (3.6) over the training sequences to
compare tree structures, and marginal log-losses (3.3) are computed over the
test sequences to evaluate predictive performance.

In the second experiment, fixed and variable-order Markov models are
considered for comparison with the PBCT. Let $D$ be the maximum tree depth, so
PBCT-$D$ denotes a parsimonious Bayesian context tree of maximum depth $D$. A
fixed order-$d$ Bayesian Markov model (FBM-$d$), $d\in\\{1,\dots,D\\}$, is
defined such that Dirichlet priors are placed on the conditional distributions
for each order-$d$ context. Similarly, a variable-order Bayesian Markov model
with maximum order $D$ (VBM-$D$; see, for example, Dimitrakakis, , 2010) is
defined with Dirichlet conditional distributions for each context. VBM
structure is inferred using a modified version of RAC, Algorithm 1: at each
vocabulary partitioning step, only two outcomes are compared: (i) all elements
in the same cluster, or (ii) $V$ singleton clusters.

Predictive performance is investigated as the length of training sequence
varies. The hyperparameter pair $\bm{\eta}=\mathbf{1}_{V}$ and $\lambda=0$ are
used for sequential prediction so that the predictive distributions associated
with each leaf are standard Dirichlet distributions. The synthetic models are
generated with vocabulary size $V=10$, CRP parameter $\alpha=1$ and maximum
depth $D=3$. Predictive log-losses are averaged over 15 model simulations.

Figure 3: Boxplots to compare structures of fitted and simulated trees.
Synthetic trees are generated with vocabulary size $V=10$, CRP parameter
$\alpha=1$ and maximum depth $D=3$. $N=10,000$ training samples are simulated
using the trees for varying distributional hyperparameters $\eta$ and
$\lambda$. Boxplots of weighted ARI for each depth $d=1,2,3$ over 15 different
model simulations.

###  4.2 Results

Figure 3 illustrates boxplots of tree similarities for different simulation
parameter configurations, and Table 1 contains the average log-losses of the
simulated and fitted trees for the first experiment. It is found that the RAC
algorithm for learning PBCTs generally performs well, and can perfectly
recover synthetic context tree structures for several parameter
configurations. The best and most consistent performance is achieved by the
hyperparameter pair $\eta=1,\,\lambda=0$. It can be seen that the presence of
random noise can lead to improved performance for small values of $\eta$.
Small $\eta$ values lead to highly sparse Dirichlet distributions, leading to
issues with model identifiability. Adding noise smooths the conditional
distributions to ensure greater variation of contexts in the simulated
sequence, which helps RAC recover the simulated tree. On the other hand, as
$\eta$ becomes large, the predictive distributions each become close to
uniform, again leading to identifiability issues.

Table 1: Average log-losses for simulated and fitted PBCT models for varying distributional hyperparameter pairs $\eta$ and $\lambda$ repeated over 15 models. Standard deviations in parentheses. Synthetic trees generated with vocabulary size $V=10$, CRP parameter $\alpha=1$ and maximum depth $D=3$. Sequences of length $N=10,000$ are simulated from each tree. Best performance is given by the hyperparameters $\eta=1$, $\lambda=0$. $\eta$ | $\lambda$ | Simulated | Fitted | Difference
---|---|---|---|---
1.0 | 0.0 | 1.9297 (0.0674) | 1.9407 (0.0667) | 0.0017 (0.0023)
1.0 | 0.25 | 2.1158 (0.0501) | 2.1269 (0.0538) | 0.0093 (0.0134)
1.0 | 0.5 | 2.2177 (0.0168) | 2.2290 (0.0211) | 0.0090 (0.0099)
0.1 | 0.0 | 0.9045 (0.1314) | 1.0249 (0.2631) | 0.1236 (0.2117)
0.1 | 0.25 | 1.5552 (0.1408) | 1.5921 (0.1988) | 0.0477 (0.1123)
0.1 | 0.5 | 1.9397 (0.0337) | 1.9436 (0.0372) | 0.0042 (0.0054)
0.5 | 0.0 | 1.6800 (0.0888) | 1.6979 (0.0993) | 0.0214 (0.0252)
0.5 | 0.25 | 1.9587 (0.0693) | 1.9656 (0.0669) | 0.0088 (0.0076)
0.5 | 0.5 | 2.1503 (0.0336) | 2.1681 (0.0456) | 0.0231 (0.0309)
5.0 | 0.0 | 2.2165 (0.0262) | 2.2271 (0.0257) | 0.0110 (0.0109)
5.0 | 0.25 | 2.2581 (0.0160) | 2.2716 (0.0161) | 0.0089 (0.0099)
5.0 | 0.5 | 2.2871 (0.0054) | 2.2959 (0.0073) | 0.0094 (0.0067)

Figure 4 displays average log-losses for a selection of fitted models, as well
as the difference in log-losses with the simulated models, as training length
varies. Figures 4(a) and 4(c) show the log-losses for training lengths up to
$N=5,000$, where the PBCT outperforms fixed and variable-order Markov models
for all training lengths. Additionally, Figures 4(b) and 4(d) illustrate
continuations of the same experimental simulations for longer training
lengths, plotted on a log-scale, from $N=5,000$ up to $N=200,000$. The PBCT
model consistently recovers the simulated trees given longer training
sequences. The fixed and variable-order Markov models, FBM-3 and VBM-3,
continue to improve with longer training lengths but do not reach the same
performance as PBCT-3.

(a)

(b)

(c)

(d)

Figure 4: Average fitted log-losses and difference in log-losses between
fitted and simulated models averaged over 15 simulations, for varying training
lengths and a fixed test length of 1000 elements. Synthetic trees are
generated with vocabulary size $V=10$, CRP parameter $\alpha=1$ and maximum
depth $D=3$. (4(a)) and (4(c)) show PBCT outperforms the other Markov models
over all training lengths. (4(b)) and (4(d)) show consistent recovery of
simulated PBCT structure for longer training lengths.

##  5 Application to real data

The parsimonious Bayesian context tree model is applied to two real-world
examples of categorical sequence data: (i) Imperial honeypot terminal sessions
and (ii) UniProt protein sequences. Comparisons are made with fixed and
variable-order Markov models.

Table 2: Sequence statistics for real-world datasets. Dataset | Mean length | Total training length | Total test length
---|---|---|---
Honeypot | 52 | 46225 | 5402
UniProt | 69 | 31167 | 3370
Table 3: Predictive performance of each model evaluated on the honeypot and UniProt datasets. Best log-loss performance is achieved by PBCT for both datasets. Model | Honeypot | UniProt
---|---|---
$L$ |  log-loss | $L$ |  log-loss
FBM-0 | 1 |  2.72327 | 1 |  2.83083
FBM-1 | 93 |  1.04483 | 21 |  2.57862
FBM-2 | 8,649 |  0.72451 | 441 |  1.89768
FBM-3 | 804,357 |  - | 9,261 |  1.48856
FBM-4 | - |  - | 194,481 |  1.52664
VBM | 93 |  1.04483 | 5,641 |  1.48510
BCT111BCTs were fitted to concatenations of sequences; the software did not accept multiple sequence inputs. | 1,312 |  0.70253 | 3,820 |  1.48539
PBCT | 654 |  0.69012 | 1,076 |  1.47870

###  5.1 Honeypot terminal sessions

In this example, the data are sessions of command-line sequences collected on
a honeypot at Imperial College London. A honeypot is a type of computer
network host designed to observe the behaviour of cyber attackers when granted
access to the network. The Imperial honeypot data are in the form of sessions,
where each session is formed of a sequence of Unix terminal commands executed
by an attacker. A collection of 1000 sessions, recorded between May 2021 and
January 2022, are considered for model training and evaluation. 90% of the
sessions are used for training and the remaining 10% are used to evaluate
predictive performance. Table 2 details the average session lengths and total
number of commands in the data.  From the honeypot sessions, each command is
considered a word in the vocabulary $\mathcal{V}$. After preprocessing
commands into common categories and key Unix commands, the vocabulary has
$V=93$ elements. The PBCT model is trained using CRP rate parameter $\alpha=1$
and Dirichlet hyperparameter $\bm{\eta}=\mathbf{1}_{V}$. Fixed Bayesian Markov
models of orders $D\in\\{0,1,2,3\\}$ are defined (FBM-$D$, where FBM-0 is a
bag-of-words model) and the variable-order Bayesian Markov (VBM) model is
described in Section 4. For comparison, a Bayesian context tree (BCT) is
fitted using the authors’ publicly available software (Kontoyiannis et al., ,
2022). For each model, the conditional distributions $\bm{\phi}_{\bm{e}}$ used
for prediction are Dirichlet with parameter $\bm{\eta}=\mathbf{1}_{V}$. A
maximum depth of $D=3$ is specified before training the variable-order models
PBCT, VBM and BCT.

Table 3 outlines the means and standard deviations of log-losses computed on
the test sequences for each model. The number of conditional distributions
defined by each model is denoted by $L\in\mathbb{N}$. The best predictive
performance is achieved by the PBCT. The fitted PBCT has $L=654$ and a maximum
depth of 2, exhibiting a significant reduction in dimensionality while
achieving better predictive performance than the fixed order-2 Markov model
(FBM-2), which has $L=8,649$. The fitted BCT, also of maximum depth $D=2$, has
$L=1,132$ and is slightly outperformed by the PBCT in terms of predictive log-
loss. The log-loss for FBM-3 is omitted due to excessive computational cost
($L=804,357$).

$\text{children}(C_{12})\equiv\mathcal{C}_{12}=\\{C_{121},\,C_{122}\\}$ A CBA BCAB CABACCBACBACBFull order-2 treeParsimonious context tree$C_{\emptyset}$A C$C_{1}$B$C_{2}$A B$C_{11}$C$C_{12}$A$C_{121}$B C$C_{122}$$\bm{\phi}_{11}$$\bm{\phi}_{121}$$\bm{\phi}_{122}$$\bm{\phi}_{2}$A$x_{1}$C$x_{2}$CAB$\cdots$ABCA$x_{N}$Sequence $\bm{x}$Predict: $x_{N+1}=\textnormal{{\color[rgb]{0.75,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.75,0,0}{{B}}}}$$\bm{x}\in\mathcal{P}_{122}^{N}=C_{1}\times C_{12}\times C_{122}\times V^{N-3}$$x_{N+1}\sim\bm{\phi}_{122}$ {B, C}{C}{A, C}ABCTree $\bm{\mathcal{T}}$ .S4Y | .LAYER | .PTMX
---|---|---
.HUMAN | .MISA | SH
WRITE_FILE | OTHER |
$x_{N-1}$ cd
---
$x_{N}$HOME_DIROTHER$0$$0.5$$1$$p(x_{N+1}\,|\,\bm{x})$Predict next command
$x_{N+1}$ Figure 5: Example of learned command contexts and prediction using
the PBCT fitted to honeypot sessions.

Figure 5 illustrates an example of a learned order-2 context in the PBCT and
the two largest predictive probabilities, calculated as posterior means. The
example demonstrates the ability to predict the behaviour of a network
intruder: an attacker first attempts installation of a MIRAI variant (such as
.PTMX or .MISA; see Sanna Passino et al., , 2023, for further discussion),
then changes directory using the cd command. Following this context, the PBCT
predicts the common behaviour of navigating back to the home directory.

###  5.2 Protein sequences

The UniProt knowledgebase (The UniProt Consortium, , 2016) is a large
collection of protein sequences. Each entry is a sequence with elements from a
vocabulary of the 20 standard amino acids. There is an extra element added to
the vocabulary in the case of a missing or erroneous amino acid, so the final
vocabulary has size $V=21$. A random sample of sequences of lengths between 10
and 100 is considered for computational feasibility. Table 2 details the
average and total protein sequence lengths used for analysis.

The PBCT model is fitted to the protein sequences using a CRP rate parameter
$\alpha=1$ and Dirichlet hyperparameter $\bm{\eta}=\mathbf{1}_{V}$. The fixed
and variable-order Bayesian Markov (FBM, VBM) and Bayesian context tree (BCT)
models are considered for comparison. Each variable-order model is specified a
maximum depth $D=6$ before training. Table 3 contains the results of fitting
each model to the UniProt data. The fitted PBCT, VBM and BCT models each have
maximum depth $D=4$. The PBCT achieves the best predictive performance in
terms of log-loss with fewer parameters than the comparable models FBM-3, VBM
and BCT.

$\text{children}(C_{12})\equiv\mathcal{C}_{12}=\\{C_{121},\,C_{122}\\}$ A CBA BCAB CABACCBACBACBFull order-2 treeParsimonious context tree$C_{\emptyset}$A C$C_{1}$B$C_{2}$A B$C_{11}$C$C_{12}$A$C_{121}$B C$C_{122}$$\bm{\phi}_{11}$$\bm{\phi}_{121}$$\bm{\phi}_{122}$$\bm{\phi}_{2}$A$x_{1}$C$x_{2}$CAB$\cdots$ABCA$x_{N}$Sequence $\bm{x}$Predict: $x_{N+1}=\textnormal{{\color[rgb]{0.75,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.75,0,0}{{B}}}}$$\bm{x}\in\mathcal{P}_{122}^{N}=C_{1}\times C_{12}\times C_{122}\times V^{N-3}$$x_{N+1}\sim\bm{\phi}_{122}$ {B, C}{C}{A, C}ABCTree $\bm{\mathcal{T}}$ C | L
---|---
S | R
$x_{N-3}$ $\neg\,$C
---
$x_{N-2}$ T | X
---|---
$x_{N-1}$ D
---
$x_{N}$RKN$0$$0.25$$0.5$$p(x_{N+1}\,|\,\bm{x})$Predict next amino acid
$x_{N+1}$ Figure 6: Example of protein motif discovery and prediction using
the PBCT fitted to UniProt sequences. Here, the shorthand $\neg\,$C denotes
all amino acids excluding C.

Figure 6 shows a learned order-4 context and the corresponding estimated
probabilities for predicting the next amino acid. Aside from prediction, this
example demonstrates an application of the PBCT to the discovery of protein
sequence motifs (see, for example, Bailey, , 2007). Sequence motifs are short
subsequences of amino acids used to characterise groups of proteins by
functional or structural similarity. Fitting PBCTs to collections of protein
sequences may lead to the discovery of new motifs, represented by learned
contexts.

##  6 Conclusion

A novel inference method has been introduced to efficiently capture complex
dependence structures and make predictions in discrete sequences. The proposed
parsimonious Bayesian context tree model has been verified using a simulation
study for a range of model configurations, and outperforms existing context
tree models in real-world examples while reducing the parameter space. Key
advantages of the PBCT model include scalability to large vocabularies and
interpretability of results. PBCTs can be applied in practice for prediction
tasks, anomaly detection and changepoint analysis, in addition to the protein
motif discovery application discussed in Section 5.2. An interesting extension
of the inference method would involve updates of tree structures given
streaming data.

Markov chain Monte Carlo (MCMC) provides another method for inference of PBCT
structure given data. MCMC via Gibbs sampling can be implemented by
iteratively resampling from the posterior distributions of parent-child node
structures. At each vocabulary clustering step, the results of agglomerative
clustering can be used to initialise a Gibbs sampler. In such an
implementation of MCMC, preliminary testing showed that the fitted tree
structures rarely improved on the results using only RAC, Algorithm 1.

SUPPLEMENTARY MATERIAL

The Python library pbct contains the code used to fit parsimonious Bayesian
context trees, available in the GitHub repository at
https://github.com/daniyarghani/pbct.

## References

  * Aldous, (1985) Aldous, D. J. (1985).  Exchangeability and related topics.  In Lecture Notes in Mathematics, pages 1–198. Springer Berlin Heidelberg.
  * Bailey, (2007) Bailey, T. L. (2007).  Discovering sequence motifs.  Methods in Molecular Biology, 395:271–292.
  * Begleiter et al., (2004) Begleiter, R., El-Yaniv, R., and Yona, G. (2004).  On prediction using variable order Markov models.  Journal of Artificial Intelligence Research, 22:385–421.
  * Bennett et al., (2023) Bennett, I., Martin, D. E. K., and Lahiri, S. N. (2023).  Fitting sparse Markov models through a collapsed Gibbs sampler.  Computational Statistics, 38:1977–1994.
  * Bourguignon and Robelin, (2004) Bourguignon, P. Y. and Robelin, D. (2004).  Modèles de Markov parcimonieux : sélection de modèle et estimation.  In Proceedings of the 5e édition des Journées Ouvertes en Biologie, Informatique et Mathématiques, Montréal.
  * Dimitrakakis, (2010) Dimitrakakis, C. (2010).  Bayesian variable order Markov models.  In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9, pages 161–168. PMLR.
  * Duda and Hart, (1973) Duda, R. O. and Hart, P. E. (1973).  Pattern classification and scene analysis, volume 3.  Wiley New York.
  * Eggeling et al., (2013) Eggeling, R., Gohr, A., Bourguignon, P.-Y., Wingender, E., and Grosse, I. (2013).  Inhomogeneous parsimonious Markov models.  In Machine Learning and Knowledge Discovery in Databases, pages 321–336. Springer Berlin Heidelberg.
  * Eggeling et al., (2019) Eggeling, R., Grosse, I., and Koivisto, M. (2019).  Algorithms for learning parsimonious context trees.  Machine Learning, 108(6):879–911.
  * Ferguson, (1973) Ferguson, T. S. (1973).  A Bayesian analysis of some nonparametric problems.  The Annals of Statistics, 1(2):209–230.
  * García and González-López, (2011) García, J. E. and González-López, V. A. (2011).  Minimal Markov models.  In Fourth Workshop on Information Theoretic Methods in Science and Engineering: Proceedings, pages 25–28. University of Helsinki.
  * García and González-López, (2017) García, J. E. and González-López, V. A. (2017).  Consistent estimation of partition Markov models.  Entropy, 19(4):160.
  * Heard et al., (2006) Heard, N. A., Holmes, C. C., and Stephens, D. A. (2006).  A quantitative study of gene regulation involved in the immune response of Anopheline mosquitoes.  Journal of the American Statistical Association, 101(473):18–29.
  * Heller and Ghahramani, (2005) Heller, K. A. and Ghahramani, Z. (2005).  Bayesian hierarchical clustering.  In Proceedings of the 22nd International Conference on Machine Learning, ICML ’05, pages 297–304.
  * Hubert and Arabie, (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of Classification, 2(1):193–218.
  * Jääskinen et al., (2014) Jääskinen, V., Xiong, J., Corander, J., and Koski, T. (2014).  Sparse Markov chains for sequence data.  Scandinavian Journal of Statistics, 41(3):639–655.
  * Kontoyiannis et al., (2022) Kontoyiannis, I., Mertzanis, L., Panotopoulou, A., Papageorgiou, I., and Skoularidou, M. (2022).  Bayesian context trees: Modelling and exact inference for discrete time series.  Journal of the Royal Statistical Society Series B, 84(4):1287–1323.
  * Lewis, (2001) Lewis, P. (2001).  A likelihood approach to estimating phylogeny from discrete morphological character data.  Systematic Biology, 50(6):913–925.
  * Mächler and Bühlmann, (2004) Mächler, M. and Bühlmann, P. (2004).  Variable Length Markov Chains: Methodology, Computing, and Software.  Journal of Computational and Graphical Statistics, 13(2):435–455.
  * Papageorgiou and Kontoyiannis, (2024) Papageorgiou, I. and Kontoyiannis, I. (2024).  Posterior Representations for Bayesian Context Trees: Sampling, Estimation and Convergence.  Bayesian Analysis, 19(2):501–529.
  * Rissanen, (1983) Rissanen, J. (1983).  A universal data compression system.  IEEE Transactions on Information Theory, 29(5):656–664.
  * Sanna Passino et al., (2023) Sanna Passino, F., Mantziou, A., Ghani, D., Thiede, P., Bevington, R., and Heard, N. A. (2023).  Nested Dirichlet models for unsupervised attack pattern detection in honeypot data.  arXiv e-prints, 2301.02505.
  * The UniProt Consortium, (2016) The UniProt Consortium (2016).  UniProt: the universal protein knowledgebase.  Nucleic Acids Research, 45(D1):D158–D169.
  * Xiong et al., (2016) Xiong, J., Jääskinen, V., and Corander, J. (2016).  Recursive learning for sparse Markov models.  Bayesian Analysis, 11(1):247–263.