#  Markov Modeling of Time-Series Data using Symbolic Analysis

Devesh K. Jha  MERL, Cambridge, MA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Markov models are often used to capture the temporal patterns of sequential
data for statistical learning applications. While the Hidden Markov modeling-
based learning mechanisms are well studied in literature, we analyze a
symbolic-dynamics inspired approach. Under this umbrella, Markov modeling of
time-series data consists of two major steps- discretization of continuous
attributes followed by estimating the size of temporal memory of the
discretized sequence. These two steps are critical for the accurate and
concise representation of time-series data in the discrete space.
Discretization governs the information content of the resultant discretized
sequence. On the other hand, memory estimation of the symbolic sequence helps
to extract the predictive patterns in the discretized data. Clearly, the
effectiveness of signal representation as a discrete Markov process depends on
both these steps. In this paper, we will review the different techniques for
discretization and memory estimation for discrete stochastic processes. In
particular, we will focus on the individual problems of discretization and
order estimation for discrete stochastic process. We will present some results
from literature on partitioning from dynamical systems theory and order
estimation using concepts of information theory and statistical learning. The
paper also presents some related problem formulations which will be useful for
machine learning and statistical learning application using the symbolic
framework of data analysis. We present some results of statistical analysis of
a complex thermoacoustic instability phenomenon during lean-premixed
combustion in jet-turbine engines using the proposed Markov modeling method.

##  1\. Motivation and Introduction

Markov models are used for statistical learning of sequential data when we are
interested in temporal behavior of data and there is a requirement to relax
the i.i.d. assumptions on measured data. Various different models have been
proposed and analyzed in literature. Some of the common approaches are
Autoregressive or AR models [1] and the Hidden Markov models [2]. In contrast,
in this paper we are interested in a non-linear symbolic analysis-based Markov
modeling of time-series data which is not so well studied in literature.
Symbolic Time Series Analysis (STSA) [3] is a non-linear technique for
representing temporal patterns in sequential data. STSA borrows concepts from
symbolic dynamics and chaos theory for non-linear analysis of irregular time
series data for systems which inherently do not seem stochastic. Symbolic time
series analysis consists of two critical steps, discretization where the
continuous attributes of the sequential data are projected onto a symbolic
space which is followed by identification of concise probabilistic patterns
that help compress the discretized data. Under this umbrella, finite-memory
Markov models have been shown to be a reasonable finite-memory approximation
(or representation) of systems with fading memory (e.g., engineering systems
that exhibit stable orbits or mixing) [4, 5]. Once the continuous data is
discretized, the memory estimate for the discretized sequence is used to
compress it as a finite-memory Markov process, which is represented by a state
transition matrix. This helps find the causal, dynamical structure intrinsic
to the symbolic process we are investigating, ideally to extract all the
patterns in it that may have any predictive power. The transition matrix could
be estimated by frequency counting under the assumption of infinite data and a
uniform prior for all the elements of the transition matrix. It is noted that
in this paper we only consider discrete Markov processes with finite memory
(or order).

It is desirable from a machine learning perspective to establish the
fundamental limits of information that could be learned from the underlying
data when it is compressed as a Markov chain using the symbolic analysis
framework. This problem is complex from the following perspectives.

  1. (1)

The underlying model for data is unknown which makes data discretization
difficult to analyze.

  2. (2)

The order of the discrete stochastic process depends on the discretization
parameters (i.e., number and locations of partitions) and this coupling is
difficult to analyze in the absence of any underlying model.

These two reasons make analysis of this problem very difficult from a
dynamical systems perspective. As the data model is unknown, there is no
metric to characterize the properties of the discrete sequence w.r.t. the
original data. Furthermore, the behavior of the discrete sequence is governed
by the discretization process. As a result, characterizing the final Markov
model w.r.t. the original data is extremely difficult due to this composite
process. Dynamical systems theory provides some consistent characterization of
partitions but they are in general very difficult to estimate numerically in
the absence of a model. On the other hand, there is no information-theoretic
or statistical formalism for the discretization process in open literature.
However, in the absence of a model, use of information-theoretic and
statistical measures for data discretization seems natural.

While the discretization step decides the information content of the symbolic
sequence, memory estimation is critical for concise, yet precise,
representation of the discretized sequence. A lot of techniques could be found
in open literature for discretization as well as memory estimation of Markov
processes. For machine learning applications, most of the methods tie a
technique with an end objective which can then be used to find a solution that
satisfies the end objective. However, most approaches consider the problems of
data discretization and memory estimation as separate; there is no unifying
theory for signal representation as a Markov model. Moreover, there are no
results on the interplay between the complexity of the discrete dynamical
system and the discretization process e.g., the cardinality of the discrete
set, the discretization technique, etc.. In particular, there is no Bayesian
inference technique that combines these steps together to estimate the model,
together with its parameters, for representation of the time series signal. It
is noted that the word symbolization is often interchangeably used for
discretization throughout the paper. The word order is often interchangeably
used for depth.

Symbolization is carried out via partitioning of the phase-space of the system
within which the system dynamics evolves. In general, there are two main lines
of thought behind the symbolization process, one inspired by dynamical systems
theory and the other inspired by machine learning objectives. Some approaches
inspired by the dynamical systems theory could be found in [6, 7, 8, 9].
Several partitioning methods have been proposed in literature such as maximum
entropy partition [10], symbolic aggregate approximation (SAX) [11],
maximally-bijective partition [12]. In these approaches, the discretization is
not tied to any performance objective expected from the data (e.g., in maximum
entropy partitioning, each partition contains equal number of data points and
thus it presents an unbiased discrete representation of the system). This is
somewhat different from machine learning-inspired techniques where an
optimization over the parameters of partitions is performed for near-optimal
results; however, these techniques are always susceptible to over-fitting of
data and may lead to overly-capable algorithms. Most of the machine learning
inspired discretization techniques are mainly based on some end objectives
like class separability or some unsupervised metrics like entropy minimization
of the resulting discrete data [13, 14]. Good reviews of partitioning methods
for time-series symbolization can be found in [15, 16]. A information theory-
based approach to select alphabet size was presented recently in [17]. Some
empirical results regarding the comparison and performance of different
methods have been reported in [18, 19, 20]. Even though a lot of techniques
have been reported in open literature, there is no standard approach to it.
This is for the reason that the effectiveness of a discretization technique
depends on a lot of factors like the nature of the dynamical system, the
topological properties of the phase space of the system, the similarity
metrics etc.. Moreover, in statistical learning problems, the underlying model
generating the data is unknown and thus, it is difficult to impose a
consistency criterion based on system behavior. Consequently, finding a
universal rule which can be used for a wide variety of systems with consistent
performance is evasive. To the best of author’s knowledge, there are no
consistency results reported in open literature for discretization of signals
for representation and learning.

The next step in the process is modeling of the symbol sequence which allows
further compression of the data as a Markov model. Working in the symbolic
domain, the task is to identify concise probabilistic models that relate the
past, present and the future states of the system under consideration. For
Markov modeling of the symbol sequence, this is achieved by first estimating
the depth (or size of memory) for the discrete symbolic process and then,
estimating the approximate the stochastic model from the observed sequence.
Various approaches have been reported in open literature for order estimation
of Markov chains. Most of these approaches are inspired by information theory
and makes use of concepts like the Minimum description length (MDL) [21, 22],
complexity theory [23] and/or some information criterion like Akaike’s
Information Criterion (AIC) [24] or Bayesian Information Criterion (BIC) [25,
22]. Some important results regarding consistency of these techniques have
also been reported for order estimation. In [26, 27], the authors show that
the minimum description length Markov estimator will converge almost surely to
the correct order if the alphabet size is bounded a priori. In [28, 29, 30,
31], the authors present a more direct way of order estimation by making use
of convergence rate of $k$-block distributions. Some other techniques
techniques based on the use of some heuristic information theory-based
criterion for order estimation could be found in [32, 33, 34]. For machine
learning application, the estimation process follows wrapper approach [2]
where a wrapper search algorithm with a certain stopping criteria calls the
main modeling module to build several temporal models with varying depths and
the search is stopped when the stopping criteria such as information gain or
entropy rate show marginal improvement with the added complexity [35]. This
approach in essence creates all models first and then chooses one model based
on some given metric and threshold. Making multiple passes through the data
searching for correct depth is computationally intensive and might be
infeasible for large data sets that are common today.

Recently, a new approach for depth estimation has been proposed based on the
spectral properties of the one-step transition matrix [36]. An upper bound on
the size of temporal memory has been derived that requires a single-pass
through the time-series data. In [37], the authors presented a rigorous
comparative evaluation of spectral property-based approach with three popular
techniques – (1) log-likelihood, (2) state-splitting using entropy-rate, and
(3) signal reconstruction error based approach. The latter three techniques
fall under the wrapper approach of depth estimation and therefore
computationally. They were found to be in close agreement about the estimated
depth.

Once these two parameters are estimated, the data is represented by a
stochastic matrix for the inferred Markov chain. The stochastic matrix could
be estimated by a Bayesian approach where the prior and posterior of the
individual elements of the matrix can be represented by Dirichlet and
Multinomial distributions, respectively. Under the assumption of finite state-
space of the Markov chain and an uniform prior, it can be shown that the
posterior estimates asymptotically converge to the Maximum Likelihood estimate
for the individual elements. Thus, through this sequential process of
discretization and order estimation, we render a generative Markov model for
the data. This representation for the signal can then be used for various
learning algorithms like modeling, classification, pattern matching etc..

Even though we see that a lot of work has been done on discretization and
approximation of discrete sequences as finite-order Markov chains, there is no
theory or approach that ties them together for signal representation.
Furthermore, as these pieces have not been studied together there is no result
on consistency of signal representation as a Markov model, where the
underlying model generating the data is unknown. As a result, statistical
modeling of data for STSA-based Markov modeling is still largely practiced in
an ad-hoc fashion following a wrapper-inspired technique where the solution
space is searched exhaustively for solutions and thus, still remains hugely
dependent on domain expertise. The challenges are both algorithmic and
computational. The algorithmic challenge is to synthesize a consistent
framework with guarantees on performance while the computational challenge is
to reduce computations required to arrive at a desired solution mainly
inspired by machine learning applications.

At this point we would like to point out that the presented formalism for
statistical learning is different from the standard Hidden Markov models in
the following ways.

  * •

The state-space of these models is inferred from the observed data while in
HMM, the state-space is never observed and thus it is difficult to infer.

  * •

Under the assumption that the observed sequence is a finite-order Markov
chain, the estimation of parameters is simplified. The sufficient statistics
could be obtained by estimating the order of the Markov chain and the symbol
emission probabilities conditioned on the memory words of the discrete
sequence.

Thus, the current framework of STSA is simplistic and can be easily used for
embedded applications for the simplicity of the inference algorithms. The
details are provided in the subsequent sections.

In this paper, we will present the state-of-the-art mathematical formalism of
these concepts which may be possibly used for signal representation. We
present a review and propose some mathematical problems which are required to
understand the symbolic time-series analysis framework. In order to do so, we
will review and discuss the concept of Markov partitions from dynamical
systems literature and order estimation of Markov chains from probability and
information theory. We will study some properties of Markov partitions and
present an estimation technique for order estimation of a Markov chain that
guarantees asymptotic convergence to the true order of the Markov chain. The
closure of the paper will discuss some possible directions of research. In
particular, we will discuss some possible ways in the underlying framework of
STSA can be modeled from an information-theoretic and statistical perspective.

The paper is organized in eight sections including the current section.
Section 2 briefly presents the underlying mathematical preliminaries and
prepares the background for presentation of the subsequent material. Section 3
presents the statement of the problem along with a list of major assumptions.
Section 4 presents some results on a particular type of well-behaved
partitions studied in dynamical systems literature. Section 5 presents an
estimation procedure for order of a discrete Markov chain with guarantees on
convergence of the estimates. Section 6 presents some details on estimation of
statistical parameters for the Markov models. In section 7.2, we present some
details of application to a complex thermoacoustic phenomena in combustion
using the presented ideas of Markov modeling. Finally the paper is summarized
and concluded in Section 8 along with recommendations for future research.

##  2\. Mathematical Preliminaries and Background

In this section, we explain briefly the related ideas and concepts for
symbolic time-series analysis-based Markov modeling of data. In particular, we
will introduce the concepts of discretization, order (memory) and Markov
modeling. Symbolic analysis of time-series data is a recent approach where
continuous sensor data are converted to symbol sequences via partitioning of
the continuous domain [11, 35, 38]. The discrete symbolic dynamic process is
then compressed as a Markov chain whose states are collection of words over a
finite alphabet with finite length. In this section, we provide briefly the
related concepts which are required to explain the subsequent material. We
present some relevant definitions and concepts which are used to describe the
subsequent materials.

######  Definition 2.1 (Pre-Image).

Given $\phi:M\rightarrow L$ (assume $\phi$ is surjective), the image of $x\in
M$ is $\phi(x)\in L$. Then, the pre-image of $y\in L$ is given by the set
$\phi^{-1}(y)=\\{x\mid\phi(x)=y\\}$.

######  Definition 2.2 (Homeomorphism).

A function $\phi:M\rightarrow N$ from one metric space to another is
continuous if, whenever $x_{n}\rightarrow x$ in $M$, then
$\phi(x_{n})\rightarrow\phi(x)$ in $N$. If $\phi$ is continuous, one-to-one,
onto and has a continuous inverse, then we call $\phi$ a homeomorphism.

######  Definition 2.3 (Dynamical System).

A dynamical system $(M,\phi)$ consists of a compact metric space $M$ together
with a continuous map $\phi:M\rightarrow M$. If $\phi$ is homeomorphism we
call $(M,\phi)$ to be an invertible dynamical system.

######  Definition 2.4 (Irreducible Dynamical System).

A dynamical system $(M,\phi)$ is said to be irreducible if for every pair of
open sets $U,V$ there exists $m\geq 0$ such that $\phi^{m}U\cap V\neq 0$.

######  Definition 2.5 (Bilaterally Transitive).

A point $p\in M$ is called to be bilaterally transitive if the forward orbit
$\\{\phi^{n}p\mid n>0\\}$ and the backward orbit $\\{\phi^{n}p\mid n<0\\}$ are
both dense in $M$.

A homeomorphism $\phi$ is said to be expansive if there exists a real number
$c>0$ such that $d(\phi^{n}p,\phi^{n}q)<c$ for all $n\in\mathds{Z}$, then
$p=q$.

######  Definition 2.6.

For two dynamical systems $(M,\phi)$ and $(L,\psi)$ we call the second a
factor of the first and the first an expansion of the second if there exists a
map $\pi$ of $M$ into $L$, which we call a factor map such that

  * •

$\psi\pi=\pi\phi$.

  * •

$\pi$ is continuous.

  * •

$\pi$ is onto.

Furthermore, we say that $\pi$ is a finite factor map or that it is bounded
one-to-one if

  * •

there is a bound to the number of pre-images.

  * •

every doubly transitive point has a unique pre-image.

We are interested in symbolic representation of dynamical systems $(M,\phi)$.
To understand the idea, consider the case where $\phi$ is invertible so that
all the iterates of $\phi$, positive and negative, are used. Then, to describe
the orbits $\\{\phi^{n}(y):n\in\mathds{Z}\\}$ of points $y\in M$, we can try
to use an approximate description constructed in the following way. Divide $M$
into a finite number of pieces $E_{0},E_{1},\dots,E_{r-1}$ which cover the
space $M$ and are mutually disjoint. Then, we can track the orbit of a point
$y\in M$ by keeping a record of which of these pieces $\phi^{n}(y)$ lands in.
This yields a corresponding point $x=\dots
x_{-1}x_{0}x_{1}\dots\in\\{0,1,\dots r-1\\}^{\mathds{Z}}=X_{[r]}$ (full
$r$-shift space) defined by

$\phi^{n}(y)\in E_{x_{n}}\text{ for }n\in\mathds{Z}$

Thus for every $y\in M$, we get a point $x$ in the full $r$-shift, and the
definition shows that the image $\phi(y)$ corresponds to the shifted point
$\sigma(x)$.
More formally, the partitions can be defined as follows.

######  Definition 2.7 (Partitions).

We call a family of sets $\mathcal{R}=\\{R_{1},R_{2},\dots,R_{N-1}\\}$ a
topological partition for a compact metric space $M$ if the following holds
true.

  1. (1)

each $R_{i}$ is open;

  2. (2)

$R_{i}\cap R_{j}=\emptyset$, $i\neq j$;

  3. (3)

$M=\overline{R}_{1}\cup\overline{R}_{2}\cup\dots\overline{R}_{N-1}$

As there are only a finite number of sets in the partition, every set is
denoted by a symbol from a finite alphabet and thus the partitioning process
for a dynamical system could also be viewed as identification of a many-to-one
projection mapping such that the continuous dynamical system is mapped onto a
discrete space. More formally, let the continuously varying physical process
be modeled as a finite-dimensional dynamical system,
$\dot{\mathbf{x}}=f(\mathbf{x}(t))$ where $t\in[0,\infty)$ and
$\mathbf{x}\in\mathbb{R}^{n}$ is the state-vector in the compact phase space
$\Omega$ of the system. Then, the partitioning process could be defined as a
map $\varphi:\Omega\rightarrow\mathcal{A}$ such that $\varphi(\mathbf{x})=a$
if $\mathbf{x}\in{R}_{a}$, where $R_{a}\subseteq\Omega$,
$\bigcup\limits_{a\in\mathcal{A}}R_{a}=\Omega$.

Let $\mathcal{A}$ denote the ordered set of $n$ symbols. The phase space of
the symbolic system is the space.

(1) $X_{[n]}=\mathcal{A}^{\mathds{Z}}=\\{a=(a_{k})_{k\in\mathds{Z}}\mid
a_{k}\in\mathcal{A}\\}$

of all bi-infinite sequences of elements from a set of $n$ symbols. The shift
transformation $\sigma$ is defined by shifting each bi-infinite sequence one
step to the left. This is expressed as

$(\sigma a)_{k}=a_{k+1}$

The symbolic system $(X_{[n]},\sigma)$ is called the full n-shift. Restricting
the shift transformation to a full-shift $X_{[n]}$ to a closed shift-invariant
subspace $X$, we get a very general kind of dynamical system $(X,\sigma)$
called a subshift. Given a symbolic phase space $X$, we call a $k$-tuple an
allowable $k$-block if it equals $a_{[m,m+k-1]}$ for some sequence $a$. Then,
we define a shift of finite type, also called the topological Markov shift, as
the subshift (i.e., shift-invariant) of a full shift restricted to the set
$X_{G}$ of bi-finite paths in a finite directed graph $G$ derived from a full
one by possibly removing some edges. The space could also be denoted by
$X_{V}$ where $V$ is a matrix and $v_{ij}$ denotes the number of edges going
out from node $i$ to node $j$. This gives the resemblance to a Markov chain as
normalizing each row of $V$ gives a stochastic matrix which defines a Markov
chain over $G$. Before we move to more details of theory of discrete ergodic
processes, which provides a statistical view to the discrete stochastic
processes, we would like to clarify that a symbolic sequence in a topological
Markov shift is bilaterally transitive if every admissible block appears in
both directions and infinitely often. In the subsequent section we use the
notation $BLT(A)$ to denote the subset of bilaterally transitive points in
$A\subseteq M$.
A (discrete-time, stochastic) process is a sequence
$S_{1},S_{2},\dots,S_{n},\dots$ of random variables defined on a probability
space $(S,\mathcal{E},P)$. The process has alphabet $\mathcal{A}$ if the range
of each $S_{i}$ is contained in $\mathcal{A}$. In this paper, we focus on
finite-alphabet processes, so, unless otherwise specified, “process” means a
discrete-time finite-alphabet process. Also, unless it is specified otherwise,
“measure” will mean “probability measure” and “function” will mean “measurable
function” w.r.t. some appropriate $\sigma$-algebra on a probability space. The
cardinality of a finite set $A$ is denoted by $|A|$. The sequence
$a_{m},a_{m+1},\dots,a_{n}$, where each $a_{i}\in\mathcal{A}$ is denoted by
$a_{m}^{n}$; the corresponding set of all such $a_{m}^{n}$ is denoted by
$\mathcal{A}_{m}^{n}$ except for $m=1$, when $\mathcal{A}^{n}$ is used. The
$k^{th}$ order joint distribution of the process $\\{S_{k}\\}$ is the measure
$P_{k}$ on $\mathcal{A}^{k}$ defined by the formula

(2)
$P_{k}(a_{1}^{k})=\text{Pr}(S_{1}^{k}=a_{1}^{k}),a_{1}^{k}\in\mathcal{A}^{k}$

######  Definition 2.8 (Stationary Process).

A process is stationary, if the joint distribution do not depend on the choice
of the time origin, i.e.,

(3) $\text{Pr}(S_{i}=a_{i},m\leq i\leq n)=\text{Pr}(S_{i+1}=a_{i},m\leq i\leq
n),$

for all $m,\quad n$ and $a_{m}^{n}$.

These stationary finite-alphabet processes serve as models in many settings of
interest like physics, data transmission, statistics, etc.. In this paper, we
are mainly interested in statistical and machine learning applications like
density estimation, pattern matching, data clustering, and estimation.

The simplest example of a stationary finite-alphabet process is an
independent, identically distributed distributed (i.i.d.) process. A sequence
of random variables, $\\{S_{n}\\}$ is independent if

(4) $\text{Pr}(S_{n}=a_{n}|S_{1}^{n-1}=a_{1}^{n-1})=\text{Pr}(S_{n}=a_{n})$

holds for all $n\geq 1$ and all $a_{1}^{n}\in\mathcal{A}^{n}$. It is
identically distributed if $\text{Pr}(S_{n}=a)=\text{Pr}(S_{1}=a)$ holds for
all $n\geq 1$ and for all $a\in\mathcal{A}$. An independent process is
stationary if and only if it is identically distributed. The simplest example
of dependent finite-alphabet processes are the Markov processes.

######  Definition 2.9 (Markov Chain).

A sequence of random variables, $\\{S_{n}\\}$, is a Markov chain if

(5)
$\text{Pr}(S_{n}=a_{n}|S_{1}^{n-1}=a_{1}^{n-1})=\text{Pr}(S_{n}=a_{n}|S_{n-1}=a_{n-1})$

holds for all $n\geq 2$ and all $a_{1}^{n}\in\mathcal{A}^{n}$. A Markov chain
is known as homogeneous or to have stationary transitions if
$\text{Pr}(S_{n}=b|S_{n-1}=a)$ does not depend on $n$, in which case the
$|A|\times|A|$ matrix $M$ defined by

(6) $M_{b|a}=\text{Pr}(S_{n}=b|S_{n-1}=a)$

is called the transition matrix for the chain.

A generalization of the Markov property allows dependence on $k$ steps in the
past. A process is called $k-$step (or $k-$order) Markov if

(7)
$\text{Pr}(S_{n}=a_{n}|S_{1}^{n-1}=a_{1}^{n-1})=\text{Pr}(S_{n}=a_{n}|S_{n-k}^{n-1}=a_{n-k}^{n-1})$

holds for all $n>k$ and for all $a_{1}^{n}$.

The symbolic time-series analysis is initialized by partitioning the phase
space of a dynamical system which is a non-linear mapping from a continuous
space to a discrete space. In machine learning literature, discretization is
generally studied as a feature extraction technique. In dynamical system
literature, the partitioning or discretization is characterized by the extent
to which a dynamical system can be represented by a symbolic one.

The data is partitioned based on the choice of a partitioning technique and
then, the dynamics of the discrete process is studied. The dynamics of the
symbols sequences are compressed as a Probabilistic Finite State Automata
(PFSA), which is defined as follows:

######  Definition 2.10 (PFSA).

A Probabilistic Finite State Automata (PFSA) is a tuple
$T=\\{Q,\mathcal{A},\delta,\bm{\tilde{\Pi}}\\}$ where

  * •

$Q$ is a finite set of states of the automata;

  * •

$\mathcal{A}$ is a finite alphabet set of symbols $a\in\mathcal{A}$;

  * •

$\delta:Q\times\mathcal{A}\rightarrow Q$ is the state transition function;

  * •

$\bm{\tilde{\Pi}}:Q\times\mathcal{A}\rightarrow[0,1]$ is the
$\left|Q\right|\times\left|\mathcal{A}\right|$ emission matrix. The matrix
$\bm{\tilde{\Pi}}=[\tilde{\pi}_{ij}]$ is row stochastic such that
$\tilde{\pi}_{ij}$ is the probability of generating symbol $a_{j}$ from state
$q_{i}$.

######  Remark 2.1.

The $\left|Q\right|\times\left|Q\right|$ state-transition matrix
$\bm{\Pi}:Q\times Q\rightarrow[0,1]$ can be constructed from the state
transition function $\delta$ and the emission matrix $\bm{\tilde{\Pi}}$. The
state transition matrix $\bm{\Pi}$ can be defined as follows.

$\bm{\Pi}_{jk}\triangleq\sum\limits_{a\in\mathcal{A}:\delta(q_{j},a)=q_{k}}\bm{\tilde{\Pi}}(q_{j},a),\forall
q_{j},q_{k}\in Q$

The matrix $\bm{\Pi}=[\pi_{ij}]$ is row stochastic and $\pi_{ij}$ is the
probability $\text{Pr}(q_{j}|q_{i})$ of visiting state $j$ given state $i$ in
the past.

A PFSA can be seen as a generative machine - probabilistically generating a
symbol string through state transitions. Ranging from pattern recognition,
machine learning to computational linguistics, PFSA has found many uses.
Certain well known classes of PFSAs are Hidden Markov Models, stochastic
regular grammars and n-grams. A comprehensive review of PFSA can be found in
[39, 40]. From the perspective of machine learning, the idea is that we infer
the regular patterns in the symbol sequence which help explain the temporal
dependence in the discrete sequence; the patterns we look for are words of
finite length over the alphabet of the discrete process. Once these patterns
are discovered, then the discrete data is compressed as a PFSA whose states
are the words of finite length; this step leads to some information loss. The
parameters of the PFSA are the symbol emission probabilities from each of the
finite length words (states of PFSA). Under the assumption of an ergodic
Markov chain, we relax requirement for initial states and thus the only
sufficient statistics [41, 42] for the inferred Markov model are the
corresponding symbol emission probabilities. For clarification of
presentation, we present the following.

######  Definition 2.11 (Sufficient Statistic).

Let us assume $S_{1},S_{2},\dots,S_{n}\sim p(s;\alpha)$. Then any function
$\mathcal{T}=\mathcal{T}(S_{1},S_{2},\dots,S_{n})$ is itself a random variable
which we will call a statistic. Then, $\mathcal{T}$ is sufficient for $\alpha$
if the conditional distribution of $S_{1},S_{2},\dots,S_{n}\mid\mathcal{T}$
does not depend on $\alpha$. Thus we have $p(s_{1},s_{2},\dots,s_{n}\mid
t;\alpha)=p(s_{1},s_{2},\dots,s_{n}\mid t).$

For the case of finite-order, finite-state Markov chains it means that the
conditional symbol emission probabilities and the initial state summarizes in
entirety the whole of the relevant information supplied by any sample [43].
Under stationarity assumptions, the initial state becomes unnecessary and
thus, the sufficient statistics is provided by the conditional symbol emission
probabilities. The details on estimation of the parameters are provided later
in the paper.

For symbolic analysis of time-series data, a class of PFSAs called the
$0pt$-Markov machine have been proposed [4] as a sub-optimal but
computationally efficient approach to encode the dynamics of symbol sequences
as a finite state machine. The main assumption (and reason for sub-optimality)
is that the symbolic process can be approximated as a $D^{th}$ order Markov
process. For most stable and controlled engineering systems that tend to
forget their initial conditions, a finite length memory assumption is
reasonable. The states of this PFSA are words over $\mathcal{A}$ of length
$0pt$ (or less); and state transitions are described by a sliding block code
of memory $0pt$ and anticipation length of one [44]. The dynamics of this PFSA
can both be described by the $\left|Q\right|\times\left|Q\right|$ state
transition matrix $\bm{\Pi}$ or the $\left|Q\right|\times 1$ state visit
probability vector $\bm{p}$. Next we present definitions of D-Markov machines
which has been recently introduced in literature as finite-memory approximate
models for inference and learning.

######  Definition 2.12.

($D$-Markov Machine [4, 5]) A $D$-Markov machine is a statistically stationary
stochastic process $S=\cdots a_{-1}a_{0}a_{1}\cdots$ (modeled by a PFSA in
which each state is represented by a finite history of $D$ symbols), where the
probability of occurrence of a new symbol depends only on the last $D$
symbols, i.e.,

(8) $P[a_{n}\mid\cdots a_{n-D}\cdots a_{n-1}]=P[a_{n}\mid a_{n-D}\cdots
a_{n-1}]$

where $D$ is called the depth of the Markov machine.

A $D-$Markov machine is thus a $D^{th}$-order Markov approximation of the
discrete symbolic process.

The PFSA model for the time-series data could be inferred as a probabilistic
graph whose nodes are words (or symbol blocks) over $\mathcal{A}$ of length
equal to $0pt$. This probabilistic model induces a Markov model over the
states of the PFSA and the parameters of the Markov model can then be
estimated from the data using a Maximum-likelihood approach. This completes
the process of model inference using STSA.

Next we define some notions from information and communication theory which
will be used in this paper. In particular, we define the notions of Kolmogorov
complexity and Minimum Description length and the related ideas.

######  Definition 2.13.

(Kolmogorov Complexity) The Kolmogorov Complexity $K_{\mathcal{U}}(x)$ of a
string $x$ w.r.t. an universal computer $\mathcal{U}$ is defined as

(9) $K_{\mathcal{U}}(x)=\min\limits_{p:\mathcal{U}(p)=x}\ell(p)$

the minimum length over all programs that print $x$ and halt. Thus,
$K_{\mathcal{U}}(x)$ is the shortest description length of $x$ over all
descriptions interpreted by computer $\mathcal{U}$. The term $\mathcal{U}(p)$
denotes the output of the universal computer $\mathcal{U}$, when presented
with a program $p$.

The Minimum Description length (MDL) [45] principle is a relatively recent
method for inductive inference that provides a generic solution to the model
selection problem. MDL is based on the following insight: any regularity in
data can be used to compress that data i.e., to describe it using fewer
symbols needed to describe the data literally. The more regularities there
are, the more the data can be compressed. Thus equating ’learning’ with
’finding regularity’, we can therefore say that the more we are able to
compress the data, the more we have learned about the data.

Next we present Borel-Cantelli lemma, which is a fundamental result in
probability theory. This is used in the subsequent sections to establish
consistency of the order estimator for Markov chain.

######  Lemma 2.1 (Borel-Cantelli Lemma).

Let $\\{B_{k}\\}$ be an arbitrary sequence of events on a probability space
$(\Omega,\mathcal{E},P)$ and $p_{k}=P(B_{k})$. Then, the following is true

$\sum_{k=1}^{\infty}P(B_{k})<\infty\implies P(\lim\sup B_{k})=0$

where $\lim\sup$ is defined as follows for any sequence of events
$\\{E_{k}\\}$:
$\lim\sup\limits_{n\rightarrow\infty}\triangleq\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}E_{k}$.

###### Proof.

We have that the following is true.

$P(\lim\sup
A_{n})=P(\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}B_{k})\leq\lim\limits_{n\rightarrow\infty}\sum_{k>n}p(B_{k}).$

Since $\sum_{k=1}^{\infty}P(B_{k})<\infty$, the tail end of the series must
sum to zero. Hence $\lim\limits_{n\rightarrow\infty}\sum_{k>n}P(B_{k})=0$.
$\square$∎

Next we introduce an information-theoretic distance which is used in the
subsequent sections to quantify changes in the Markov models during a physical
process.

######  Definition 2.14 (Kullback-Leibler Divergence).

The Kullback-Leibler (K-L) divergence of a discrete probability distribution
$P$ from another distribution $\tilde{P}$ is defined as follows.

$D_{\textrm{KL}}(P\|\tilde{P})=\sum_{x\in
X}{p}_{X}(x)\frac{\log({p}_{X}(x))}{\log(\tilde{p}_{X}(x))}$

It is noted that K-L divergence is not a proper distance as it is not
symmetric. However, to treat it as a distance it is generally converted into
symmetric divergence as follows,
$d(P,\tilde{P})=D_{\textrm{KL}}(P\|\tilde{P})+D_{\textrm{KL}}(\tilde{P}\|P)$.
This is defined as the K-L distance between the distributions $P$ and
$\tilde{P}$.

##  3\. Problem Formulation

As explained in section 2, there are two critical steps for Markov modeling of
data using symbolic dynamics. Once the model structure is fixed, we need to
estimate the parameters of the model from the data and thus, the following
steps are followed during model inference.

  * •

Symbolization or Partitioning

  * •

Order estimation of the discrete stochastic process

  * •

Parameter estimation of the Markov model inferred with the composite method of
discretization and order estimation.

The overall idea of modeling the signal as a Markov chain is also represented
in Figure 1 where the various factors to be considered during Markov modeling
are emphasized. For example, how the original time-series data is related to
the discrete symbol sequence and then, how can we define a metric to relate
the representation in the original space and the discrete space. Another
important question is to analyze the relation of the original data with the
discrete Markov model inferred using the parameters estimated during the
modeling process.

Figure 1. Markov modeling of time series data and the associated decisions of
partitioning and parameter estimation

In general, for most of the machine learning applications, no description of
the underlying system model is available. This in turn makes evaluation of
different models for signal representation very difficult as the topological
structure of the underlying dynamical system is unknown. In the presence of
some end objectives like class separability, anomaly detection, etc. a
Bayesian inference rule could be used for selection of the optimal model [14].
However, the problem is ill-posed in the absence of a well defined end
objective when the aim is to just infer the precise model for signal
representation which can be later used for various different operations like
clustering or classification. The problem of partitioning the data is
particularly ill-posed in the absence of a model, as it is difficult to
specify a metric to define a consistent criterion as the process terminates
with a Markov model of the data when the true statistical nature of data is
unknown. In general, the order of the discrete data depends on the parameters
of the partitioning process i.e., the cardinality of the discrete set and the
location of the partitions. Thus, although, the precision of the Markov model
depends on the partitioning process but in the absence of any model for the
data, it is hard to evaluate performance or find a criterion for consistency.

In this paper, we try to present some results available in open literature on
the mathematical formalism of partitions and order estimation of discrete
stationary Markov processes. In particular, we try to answer the following
questions.

  * •

Problem 1. Given a compact phase space $\Omega\subseteq\mathbb{R}^{n}$ of a
dynamical system $\dot{\mathbf{x}}=f(\mathbf{x})$, under what conditions a
topological Markov shift representation for the dynamical system yield a
partitioning for the system? Furthermore, how are such partitions
characterized?

  * •

Problem 2. Synthesis of an estimator for estimating the order of a discrete
stochastic process that estimates the right order of the Markov process with
probability $1$.

  * •

Problem 3. For practical applications, it is more suitable if the parameters
of the final Markov model could be estimated in a recursive manner. How to
synthesize the recursive update rule for the Markov model?

In the following sections, we try to answer the above questions and then,
demonstrate the proposed methods of data-sets from some engineering systems.
We will also present how the two steps of discretization and order estimation
can be tied together for inferring a Markov model from the time series where
we try to explore the interplay between discretization and order estimation.

##  4\. Partitioning of Phase Space

In this section, we will present some concepts and results on some special
kind of well-behaved partitioning. Most the results in this section are based
on earlier work presented in [8] [44]. Specifically, we will focus on a
special class of partitions called Markov partitions. Loosely speaking, by
using a Markov partition, a dynamical system can be made to resemble a
discrete-time Markov process, with the long-term dynamical characteristics of
the system represented by a Markov shift. Before providing a formal definition
for Markov partitions we need to visit some properties of a partition from
dynamical systems literature. For brevity proofs of some simple properties are
skipped. Interested readers are referred to [8] for more details. It is noted
that in the subsequent material interior of an arbitrary set $\mathcal{X}$ is
denoted by $\mathcal{X}^{o}$.

######  Definition 4.1.

Given two topological partitions $\mathcal{R}=\\{R_{0},R_{1},\dots,R_{n-1}\\}$
and $\mathcal{P}=\\{P_{0},P_{1},\dots,P_{n-1}\\}$, we define their topological
refinement $\mathcal{R}\lor\mathcal{P}$ as

$\mathcal{R}\lor\mathcal{P}=\\{R_{i}\cap
P_{j}:R_{i}\in\mathcal{R},P_{j}\in\mathcal{P}\\}$

Based on the definition 4.1, the following are true.

  * •

The common topological refinement of two topological partitions is a
topological partition.

  * •

For a dynamical system $(M,\phi)$ with topological partition $\mathcal{R}$ of
$M$, the set defined by
$\phi^{n}\mathcal{R}=\\{\phi^{n}R_{0},\phi^{n}R_{1},\dots,\\\
\phi^{n}R_{n-1}\\}$ is also a topological partition.

From the above two results we have that for $m\leq n$,
$\bigvee_{m}^{n}\phi^{k}\mathcal{R}=\phi^{m}\mathcal{R}\lor\phi^{m-1}\mathcal{R}\lor\dots\lor\phi^{n}\mathcal{R}$
is again a topological partition. In the subsequent material, we use the
following notations.

$\mathcal{R}^{n}\equiv\bigvee\limits_{k=0}^{n-1}\phi^{-k}\mathcal{R}$

The diameter of a partition is defined as
$d(\mathcal{R})=\max\limits_{R_{i}\in\mathcal{R}}d(R_{i})$, where
$d(R_{i})\equiv\sup\limits_{x,y\in R_{i}}d(x,y)$. Next, we define cylinder set
partition for symbol sequences which will be useful to characterize sufficient
conditions for existence of Markov partitions for a dynamical system. Consider
a topological Markov shift $(X_{V},\sigma)$ where the vertex set is labeled by
the alphabet $\mathcal{A}$. Then, we form the partition
$\mathcal{C}=\\{C_{a}:a\in\mathcal{A}\\}$ of the elementary cylinder sets
determined by fixing the $0th-$coordinate i.e., $\mathcal{C}_{a}=\\{x\in
X_{V}:x_{0}=a\\}$. Then, this has the following properties.

  1. (1)

$\bigcap_{m=0}^{\infty}\bigcap_{-m}^{m}\phi^{-k}C_{a_{k}}=\\{x\\}$.

  2. (2)

If $x\in C_{a}\cap\sigma^{-1}C_{b}$, then $x\in C_{a}$ and $\sigma x\in
C_{b}$, i.e., $x_{0}=a,x_{1}=b$ (or it means an edge from $a$ to $b$ in terms
of a graph.)

  3. (3)

If $C_{a}\cap\sigma^{-1}C_{b}\neq\emptyset$ and
$C_{b}\cap\sigma^{-1}C_{c}\neq\emptyset$, then
$C_{a}\cap\sigma^{-1}C_{b}\cap\sigma^{-2}C_{c}\neq\emptyset$. This can be
extended to length $n$ for arbitrary $n$. We shall call such a countable set
of conditions for $n=2,3,\dots$ the Markov property: this turns out to be one
of the key properties in getting the desired symbolic representation from a
partition.

######  Definition 4.2 (Generator).

We call a topological partition a generator for a dynamical system $(M,\phi)$
if $\lim\limits_{n\rightarrow\infty}d(\bigvee_{-n}^{n}\phi^{k}\mathcal{R})=0$.

For an expansive dynamical system $(M,\phi)$ let $\mathcal{R}$ be a
topological partition such that $d(\mathcal{R})<c$ where $c$ is the expansive
constant. Then, $\mathcal{R}$ is a generator.

For a topological partition $\mathcal{R}$, let us call $X_{\mathcal{R},\phi}$
the symbolic dynamical representation of $(M,\phi)$. Then, for each $x\in
X_{\mathcal{R},\phi}$ we define
$E_{n}(x)=\bigcap\limits_{k=-n}^{n}\phi^{-k}R_{x_{k}}$. Markov partitions are
defined more formally next.

######  Definition 4.3 (Markov Partitions).

Let $(M,\phi)$ define be an invertible dynamical system. A topological
partition $\mathcal{R}=\\{R_{1},R_{2},\dots,R_{n-1}\\}$ of $M$ gives a
symbolic representation of $(M,\phi)$ if for every $x\in X_{\mathcal{R},\phi}$
the intersection of $\cap_{n=0}^{\infty}\bar{E}_{n}(x)$ contains exactly one
point. We call $\mathcal{R}$ to be a Markov partition for $(M,\phi)$ if
$\mathcal{R}$ gives a symbolic representation of $(M,\phi)$ and furthermore,
$X_{\mathcal{R},\phi}$ is a shift of finite type.

There are natural extensions of Markov partitions for the cases where the
dynamical system is not invertible where we define one-sided symbolic
representation of $(M,\phi)$ which is denoted as $X_{\mathcal{R},\phi}^{+}$.
Next we present a theorem which specifies sufficient condition for existence
of Markov partitions. In the following we prove that the partition
$\mathcal{C}=\\{C_{i}:i=0,1,\dots,n-1$ consisting of elementary cylinder sets
$C_{i}=\\{x\in X_{V}:x_{0}=i\\}$ for a dynamical system $(X_{V},\sigma)$ is a
topological Markov generator.

######  Theorem 4.1.

[8] Let $(M,\phi)$ be a dynamical system, $(X_{G},\pi)$ an irreducible shift
of finite type based on $n$ symbols, and suppose there exists an essentially
one-to-one factor map $\pi$ from $X_{G}$ to $M$. Then the partition
$\mathcal{R}=\\{R_{i}=\pi(C_{i})^{o}:i=0,1,\dots,n-1\\}$ is a topological
Markov partition.

###### Proof.

We need to prove the following in order to prove that $\mathcal{R}$ is the
Markov partition.

  1. (1)

Elements of $\mathcal{R}$ are disjoint.

  2. (2)

The closure of elements of $\mathcal{R}$ cover $M$.

  3. (3)

$R$ is a generator.

  4. (4)

$\mathcal{R}$ satisfies the Markov property.

In the following we prove the four items listed above in order to prove that
$\mathcal{R}$ is a Markov partitions.

  1. (1)

Elements of $\mathcal{R}$ are disjoint, i.e., $R_{i}\cap R_{j}=\emptyset$,
$i\neq j$.
The idea of the proof is to use bilaterally transitive points to overcome a
difficulty that, in general, maps do not enjoy the property that the image of
an intersection is equal to the intersection of images, but one-to-one maps
do. Suppose $R_{i}\cap R_{j}\neq 0$ for $i\neq j$. Then, by proposition ..
$BLT(R_{i}\cap R_{j})\neq\emptyset$. By proposition $4.8$ and $4.9$,
$\pi^{-1}$ maps $BLT(R_{i})$ and $BLT(R_{j})$ homeomorphically onto
$BLT(C_{i})$ and $BLT(C_{j})$ respectively. Therefor $\pi^{-1}$ maps
$BLT(R_{i}\cap R_{j})=BLT(R_{i})\cap BLT(R_{j})$ homeomorphically onto
$BLT(C_{i}\cap C_{j})=BLT(C_{i})\cap BLT(C_{j})$, which implies that
$\emptyset\neq BLT(C_{i}\cap C_{j})\subset C_{i}\cap C_{j}$, which is a
contradiction.

  2. (2)

$M=\cup_{i=0}^{n-1}\bar{R}_{i}.$
$M=\pi(\Sigma_{G})$ = $\pi(\cup_{i=0}^{n-1}C_{i})$ =
$\cup_{i=0}^{n-1}\pi(C_{i})$ = $\cup_{i=0}^{n-1}\bar{R}_{i}$, the last
equality follows from proposition …

To prove the next two items, we need another lemma which is presented next.

######  Lemma 4.1.

Under the hypothesis of Theorem 4.1,
$\pi(\cap_{m}^{n}\sigma^{-k}C_{a_{k}})=\overline{\cap_{m}^{n}\phi^{-k}R_{a_{k}}}$.

###### Proof.

Once again we use the bilaterally transitive points to deal with images of
intersections. We have the following equalities.

$\pi(\bigcap_{m}^{n}\sigma^{-k}C_{a_{k}})=\pi\big{(}\overline{BLT(\bigcap_{m}^{n}\sigma^{-k}C_{a_{k}})}\big{)}=\overline{\pi(\bigcap_{m}^{n}BLT(\sigma^{-k}C_{a_{k}})\big{)}}$

which by injectivity of $\pi$ and shift-invariance of bilateral transitive
points.

$=\overline{\bigcap_{m}^{n}\pi(BLT(\sigma^{-k}C_{a_{k}}))}=\overline{\bigcap_{m}^{n}\pi\sigma^{-k}BLT(C_{a_{k}})}$

which by commutativity of $\pi$ can be reduced to the following.

$=\overline{\bigcap_{m}^{n}\phi^{-k}\pi(BLT(C_{a_{k}}))}=\overline{\bigcap_{m}^{n}\phi^{-k}(BLT(\pi(C_{a_{k}})))}$

This can be further reduced to the following.

$=\overline{\bigcap_{m}^{n}\phi^{-k}(BLT(R_{a_{k}}))}=\overline{BLT(\bigcap_{m}^{n}\phi^{-k}(R_{a_{k}}))}=\overline{\bigcap_{m}^{n}\phi^{-k}(R_{a_{k}})}.$

$\square$∎

  3. (3)

$\mathcal{R}$ is a generator.
Because $\mathcal{C}$ is a generator,
$d(\cap_{-n}^{n}\sigma^{-k}C_{a_{k}})\rightarrow 0$. By Lemma,
$\pi(\cap_{-n}^{n}\sigma^{-k}C_{a_{k}})=\overline{\cap_{-n}^{n}\phi^{-k}R_{a_{k}}}$.
So, by continuity of $\pi$, we get.

(10)
$d\big{(}\overline{\cap_{-n}^{n}\phi^{-k}R_{a_{k}}}\big{)}=d\big{(}\cap_{-n}^{n}\phi^{-k}R_{a_{k}}\big{)}\rightarrow
0$

This proves that $\mathcal{R}$ is a generator.

  4. (4)

$\mathcal{R}$ satisfies the Markov property.
Suppose $R_{a_{i}}\cap\phi^{-1}R_{a_{i}+1}\neq\emptyset$, $1\leq k\leq n-1$.
By Lemma, we have
$\pi[C_{a_{i}}\cap\sigma^{-1}C_{a_{i}+1}]=\overline{R_{a_{i}}\cap\phi^{-1}R_{a_{i}+1}}\neq\emptyset$,
$1\leq k\leq n-1$. Thus $C_{a_{i}}\cap\sigma^{-1}C_{a_{i}+1}\neq\emptyset$,
$1\leq k\leq n-1$. Since $\mathcal{C}$ satisfies the Markov property,
$\cap_{k=1}^{n}\phi^{-k}C_{a_{k}}\neq\emptyset$ for all $n>1$. Therefore,
$\cap_{k=1}^{n}\phi^{-k}R_{a_{k}}\neq\emptyset$ for all $n>1$.

This completes the proof of this theorem. $\square$∎

The above theorem establishes the a sufficient condition for existence of
Markov partitions. However, except for some well-behaved dynamical systems,
there are no results on existence of Markov partitioning for a dynamical
system. Thus, for most of the machine learning applications we are interested
in, it would be very difficult to estimate the partitioning which can be
proved to induce Markov dynamics. However, given that the dynamics of a
discrete stochastic process is Markovian, we need to be able to estimate the
order of the Markovian dynamics. This is discussed in the next section.

##  5\. Order Estimation of Markov Chains

In this section, we will present an order estimation technique for Markov
chains which converges almost surely to $k$ given that the underlying
stationary and ergodic stochastic process is a $k^{th}$ order Markov chain and
infinity otherwise. The work on order estimation presented in this section is
based on earlier work presented in [28, 29, 30, 31].

Let $\\{A\\}_{n=0}^{\infty}$ be a stationary and ergodic time-series taking
values from a discrete (finite or infinite) alphabet $\mathcal{A}$. In order
to estimate the order we need to define some explicit statistics. For $k\geq
0$ let $\mathcal{S}_{k}$ denote the support of the distribution of
$S_{-k}^{0}$ as

$\displaystyle\mathcal{S}_{k}=\\{a_{-k}^{0}\in\mathcal{A}^{k+1}:p(a_{-k}^{0})>0\\}$

Next we define

$\displaystyle\Delta_{k}=\sup\limits_{1\leq
i}\sup\limits_{(b_{-k-i+1}^{0},a)\in\mathcal{S}_{k+i}}|p(a|b_{-k+1}^{0})-p(a|b_{-k-i+1}^{0})|$

We divide the data segment $S_{0}^{n}$ into two parts:
$S_{0}^{\lceil\frac{n}{2}\rceil-1}$ and $S_{\lceil\frac{n}{2}\rceil}^{n}$. Let
$\mathcal{S}_{n,k}^{(1)}$ denote the set of strings with length $k+1$ which
appear at all in $S_{0}^{\lceil\frac{n}{2}\rceil-1}$, i.e.,

$\mathcal{S}_{n,k}^{(1)}=\\{a_{-k}^{0}\in\mathcal{A}^{k+1}:\exists k\leq
t\leq\lceil\frac{n}{2}\rceil-1:S_{t-k}^{t}=a_{-k}^{0}\\}$

For a fixed $0<\gamma<1$, let $\mathcal{S}_{n,k}^{(2)}$ denote the set of
strings with length $k+1$ which appear more than $n^{1-\gamma}$ times in
$S_{\lceil\frac{n}{2}\rceil}^{n}$. That is,

$\mathcal{S}_{n,k}^{(2)}=\\{a_{-k}^{0}\in\mathcal{A}^{k+1}:N\\{\lceil\frac{n}{2}\rceil+k\leq
t\leq\lceil\frac{n}{2}\rceil-1:S_{t-k}^{t}=a_{-k}^{0}\\}>n^{1-\gamma}\\}$

where $N\\{\bullet\\}$ denotes the count of $\bullet$. Let

$\mathcal{S}_{k}^{n}=\mathcal{S}_{n,k}^{(1)}\cap\mathcal{S}_{n,k}^{(2)}$

For the sake of notational convenience, let $C(a\mid
b_{-k+1}^{0}:[n_{1},n_{2}])$ denote the empirical conditional probability of
$A_{1}=a$ given $A_{-k+1}^{0}=b^{0}_{-k+1}$ from the samples
$(A_{n_{1}},\dots,A_{n_{2}})$, that is,

(11) $C(a\mid b_{-k+1}^{0}:[n_{1},n_{2}])=\frac{N\\{n_{1}+k\leq k\leq
n_{2}:A_{t-k}^{t}=(b_{-k+1}^{0},a)\\}}{\\{n_{1}+k-1\leq t\leq
n_{2}-1:A_{t-k+1}^{t}=b_{-k+1}^{0}\\}}$

where $0/0$ is defined as $0$.

We define the empirical version of $\Delta_{k}$ as follows:

(12) $\hat{\Delta}_{k}^{n}=\max\limits_{1\leq i\leq
n}\max\limits_{(b_{-k-i+1}^{0},a)\in\mathcal{S}_{k+i}^{n}}\mid C(a\mid
b_{-k+1}^{0}:[\lceil\frac{n}{2}\rceil,n])-C(a\mid
b_{-k-i+1}^{0}:[\lceil\frac{n}{2}\rceil,n])\mid$

Then we observe that by ergodicity, for any fixed $k$, the following is true.

(13) $\lim\inf\limits_{n\rightarrow\infty}\hat{\Delta}_{k}^{n}\geq\Delta_{k}$

######  Definition 5.1 (Markov Chain Order Estimator).

We define an estimate $D_{n}$ for the order of samples from $S_{0}^{n}$ as
follows. Let the $0<\beta<\frac{1-\gamma}{2}$ be arbitrary. Set $D_{0}=0$, and
for $n\geq 1$ let $D_{n}$ be the smallest $0\leq k_{n}<n$ such that
$\hat{\Delta}^{n}_{k_{n}}\leq n^{-\beta}$.

The consistency of the Markov chain estimator defined in Definition 5.1 is
proved in a theorem which is presented next.

######  Theorem 5.1.

[28] If the stationary and ergodic discrete time series $\\{S_{n}\\}$ taking
values from a discrete alphabet happens to be Markov with any finite order
then $D_{n}$ equals the order eventually almost surely, and if it is not
Markov with any finite order then, $D_{n}\rightarrow\infty$ almost surely.

###### Proof.

If the process is Markov, it is immediate that for all $k$ greater than or
equal to the order $\Delta_{k}=0$. For $k$ less than the order $\Delta_{k}>0$.
If the process is not a Markov chain with any finite order then $\Delta_{k}>0$
for all $k$. This by (13) if the process is not Markov then
$D_{n}\rightarrow\infty$ and if it is Markov then $D_{n}$ is greater or equal
to the order eventually almost surely. We have to show that $D_{n}$ is less or
equal the order eventually almost surely provided that the process is a Markov
chain.

Assume that the process is a Markov chain with order $k$. Let $n\geq k$. We
will estimate the probability of the undesirable event as follows.

$\displaystyle\text{Pr}(\hat{\Delta}_{k}^{n}>n^{-\beta}|S_{0}^{\lceil\frac{n}{2}\rceil})\leq$
(14)
$\displaystyle\sum_{i=1}^{n}\text{Pr}(\max\limits_{(b_{-k-i+1,a}\in\mathcal{S}_{k+i}^{n})}\mid
C(a\mid b_{-k+1}^{0}:[\lceil\frac{n}{2}\rceil,n])-C(a\mid
b_{-k-i+1}^{0}:[\lceil\frac{n}{2},n])\mid>n^{-\beta}\mid
S_{0}^{\lceil\frac{n}{2}\rceil})$

We can estimate each probability in the sum as the sum of the following two
terms:

$\displaystyle\sum_{i=1}^{n}\text{Pr}(\max\limits_{(b_{-k-i+1,a}\in\mathcal{S}_{k+i}^{n})}\mid
C(a\mid b_{-k+1}^{0}:[\lceil\frac{n}{2}\rceil,n])-C(a\mid
b_{-k-i+1}^{0}:[\lceil\frac{n}{2}\rceil,n])\mid>n^{-\beta}\mid
S_{0}^{\lceil\frac{n}{2}\rceil})$
$\displaystyle\leq\text{Pr}(\max\limits_{(b_{-k-i+1}^{0},a)\in\mathcal{S}_{k+i}^{n}}\mid
C(a\mid b_{-k+1}^{0}:[\lceil\frac{n}{2}\rceil,n])-p(a\mid
b_{-k+1}^{0}))\mid>0.5n^{-\beta}\mid S_{0}^{\lceil\frac{n}{2}\rceil})$
$\displaystyle+\text{Pr}(\max\limits_{(b_{-k-i+1}^{0},a)\in\mathcal{S}_{k+i}^{n}}\mid
p(a\mid b_{-k+1}^{0})-C(a\mid
b_{-k-i+1}^{0}:[\lceil\frac{n}{2}\rceil,n])\mid>0.5n^{-\beta}\mid
S_{0}^{\lceil\frac{n}{2}\rceil})$

We overestimate these probabilities. For any $m\geq 0$ and $a_{-m}^{0}$ define
$\sigma_{i}^{m}(a_{-m}^{0})$ as the time of the $i-$th occurrence of the
string $a_{-m}^{0}$ in the data segment $S_{\lceil\frac{n}{2}\rceil}^{n}$,
that is, let $\sigma_{0}^{m}(a_{-m}^{0})=\lceil\frac{n}{2}\rceil+m-1$ and
$i\geq 1$ define,

$\sigma_{i}^{m}(a_{-m}^{0})=\min\\{t>\sigma_{i-1}^{m}(a_{-m}^{0}):S_{t-m}^{t}:a_{-m}^{0}\\}$

Then,

$\displaystyle\text{Pr}(\max\limits_{(b_{-k-i+1}^{0},a)\in\mathcal{S}_{k+i}^{n}}\mid
C(a\mid b_{-k+1}^{0}:[\lceil\frac{n}{2}\rceil,n])-C(a\mid
b_{-k-i+1}^{0}:[\lceil\frac{n}{2},n])\mid>n^{-\beta}\mid
S_{0}^{\lceil\frac{n}{2}\rceil})$
$\displaystyle\leq\text{Pr}(\max\limits_{(b_{-k+1}^{0},a)\in\mathcal{S}_{n,k}^{(1)}}\sup\limits_{j>n^{1-\gamma}}\mid\frac{1}{j}\sum_{r=1}^{j}1_{\\{S_{\sigma_{r}^{k-1}(b_{-k+1}^{0})}=a\\}}-p(a\mid
b_{-k+1}^{0})\mid>0.5n^{-\beta}\mid S_{0}^{\lceil\frac{n}{2}\rceil})$
$\displaystyle+\text{Pr}(\max\limits_{(b_{-k-i+1}^{0},a)\in\mathcal{S}_{n,k+i}^{(1)}}\sup\limits_{j>n^{1-\gamma}}\mid\frac{1}{j}\sum_{r=1}^{j}1_{\\{S_{\sigma_{r}^{k+i-1}(b_{-k-i+1}^{0})}=a\\}}-p(a\mid
b_{-k+1}^{0})\mid>0.5n^{-\beta}\mid S_{0}^{\lceil\frac{n}{2}\rceil})$

Since both $S_{n,k}^{1}$ and $S_{n,k+i}^{1}$ depend solely on
$S_{0}^{\lceil\frac{n}{2}\rceil}$, we get

$\displaystyle\text{Pr}\big{(}\max\limits_{(b_{-k-i+1}^{0},a)\in
S_{k+i}^{n}}\mid C(a\mid b_{-k+1}^{0}:[\lceil\frac{n}{2}\rceil,n])-C(a\mid
b_{-k-i+1}^{0}:[\lceil\frac{n}{2}\rceil,n])\mid>n^{-\beta}\mid
S_{0}^{\lceil\frac{n}{2}\rceil}\big{)}$
$\displaystyle\leq\sum_{(b_{-k+1}^{0},a)\in S_{n,k}^{(1)}}\sum_{j=\lceil
n^{1-\gamma}\rceil}^{\infty}\text{Pr}(\frac{1}{j}\sum_{r=1}^{j}1_{\\{S_{\sigma_{r}^{k-1}(b_{-k+1}^{0})=a}\\}}-p(a\mid
b_{-k+1}^{0})\mid>0.5n^{-\beta}\mid S_{0}^{\lceil\frac{n}{2}\rceil})$
$\displaystyle+\sum_{(b_{-k-i+1},a)\in\mathcal{S}_{n,k+i}^{(1)}}\sum_{j=\lceil
n^{1-\gamma}\rceil}^{\infty}\text{Pr}(\mid\frac{1}{j}\sum_{r=1}^{j}1_{\\{S_{\sigma_{r}^{k+i-1}}(b_{-k-i+1}^{0})=a\\}}-p(a\mid
b_{-k+1}^{0})\mid>0.5n^{-\beta}\mid S_{0}^{\lceil\frac{n}{2}\rceil})$

Each of these terms represents the deviation of an empirical count from its
mean. The variables under consideration are independent since whenever the
block $b_{-k+1}^{0}$ occurs the next term is chosen using the same
distribution $p(a\mid b_{-k+1}^{0})$. Thus by Hoeffding’s inequality for sums
of bounded independent random variables and since the cardinality of both
$\mathcal{S}_{n,k}^{(1)}$ and $S_{n,k+i}^{(1)}$ is not greater than $(n+2)/2$,
we have the following.

$\displaystyle\text{Pr}\big{(}\max\limits_{(b_{-k-i+1}^{0},a)\mathcal{S}_{k+i}^{n}}\mid
C(a\mid b_{-k+1}^{0}:[\lceil\frac{n}{2}\rceil,n])-C(a\mid
b_{-k-i+1}^{0}:[\lceil\frac{n}{2}\rceil,n])\mid>n^{-\beta}\mid
S_{0}^{\lceil\frac{n}{2}\rceil}\big{)}$ $\displaystyle\leq
2\frac{n+2}{2}\sum\limits_{j=\lceil n^{1-\gamma}\rceil}^{\infty}2e^{-n-2\beta
j}$

Thus we obtain the following.

$\text{Pr}(\hat{\Delta}_{k}^{n}>n^{-\beta}\mid
S_{0}^{\lceil\frac{n}{2}\rceil})\leq n(n+2)4e^{-2n-2\beta+1-\gamma}$

Integrating both sides we get

(15) $\text{Pr}(\hat{\Delta}_{k}^{n}>n^{-\beta}\leq
n(n+2)4e^{-2n-2\beta+1-\gamma}$

The right hand side of equation (15) is summable provided $2\beta+\gamma<1$
and the Borel-Cantelli Lemma (see Lemma 2.1) yields that
$\text{Pr}(\hat{\Delta}\leq n^{-\beta})=1$. Thus $D_{n}\leq k$ eventually
almost surely provided the process is Markov with order $k$. The proof is now
complete. $\square$∎

The above theorem establishes a consistent estimator for the order of an
unknown discrete Markov process. Once the order of the discrete Markov process
is estimated, we estimate the sufficient statistics for the Markov chain from
the data. This is described in detail in next section.

##  6\. Parameter Estimation for Markov Chains

In this section, we present some results on estimation of parameters of the
Markov chain once the structure of the chain is inferred after discretization
and order estimation. Given a finite-length symbol string over a (finite)
alphabet $\mathcal{A}$, there exists several PFSA construction algorithms to
discover the underlying irreducible PFSA model, such as the D-Markov machines
[4][10]. Once the order of the discrete data is estimated, the states of a
D-Markov machines are memory words of length equal to the order over
$\mathcal{A}$. The sufficient statistics for the D-Markov machine are the
symbol emission probabilities conditioned on the memory words. These
statistics could be estimated using a Bayesian approach. In general, the
symbol emission probabilities conditioned on the individual memory words can
be modeled as a multinomial random variable which can be modeled by a
Dirichlet distribution. However, we skip those details here and present a very
simple estimation process where every random variable has a uniform prior.

Let $N(a_{j}|q_{k})$ denote the number of times that a symbol $a$ is generated
from the state $q_{k}$ as the symbol string evolves. The maximum a posteriori
probability (MAP) estimate of the morph probability (see Definition 2.10) for
the PFSA $T$ is estimated by frequency counting as

(16)
$\hat{\pi}(a_{j}|q_{k})\triangleq\frac{C(a_{j}|q_{k})}{\sum_{\ell}C(a_{\ell}|q_{k})}=\frac{1+N(a_{j}|q_{k})}{|\mathcal{A}|+\sum_{\ell}N(a_{\ell}|q_{k})}$

The rationale for initializing each element of the count matrix $C$ to $1$ is
that if no event is generated at a state $q\in Q$, then there should be no
preference to any particular symbol and it is logical to have
$\hat{\pi}_{MAP}(q,a)=\frac{1}{|\mathcal{A}|}\ \forall a\in\mathcal{A}$ (i.e.,
the uniform distribution of event generation at the state $q$). The above
procedure guarantees that the PFSA, constructed from a (finite-length) symbol
string, must have an (elementwise) strictly positive morph map $\pi$.

Having computed the probabilities $\hat{\pi}\left(a_{j}|q_{k}\right)$ for all
$a_{j}\in\mathcal{A}$ and $q_{k}\in Q$, the estimated emission probability
matrix of the PFSA is obtained as

(17)
$\widehat{\mathbf{\Pi}}\triangleq\left[\begin{array}[]{lll}\hat{\pi}\left(a_{1}\mid
q_{1}\right)&\dots&\hat{\pi}\left(a_{|\mathcal{A}|}\mid q_{1}\right)\\\
\vdots&\ddots&\vdots\\\ \hat{\pi}\left(a_{1}\mid
q_{|Q|}\right)&\cdots&\hat{\pi}\left(a_{|\mathcal{A}|}\mid
q_{|Q|}\right)\end{array}\right].$

The stochastic matrix $\widehat{\mathbf{\Pi}}$, estimated from a symbol
string, is then treated as a feature of the times series data that represents
the behavior of the dynamical system. These features can then be used for
different learning applications like pattern matching, classification, and
clustering. The estimated stochastic matrices $\widehat{\mathbf{\Pi}}$ are
converted into row vectors for feature divergence measurement by sequentially
concatenating the $|Q|$ rows (i.e., the $|Q|\times|\Sigma|$ matrix
$\widehat{\mathbf{\Pi}}$ is vectorized as a $(1\times|Q||\mathcal{A}|)$
vector). The convergence of the parameters of the Markov model could be
established by Glivenko-Cantelli theorem [46, 47] which guarantees the
convergence of empirical distribution of random variables under independent
observations.

######  Remark 6.1.

As compared to the Hidden Markov model-based approaches, this approach
presents a much simpler modeling technique where the structure of the model is
inferred based on the observations. As the structure of the Markov model is
also fixed by the observations, this makes inference of the parameters of the
Markov chain much more easier than the standard Hidden Markov model approaches
like the Viterbi algorithm [2] where Dynamic programming is used to estimate
the parameters that maximize the likelihood of the observations.

The sufficient statistics for the inferred Markov process serve as a compact
stochastic representation of the signal and can be used for different machine
learning applications like pattern matching, classification, clustering,
anomaly detection. As explained earlier, the symbol emission probabilities
conditioned on the memory words for the discrete process could be modeled as
random variables. A metric to measure the information gain over the marginal
symbol emission probabilities could be defined as follows:

(18) $\mathfrak{d}_{\theta}=\sum_{q\in Q}\text{Pr}(q)D_{KL}(P(\mathcal{A}\mid
q)\|\tilde{P}({\mathcal{A}}))$

where, $\theta$ represents the parameters of the Markov model inference, i.e.,
the partitioning map and the estimated order of the memory words. The term in
equation (18) measures the measures the discrepancy in the statistics of the
symbol emission probabilities when they are conditioned on the memory words as
compared to the unconditional statistics for the symbol emission
probabilities.

##  7\. Statistical Learning Applications

In this section, we present a case study using the present framework for
statistical learning for applications like anomaly detection, classification
and prognostics in a complex engineering system. While there are several
examples presented in open literature [48, 49, 50, 51, 52, 53, 54, 55, 56], we
will only show one example here for clarity. We present an example of complex
thermoacoustic instability in gas turbine engines during combustion.
Combustion instability is a highly nonlinear and complex phenomena which
results in severe structural degradation in jet turbine engines. Some good
surveys on the current understanding of the mechanisms for the combustion
instability phenomena could be found in [57, 58, 59, 60, 61]. Active
combustion instability control (ACIC) with fuel modulation has proven to be an
effective approach for reducing pressure oscillations in combustors [62, 63].
Based on the work available in literature, one can conclude that the
performance of ACIC is primarily limited by the large delay in the feedback
loop and the limited actuator bandwidth [62, 63]. From the perspective of
active control of the unstable phenomena, it is necessary to accurately detect
and, desirably, predict the states of the combustion process. In this paper,
we present Markov modeling of pressure time-series during combustion and
present results on changes in the underlying model structure and complexity as
the process undergoes some changes in its physical properties. The goal is to
be able to design a statistical model for the instability phenomenon during
combustion which could be used to design a statistical filter to accurately
predict with high confidence the system states. This can potentially alleviate
the problems with delay in the ACIC feedback loop and thus possibly improve
the performance. We first present some experimental details of the set-up that
was used to collect the experimental data and then, show the Markov modeling
of the pressure time-series data which is able to capture the changes in
system behavior.

###  7.1. Experimental Details

In this section we present the experimental details for collecting data to
analyze the complex non-linear phenomena that occurs during the instability
phenomena, in a laboratory-scale combustor. A swirl-stabilized, lean-premixed,
laboratory-scale combustor was used to perform the experimental study. Figure
2 shows a schematic drawing of the variable-length combustor. The combustor
consists of an inlet section, an injector, a combustion chamber, and an
exhaust section. The combustor chamber consists of an optically-accessible
quartz section followed by a variable length steel section.

Figure 2. Schematic drawing of test facility Table 1. Operating conditions Parameters | Value
---|---
Equivalence Ratio | 0.525, 0.55, 0.60, 0.65
Inlet Velocity | 25-50 m/s in 5 m/s increments
Combustor Length | 25-59 inch in 1 inch increments

High pressure air is delivered to the experiment from a compressor system
after passing through filters to remove any liquid or particles that might be
present. The air supply pressure is set to 180 psig using a dome pressure
regulator. The air is pre-heated to a maximum temperature of
$250\,^{\circ}{\rm C}$ by an 88kW electric heater. The fuel for this study is
natural gas (approximately $95$% methane). It is supplied to the system at a
pressure of 200 psig. The flow rates of the air and natural gas are measured
by thermal mass flow meters. The desired equivalence ratio and mean inlet
velocity is set by adjusting these flow rates with needle valves. For fully
pre-mixed experiments (FPM), the fuel is injected far upstream of a choke
plate to prevent equivalence ratio fluctuations. For technically pre-mixed
experiments (TPM), fuel is injected in the injector section near the swirler.
It mixes with air over a short distance between the swirler and the injector
exit. Tests were conducted at a nominal combustor pressure of 1 atm over a
range of operating conditions, as listed in Table 1.

In each test, the combustion chamber dynamic pressure and the global OH and CH
chemiluminescence intensity were measured to study the mechanisms of
combustion instability. The measurements were made simultaneously at a
sampling rate of 8192 Hz (per channel), and data were collected for 8 seconds,
for a total of 65536 measurements (per channel).

###  7.2. Markov Modeling and Results

In this section, we present details of the analyses done using the pressure
time-series data to infer the underlying Markov model. Time-series data is
first normalized by subtracting the mean and dividing by the standard
deviation of its elements; this step corresponds to bias removal and variance
normalization. Data from engineering systems is typically oversampled to
ensure that the underlying dynamics can be captured. Due to coarse-graining
from the symbolization process, an over-sampled time-series may mask the true
nature of the system dynamics in the symbolic domain (e.g., occurrence of self
loops and irrelevant spurious transitions in the Markov chain). Time-series is
first down-sampled to find the next crucial observation. The first minimum of
auto-correlation function generated from the observed time-series is obtained
to find the uncorrelated samples in time. The data sets are then down-sampled
by this lag. To avoid discarding significant amount of data due to
downsampling, down-sampled data using different initial conditions is
concatenated. Further details of this preprocessing can be found in [36].

The continuous time-series data set is then partitioned using maximum entropy
partitioning (MEP), where the information rich regions of the data set are
partitioned finer and those with sparse information are partitioned coarser.
In essence, each cell in the partitioned data set contains (approximately)
equal number of data points under MEP. A ternary alphabet with
$\mathcal{A}=\\{0,1,2\\}$ has been used to symbolize the continuous combustion
instability data. As discussed in section 7.1, we analyze data sets from
different phases, as the process goes from stable through the transient to the
unstable region (the ground truth is decided using the RMS-values of
pressure).

In figure 3a, we show the observed changes in the behavior of the data as the
combustion operating condition changes from stable to unstable. A change in
the empirical distribution of data from unimodal to bi-modal is observed as
the system moves from stable to unstable. We selected $150$ samples of
pressure data from the stable and unstable phases each to analyze and compare.
The temporal memory of the individual series is estimated by the spectral
decomposition method presented earlier in [36, 37]. First, we compare the
expected size of temporal memory during the two stages of operation. There are
changes in the Eigen value decomposition rate for the 1-step stochastic matrix
calculated from the data during the stable and unstable behavior, irrespective
of the combustor length and inlet velocity. During stable conditions, the
Eigen values very quickly go to zero as compared to the unstable operating
condition (see Figure 3b). This suggests that the size of temporal memory of
the discretized data increases as we move to the unstable operating condition.
This indicates that under the stable operating condition, the discretized data
behaves as symbolic noise as the predictive power of Markov models remain
unaffected even if we increase the order of the Markov model. On the other
hand, the predictive power of the Markov models can be increased by increasing
the order of the Markov model during unstable operating condition, indicating
more deterministic behavior. An $\epsilon=0.05$ is chosen to estimate the
depth of the Markov models for both the stable and unstable phases.
Correspondingly, the depth was calculated as $2$ and $3$ for the stable and
unstable conditions (see Figure 3). The corresponding $D(\epsilon)$ is used to
construct the Markov models next. First a PFSA whose states are words over
$\mathcal{A}$ of length $0pt(\epsilon)$ is created and the corresponding
maximum-likely parameters ($\bm{\tilde{\Pi}}$ and $\bm{\Pi}$) are estimated.

(a) Probability density function for the pressure time series data

(b) Spectral decomposition of the stochastic matrix for 1-step Markov model

Figure 3. The first plate in the above figure shows the change in the
empirical density calculated for the pressure time-series data as the process
deviates from the stable operating condition to unstable operating condition.
The second plate shows the spectral decomposition of the 1-step stochastic
matrix for the data under stable and unstable operating conditions. Figure 4.
Class-separability of the stable and unstable classes independent of any other
related variables.

Once the model structure is inferred and the model parameters are estimated
from the data, we estimate certain model-specific metrics to see changes in
the inferred models. In particular, we are interested in inferring the changes
in the model complexity as the process moves from stable to unstable. The size
of temporal memory can also be treated as a metric for complexity of the
underlying Markov model. However, we define another metric based on the KL-
distance between the states of the Markov model. In particular, we estimate
the following.

(19) $d_{\mathcal{M}}=\max\limits_{q_{i},q_{j}\in Q}d(q_{i},q_{j})$

where
$d(q_{i},q_{j})=D_{\textrm{KL}}(\text{Pr}(\mathcal{A}|q_{i})\|\text{Pr}(\mathcal{A}|q_{j}))+D_{\textrm{KL}}(\text{Pr}(\mathcal{A}|q_{j})\|\text{Pr}(\mathcal{A}|q_{i}))$,
i.e., the symmetric KL distance between states $q_{i}$ and $q_{j}$ based on
the conditional symbol emission probabilities. In equation (19), we measure
the maximum divergence in the set of states; however, another possible metric
could be an expectation over the set of states. The behavior of the metric
described in equation (19) is shown in figure 4. It is clear that the proposed
metric is able to achieve a clear separation of the two classes of interest
(see Figure 4).

Next we show the gradual change in the model characteristics and compare it
with the behavior of the RMS of pressure during this process. It is noted at
this point that all the analysis is done in an unsupervised fashion, so this
falls under the broad category of unsupervised learning or anomaly detection.
As shown in Figure 5, we show the different behaviors of the model as the
system moves from the stable to unstable behavior. It shows presence of a
transient phase before sticking to the unstable region. We present results
which is calculated using equation (18) in figure 6. It is interesting to see
that the symbol emission probabilities are quite close to the conditionals
during the stable operating condition, indicating symbolic noise behavior.
While gradually there is an observed discrepancy between the conditional and
marginal distribution of the symbol emissions that saturates as the system
moves to the unstable behavior. Thus the current formalism is able to infer
the underlying changes in the model structures as well as the parameters in an
unsupervised fashion. This is encouraging as these changes can be associated
with the changes in the physical dynamics during the complex process and thus,
gives a high-fidelity statistical model for the process.

(a) Chnages in the complexity of the model inferred from data

(b) The RMS of the pressure time series corresponding to the different models
showing the range of stable and unstable behavior

Figure 5. Change in the Markov models of the pressure time-series data during
the combustion process Figure 6. Change in the discrepancy of the symbol
emission probabilities indicating the information gain by the Markov models

The results presented in the above plots are inspiring as we see that the
models corresponding to different behaviors are clustered together,
irrespective of the other variables (like equivalence ratio, length of
combustor, etc.). Thus the Markov models are able to capture features
intrinsic to the governing physical dynamics independent of the other
variables and thus serves as a good representation of the pressure signals
obtained during the process.

##  8\. Conclusions and Future Work

In this paper we presented the concept of Markov modeling of time-series data
using Symbolic analysis. Symbolic analysis-based Markov modeling is a recently
introduced statistical modeling approach where the data is first discretized
and then, the discrete time-series data is approximated as a finite order
Markov model. Compared to the widely-used Hidden Markov models, the present
concept is algorithmically simple and thus could be easily inferred for
machine learning applications. This is specially useful for embedded
applications where using Hidden Markov model inference methods like Viterbi
algorithm is computationally expensive and thus might be infeasible to use.

We discussed that the efficacy of the proposed approach depends on the
characteristics of the two processes of discretization and Markov modeling of
discrete processes. In this paper, we focused on the two critical steps in the
modeling process namely, discretization and order estimation. We visited the
key concepts of Markov partitioning from the dynamical systems literature and
discussed some properties of the same. We also presented an order estimation
technique for discrete symbolic process and showed the consistency of the
estimation which guarantees almost sure convergence to the true order of the
Markov chain. Use of current mathematical formalism for partitioning or
discretization of data with unknown model is a very challenging problem. The
results in dynamical systems theory is insufficient for use in machine
learning applications. On the other hand, there is no rigorous statistical
analysis of the discretization problem available in open literature. The order
estimation problem for discrete stochastic processes is more rigorously
studied in information theory literature. Consistency of some approaches have
been rigorously established. However, for the statistical learning problems
using STSA framework, the composite process of discretization and modeling of
discrete process needs to be studied together. We also presented a case study
of statistical learning for prognostics and anomaly detection in combustion
process which is a complex thermoacoustic phenomena with undesirable effects.
The Markov modeling approach is able to identify and represent changes in the
signals of pressure time-series as the underlying physical process undergoes
some change.

One possible direction for research is to find a discretization technique
which results in order one discrete Markov process. Proving such a
discretization always exists might be possible under stationarity conditions
for systems with fading memory. For example, the metric introduced in equation
(18) provides a measure to describe the discrepancy between the independent
statistics of the symbols versus the statistics when conditioned on the states
of a probable Markov chain. Maximizing such a discrepancy provides a way to
maximize the information gain obtained by $1^{st}$-order Markov model for the
underlying data. However, further investigation is required to answer various
related questions in this regard. We state the problem more formally next.
Problem 1. Let the time-series data be denoted as the sequence
$\\{X_{t}\\}_{t\in\mathbb{N}}$ where $X_{t}\in\Omega\subseteq\mathbb{R}^{d}$
and let $\varphi$ represent the partitioning function such that
$\varphi:X_{t}\mapsto a_{t}$ where $a_{t}\in\mathcal{A}$ for all
$t\in\mathbb{N}$ and $\mid\mathcal{A}\mid\in\mathbb{N}$ is known and fixed.
Then the partitioning function is completely determined by the set
$\mathcal{R}_{\varphi}=\\{R_{0},R_{1},\dots,R_{|\mathcal{A}|}\\}$ such that
$\Omega=\overline{R_{0}}\cup\overline{R_{1}}\cup\dots\cup\overline{R_{|\mathcal{A}|}}$.
Let us assume that the possible family of sets lies in a set denoted by
$\Omega_{\varphi}$. An information gain for the Markov model could be measured
by the following equation.

(20) $d_{\mathcal{R}_{\varphi}}=\sum_{q\in
Q}\text{Pr}(q)D_{KL}(P(\mathcal{A}\mid q)\|\tilde{P}({\mathcal{A}}))$

where the measure is parameterized by the set $\mathcal{R}_{\varphi}$ which
depends on the partitioning function $\varphi$ and the set $Q$ represents the
finite memory-words of the discrete symbol sequence. The above equation
measures the information gain by creating a Markov model, where we measure the
symbol emission probabilities conditioned on the memory words (or states) of
the Markov model, over the independent or marginal symbol emission
probabilities. Another interpretation is that equation (20) describes the
discrepancy between the statistics of $\\{s_{t}\\}_{t\in\mathbb{N}}$ when
modeled as independent sequence versus when modeled as a stationary Markov
process. Then, a partitioning to optimize this measure may capture the true
temporal behavior of $\\{X_{t}\\}_{t\in\mathbb{N}}$. The problem is to obtain
the parameters of the related optimization problem.

$\mathcal{R}_{\varphi}^{\star}=\arg\max\limits_{\mathcal{R}_{\varphi}\in\Omega_{\varphi}}\sum_{q\in
Q}\text{Pr}(q)D_{KL}(P(\mathcal{A}\mid q)\|\tilde{P}({\mathcal{A}}))$

The following questions need to be answered to characterize
$\mathcal{R}_{\varphi}^{\star}$.

  * •

Is $\mathcal{R}_{\varphi}^{\star}$ unique? Under what conditions of the
underlying process $\\{X\\}_{t\in\mathbb{N}}$, can we get an unique solution?

  * •

What is the order of the corresponding discrete time-series obtained for
$\mathcal{R}_{\varphi}^{\star}$?

  * •

Let us assume that the pre-image of a symbol is represented by the centroid of
the set in $\mathcal{R}_{\varphi}^{\star}$ in the original phase-space. Then,
how can we characterize the metric
$\zeta_{\varphi}=\sum_{t\in\mathbb{N}}d({\varphi}^{-1}(s_{t}),X_{t})$ for
signal representation? Is $\mathcal{R}_{\varphi}^{\star}$ able to minimize
this metric over $\Omega_{\varphi}$?

  * •

Now imagine that the size of partitioning set is allowed to vary. Then, how to
assign a MDL score to the individual models for different sizes of the
partitioning set and how do we select a final model for signal representation?

The above problem would, thus, try to formulate and characterize the
properties of a partition for a data set for Markov representations. In the
next question, we will try to study the composite problem of signal
representation by partitioning followed by order estimation.
Problem 2. Let the time-series data be denoted as the sequence
$\\{X_{t}\\}_{t\in\mathbb{N}}$ where $X_{t}\in\Omega\subseteq\mathbb{R}^{d}$
and let $\varphi$ represent the partitioning function such that
$\varphi:X_{t}\mapsto a_{t}$ where $a_{t}\in\mathcal{A}$ for all
$t\in\mathbb{N}$ and $\mid\mathcal{A}\mid\in\mathbb{N}$ is unknown. A
desirable way to characterize a partitioning process is by predicting its
effect on the size of temporal memory of the system. Is it possible to
synthesize a partitioning function $\varphi_{I}$ which preserves the memory of
the discrete system under the transformation $\varphi_{I}$, i.e., if
$\text{Pr}(X_{t}\mid X_{t-D},\dots,X_{t-1})=\text{Pr}(X_{t}\mid
X_{t-1},\dots)$ then we have
$\text{Pr}(\varphi_{I}(X_{t})\mid\varphi_{I}(X_{t-D}),\dots,\varphi_{I}(X_{t-1}))=\text{Pr}(\varphi_{I}(X_{t})\mid\varphi_{I}(X_{t-1}),\dots)$.
Then, how can we characterize a system for which such a discretization is
guaranteed to exist? It is noted that the size of partitioning set is unknown
and not fixed.

These two problems could be treated as fundamental problems that need to be
studied for mathematical characterization of data-driven modeling of systems
from a symbolic analysis perspective. They are mainly concerned about
inference of model structure for statistical learning. There are some problems
which are of interest from applications perspective and related to estimation
of various parameters during modeling and inference. However, they are not
being presented here.

## References

  * [1] G. E. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung, Time series analysis: forecasting and control.  John Wiley & Sons, 2015.
  * [2] C. M. Bishop, Pattern recognition and machine learning.  Springer, 2006.
  * [3] C. S. Daw, C. E. A. Finney, and E. R. Tracy, “A review of symbolic analysis of experimental data,” Review of Scientific Instruments, vol. 74, no. 2, pp. 915–930, 2003.
  * [4] A. Ray, “Symbolic dynamic analysis of complex systems for anomaly detection,” Signal Processing, vol. 84, pp. 1115–1130, July 2004.
  * [5] K. Mukherjee and A. Ray, “State splitting and merging in probabilistic finite state automata for signal representation and analysis,” Signal Processing, vol. 104, pp. 105–119, 2014.
  * [6] Y. Hirata, K. Judd, and D. Kilminster, “Estimating a generating partition from observed time series: Symbolic shadowing,” Physical Review E, vol. 70, no. 1, p. 016215, 2004.
  * [7] M. Buhl and M. Kennel, “Statistically relaxing to generating partitions for observed time-series data,” Physical Review E, vol. 71, no. 4, p. 046213, 2005.
  * [8] R. Adler, “Symbolic dynamics and Markov partitions,” Bulletin of the American Mathematical Society, vol. 35, no. 1, pp. 1–56, 1998.
  * [9] M. B. Kennel and M. Buhl, “Estimating good discrete partitions from observed data: Symbolic false nearest neighbors,” Physical Review Letters, vol. 91, no. 8, p. 084102, 2003.
  * [10] V. Rajagopalan and A. Ray, “Symbolic time series analysis via wavelet-based partitioning,” Signal Processing, vol. 86, pp. 3309–3320, November 2006\.
  * [11] J. Lin, E. Keogh, L. Wei, and S. Lonardi, “Experiencing SAX: a novel symbolic representation of time series,” Data Mining and Knowledge Discovery, vol. 15, pp. 107–144, October 2007.
  * [12] S. Sarkar, A. Srivastav, and M. Shashanka, “Maximally bijective discretization for data-driven modeling of complex systems,” American Control Conference, Washington, D.C., USA, pp. 2680–2685, June 2013.
  * [13] Y. Yang and G. I. Webb, “A comparative study of discretization methods for naive-bayes classifiers,” in Proceedings of PKAW, vol. 2002, 2002.
  * [14] F. Fleuret, “Fast binary feature selection with conditional mutual information,” The Journal of Machine Learning Research, vol. 5, pp. 1531–1555, 2004.
  * [15] H. Liu, F. Hussain, C. Tan, and M. Dash, “Discretization: An enabling technique,” Data Mining and Knowledge Discovery, vol. 6, pp. 393–423, 2002\.
  * [16] S. Garcia, J. Luengo, J. A. Saez, V. Lopez, and F. Herrera, “A survey of discretization techniques: Taxonomy and empirical analysis in supervised learning,” IEEE Transactions on Knowledge and Data Engineering, vol. 99, no. PrePrints, 2012.
  * [17] S. Sarkar, P. Chattopdhyay, A. Ray, S. Phoha, and M. Levi, “Alphabet size selection for symbolization of dynamic data-driven systems: An information-theoretic approach,” in American Control Conference (ACC), 2015, pp. 5194–5199, IEEE, 2015.
  * [18] R. Steuer, L. Molgedey, W. Ebeling, and M. A. Jimenez-Montaño, “Entropy and optimal partition for data analysis,” The European Physical Journal B-Condensed Matter and Complex Systems, vol. 19, no. 2, pp. 265–269, 2001.
  * [19] J. Dougherty, R. Kohavi, and M. Sahami, “Supervised and unsupervised discretization of continuous features,” in Machine learning: proceedings of the twelfth international conference, vol. 12, pp. 194–202, 1995\.
  * [20] M. J. Pazzani, “An iterative improvement approach for the discretization of numeric attributes in bayesian classifiers.,” in KDD, pp. 228–233, 1995\.
  * [21] A. Barron, J. Rissanen, and B. Yu, “The minimum description length principle in coding and modeling,” Information Theory, IEEE Transactions on, vol. 44, no. 6, pp. 2743–2760, 1998.
  * [22] I. Csiszár and Z. Talata, “Context tree estimation for not necessarily finite memory processes, via BIC and MDL,” Information Theory, IEEE Transactions on, vol. 52, no. 3, pp. 1007–1016, 2006.
  * [23] T. M. Cover and J. A. Thomas, Elements of information theory.  John Wiley & Sons, 2012.
  * [24] H. Tong, “Determination of the order of a Markov chain by Akaike’s information criterion,” Journal of Applied Probability, pp. 488–497, 1975\.
  * [25] R. W. Katz, “On some criteria for estimating the order of a Markov chain,” Technometrics, vol. 23, no. 3, pp. 243–249, 1981.
  * [26] I. Csiszár and P. C. Shields, “The consistency of the BIC Markov order estimator,” The Annals of Statistics, vol. 28, no. 6, pp. 1601–1619, 2000\.
  * [27] I. Csiszár, “Large-scale typicality of markov sample paths and consistency of MDL order estimators,” Information Theory, IEEE Transactions on, vol. 48, no. 6, pp. 1616–1628, 2002.
  * [28] G. Morvai and B. Weiss, “On estimating the memory for finitarily Markovian processes,” in Annales de l’Institut Henri Poincare (B) Probability and Statistics, vol. 43, pp. 15–30, Elsevier, 2007.
  * [29] G. Morvai and B. Weiss, “On estimating the memory for finitarily markovian processes,” in Annales de l’Institut Henri Poincare (B) Probability and Statistics, vol. 43, pp. 15–30, Elsevier, 2007.
  * [30] G. Morvai and B. Weiss, “Estimating the lengths of memory words,” Information Theory, IEEE Transactions on, vol. 54, no. 8, pp. 3804–3807, 2008\.
  * [31] G. Morvai and B. Weiss, “Universal tests for memory words,” Information Theory, IEEE Transactions on, vol. 59, no. 10, pp. 6873–6879, 2013.
  * [32] M. Papapetrou and D. Kugiumtzis, “Markov chain order estimation with conditional mutual information,” Physica A: Statistical Mechanics and its Applications, vol. 392, no. 7, pp. 1593–1601, 2013.
  * [33] L. Zhao, C. Dorea, and C. Gonçalves, “On determination of the order of a Markov chain,” Statistical inference for stochastic processes, vol. 4, no. 3, pp. 273–282, 2001.
  * [34] N. Merhav, M. Gutman, and J. Ziv, “On the estimation of the order of a Markov chain and universal data compression,” Information Theory, IEEE Transactions on, vol. 35, pp. 1014–1019, Sep 1989.
  * [35] V. Rajagopalan, A. Ray, R. Samsi, and J. Mayer, “Pattern identification in dynamical systems via symbolic time series analysis,” Pattern Recogn., vol. 40, no. 11, pp. 2897–2907, 2007.
  * [36] A. Srivastav, “Estimating the size of temporal memory for symbolic analysis of time-series data,” American Control Conference, Portland, OR, USA, pp. 1126–1131, June 2014.
  * [37] D. K. Jha, A. Srivastav, K. Mukherjee, and A. Ray, “Depth estimation in Markov models of time-series data via spectral analysis,” in American Control Conference (ACC), 2015, pp. 5812–5817, IEEE, 2015.
  * [38] S. Chakraborty, A. Ray, A. Subbu, and E. Keller, “Analytic signal space partitioning and symbolic dynamic filtering for degradation monitoring of electric motors,” Signal, Image and Video Processing, vol. 4, no. 4, pp. 399–403, 2010.
  * [39] E. Vidal, F. Thollard, C. de la Higuera, F. Casacuberta, and R. Carrasco, “Probabilistic finite-state machines- Part I,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 27, no. 7, pp. 1013–1025, 2005.
  * [40] E. Vidal, F. Thollard, C. de la Higuera, F. Casacuberta, and R. Carrasco, “Probabilistic finite-state machines- Part II,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 27, no. 7, pp. 1026–1039, 2005.
  * [41] H. L. Van Trees, Detection, estimation, and modulation theory.  John Wiley & Sons, 2004.
  * [42] H. V. Poor, An introduction to signal detection and estimation.  Springer Science & Business Media, 2013.
  * [43] B. Laurence, “A complete sufficient statistic for finite-state markov processes with application to source coding,” IEEE transactions on information theory, vol. 39, no. 3, p. 1047, 1993.
  * [44] D. Lind and B. Marcus, An Introduction to Symbolic Dynamics and Coding.  Cambridge University Press, 1995.
  * [45] P. Grunwald, “A tutorial introduction to the minimum description length principle,” arXiv preprint math/0406077, 2004.
  * [46] V. N. Vapnik, Statistical learning theory, vol. 1.  Wiley New York, 1998.
  * [47] V. Vapnik, The nature of statistical learning theory.  Springer Science & Business Media, 2013.
  * [48] D. K. Jha, N. Virani, J. Reimann, A. Srivastav, and A. Ray, “Symbolic analysis-based reduced order markov modeling of time series data,” Signal Processing, vol. 149, pp. 68–81, 2018.
  * [49] Y. Li, D. K. Jha, A. Ray, and T. A. Wettergren, “Information fusion of passive sensors for detection of moving targets in dynamic environments,” IEEE transactions on cybernetics, vol. 47, no. 1, pp. 93–104, 2016.
  * [50] S. Sarkar, D. K. Jha, A. Ray, and Y. Li, “Dynamic data-driven symbolic causal modeling for battery performance & health monitoring,” in 2015 18th International Conference on Information Fusion (Fusion), pp. 1395–1402, IEEE, 2015.
  * [51] Y. Seto, N. Takahashi, D. K. Jha, N. Virani, and A. Ray, “Data-driven robot gait modeling via symbolic time series analysis,” in 2016 American Control Conference (ACC), pp. 3904–3909, IEEE, 2016.
  * [52] N. Virani, D. K. Jha, A. Ray, and S. Phoha, “Sequential hypothesis tests for streaming data via symbolic time-series analysis,” Engineering Applications of Artificial Intelligence, vol. 81, pp. 234–246, 2019.
  * [53] D. K. Jha, A. Ray, K. Mukherjee, and S. Chakraborty, “Classification of two-phase flow patterns by ultrasonic sensing,” Journal of Dynamic Systems, Measurement, and Control, vol. 135, no. 2, 2013.
  * [54] N. Virani, S. Marcks, S. Sarkar, K. Mukherjee, A. Ray, and S. Phoha, “Dynamic data driven sensor array fusion for target detection and classification,” Procedia Computer Science, vol. 18, pp. 2046–2055, 2013.
  * [55] Y. Li, D. K. Jha, A. Ray, and T. A. Wettergren, “Information-theoretic performance analysis of sensor networks via markov modeling of time series data,” IEEE transactions on cybernetics, vol. 48, no. 6, pp. 1898–1909, 2017.
  * [56] Y. Li, D. K. Jha, A. Ray, and T. A. Wettergren, “Feature level sensor fusion for target detection in dynamic environments,” in 2015 American Control Conference (ACC), pp. 2433–2438, IEEE, 2015.
  * [57] J. O’Connor, V. Acharya, and T. Lieuwen, “Transverse combustion instabilities: Acoustic, fluid mechanic, and flame processes,” Progress in Energy and Combustion Science, vol. 49, pp. 1–39, 2015.
  * [58] Sé, b. Ducruix, T. Schuller, D. Durox, Sé, and b. Candel, “Combustion dynamics and instabilities: Elementary coupling and driving mechanisms,” Journal of Propulsion and Power, vol. 19, no. 5, pp. 722–734, 2003.
  * [59] S. Candel, D. Durox, T. Schuller, J.-F. Bourgouin, and J. P. Moeck, “Dynamics of swirling flames,” Annual review of fluid mechanics, vol. 46, pp. 147–173, 2014.
  * [60] Y. Huang and V. Yang, “Dynamics and stability of lean-premixed swirl-stabilized combustion,” Progress in Energy and Combustion Science, vol. 35, no. 4, pp. 293–364, 2009.
  * [61] J. P. Moeck, J.-F. Bourgouin, D. Durox, T. Schuller, and S. Candel, “Nonlinear interaction between a precessing vortex core and acoustic oscillations in a turbulent swirling flame,” Combustion and Flame, vol. 159, no. 8, pp. 2650–2668, 2012.
  * [62] A. Banaszuk, P. G. Mehta, C. A. Jacobson, and A. I. Khibnik, “Limits of achievable performance of controlled combustion processes,” Control Systems Technology, IEEE Transactions on, vol. 14, no. 5, pp. 881–895, 2006\.
  * [63] A. Banaszuk, P. G. Mehta, and G. Hagen, “The role of control in design: From fixing problems to the design of dynamics,” Control Engineering Practice, vol. 15, no. 10, pp. 1292–1305, 2007.