Bayesian Nonparametric Modeling of Categorical Data for Information Fusion and Causal Inference ^†

Abstract

This paper presents a nonparametric regression model of categorical time series in the setting of conditional tensor factorization and Bayes network. The underlying algorithms are developed to provide a flexible and parsimonious representation for fusion of correlated information from heterogeneous sources, which can be used to improve the performance of prediction tasks and infer the causal relationship between key variables. The proposed method is first illustrated by numerical simulation and then validated with two real-world datasets: (1) experimental data, collected from a swirl-stabilized lean-premixed laboratory-scale combustor, for detection of thermoacoustic instabilities and (2) publicly available economics data for causal inference-making.

1. Introduction

Modeling and decision-making in complex dynamical systems (e.g., distributed physical processes [1], macro-economy [2] and human brain [3]) often rely on time series collected from heterogeneous sources. Fusion of the information extracted from an ensemble of time series is a critical ingredient for better prediction and causal inference.

In many dynamical systems, the characteristic time of the physical process under consideration is small (e.g., around 2 ms in a typical combustion process) relative to the time-scale of respective decision-making (e.g., tenths of a second for active combustion control). Therefore, fast and accurate prediction of the system states and estimation of the associated parameters is essential for online monitoring and active control of the dynamical system; for example, real-time prediction of future states can significantly improve active control of thermoacoustic instabilities [4]. One way to achieve this is to make predictions based on different but correlated information sources. Although several methods have been proposed for prediction based on fusion of heterogeneous time series (e.g., [5,6,7]), they lack a coherent probabilistic interpretation and may not be able to accommodate more general interactions between current measurements and the measurement history. Furthermore, these methods may not be sequentially implementable and hence they may not be very useful for real-time applications.

Identification of causal relationships is essential for understanding the consequences of transitions from empirical findings to actions and thus forms a significant part of knowledge discovery. Various analytical techniques (e.g., [8,9,10]) have been proposed for causal inference-making; among these techniques, the concept of causality introduced by Granger [11], hereafter called Granger causality, is apparently one of the most widely used in time series analysis [12]. Granger causality does not rely on the specification of a scientific model and thus is particularly applicable to investigation of empirical cause-effect relationships. It is noted that Granger causality is especially suited for continuous-valued data based on frequentist hypothesis testing.

The goal of this paper is to develop a flexible and parsimonious model of categorical time series in a Bayesian nonparametric setting for fusion of correlated information from heterogeneous sources (e.g., sensors of possibly different modalities), which can be used for sequential classification and causal inference. From this perspective, major contributions of the paper are delineated as follows:

• By introducing latent variables and sparsity inducing priors, a flexible and parsimonious model is developed for fusion of correlated information from heterogeneous sources (e.g., sensors of possibly different modalities), which can be used to improve the performance of sequential classification tasks.
• By testing the dimension of latent variables in the setting of Bayes factor analysis [13], Granger causality [11] is extended to categorical time series.
• Validation of the above concept with experimental data, generated from a swirl-stabilized lean-premixed laboratory-scale combustor [14], for real-time detection of thermoacoustic instabilities.
• Testing of the underlying algorithm with public economics data to infer the causal relationship between two categorical time series.

The paper is organized into eight sections including the current section. Section 2 introduces the concept of Granger causality and develops the model. Section 3 discusses the algorithm for posterior computation using Gibbs sampling, and hypothesis testing using Bayes factor analysis. Section 4 presents the sequential classification algorithm with the proposed model. The underlying algorithms are tested with simulation data in Section 5 while Section 6 validates the proposed method with some experimental data, collected from a swirl-stabilized lean-premixed laboratory-scale combustor, for thermoacoustic instabilities early detection. Section 7 validates the proposed concept on publicly available economics data. Section 8 concludes the paper and provides a few recommendations for future research. The nomenclature and list of acronyms are provided at the end before the list of references.

2. Model Development

This section first introduces the concept of Granger causality and the corresponding regression model. Next, the underlying model’s algebraic and statistical specifications are elaborated.

Definition 1.

(Granger Causality) Let ${\left\{y_t\right\}}_{t=1}^T$ and ${\left\{{\theta }_t\right\}}_{t=1}^T$ be two (statistically) stationary categorical time series. Then, the variable θ Granger-causes the variable y if the past values of θ contain statistically significant information for predictions of y besides those contained in the past values of y. Similarly, y Granger-causes θ if the past values of y contain statistically significant information for predictions of θ besides those contained in the past values of θ.

Remark 1.

The following are four types of Granger causality relationship between θ and y:

1.

θ Granger-causes y but not the vice versa;

2.

y Granger-causes θ but not the vice versa;

3.

θ and y Granger-cause each other;

4.

θ does not Granger-cause y and vice versa.

However, in practice, only finitely many past values of y and $\theta$ are considered. To test the null hypothesis that $\theta$ does not Granger-cause y, the following regression model is constructed:

$$p(y_t\mid y_{t-1},\cdots ,y_{t-D_y},{\theta }_{t-1},\cdots ,{\theta }_{t-D_{\theta }})$$

(1)

where in this model, predictors $y_{t-1}$ to $y_{t-D_y}$ represent variable y’s time lags; and predictors ${\theta }_{t-1}$ to ${\theta }_{t-D_{\theta }}$ represent variable $\theta$’s time lags. In the sequel, for simplicity of notations, predictors ${\mathit{z}}_t\equiv (z_{1,t},\cdots ,z_{q,t})$ are substituted for $(y_{t-1},\cdots ,y_{t-D_y},{\theta }_{t-1},\cdots ,{\theta }_{t-D_{\theta }})$.

Remark 2.

If the explanatory power of ${\theta }_{t-1},\cdots ,{\theta }_{t-D_{\theta }}$ to the regression is significant, then the null hypothesis (that θ does not Granger-cause y) is rejected and the alternative hypothesis (that θ Granger-causes y) is accepted. Hypothesis tests on the significance of time-lags are elaborated later in Equation (15) (see Section 3.2).

Remark 3.

If y and θ are correlated in the sense of Granger causality, the information contained in one source can be used to predict the future values in another source. Accordingly, information fusion of different sources enables fast and accurate prediction because of Granger-causality. It is noted that if the information contained in two sources is statistically independent, then information fusion cannot enhance prediction accuracy.

2.1. Conditional Tensor Factorization

This subsection addresses fusion of different sources of information by making use of the concept of conditional probability tensor that was first reported in [15], a formal definition of conditional probability tensor follows.

Definition 2.

(Conditional probability tensor) Let $C_0$ denote the number of categories of the (one-dimensional) variable, $y_t$, and let $C_j$ denote the number of categories of $z_{j,t}$ for $j=1,\cdots ,q$, where is the number of predictors. The quantity $p(y_t\mid {\mathit{z}}_t)$ is treated as a $(q+1)th$ order tensor in the $C_0\times C_1\cdots \times C_q$ dimensional space, hereafter called the conditional probability tensor.

Let $C_y$ and $C_{\theta }$ respectively denote the numbers of categories of the variables, y and $\theta$. It follows from Definition 2 that $C_1=\cdots =C_{D_y}=C_y$ and $C_{D+1}=\cdots =C_q=C_{\theta }$. Then, each one of these conditional probability tensors has a higher order singular value decomposition (HOSVD) of the following form [15]:

$$p(y_t\mid {\mathit{z}}_t)=\sum _{s_1=1}^{k_1}\cdots \sum _{s_q=1}^{k_q}{\lambda }_{s_1,\cdots ,s_q}\left(y_t\right)\prod _{j=1}^q{\omega }_{s_j}^{\left(j\right)}\left(z_{j,t}\right)$$

(2)

where $1\le k_j\le C_j$ for $j=1,\cdots ,q$; and each of the parameters ${\lambda }_{s_1\cdots s_q}\left(y_t\right)$ and ${\omega }_{s_j}^{\left(j\right)}\left(z_{j,t}\right)$ is non-negative while the following constraints are satisfied:

$$\begin{array}{cc}\hfill \sum _{y_t=1}^{C_0}{\lambda }_{s_1,\cdots ,s_q}\left(y_t\right)& =1,\mathrm{for}\mathrm{each}(s_1,\cdots ,s_q)\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \sum _{s_j=1}^{k_j}{\omega }_{s_j}^{\left(j\right)}\left(z_{j,t}\right)& =1,\mathrm{for}\mathrm{each}(j,z_{j,t})\hfill \end{array}$$

(4)

Remark 4.

Since there exists a factorization as in Equation (2) for each one of the conditional probability tensors, the two constraints Equations (3) and (4) are not restrictive. Furthermore, it is ensured that ${\sum }_{y_t=1}^{C_0}p(y_t\mid {\mathit{z}}_t)=1$.

2.2. Bayesian Nonparametric Modeling

In order to build a statistically interpretable model, two techniques can be used to convert the tensor factorization in Equation (2) to a Bayes network, i.e., (1) introduce latent allocation-class variables; (2) assign sparsity-inducing priors. To this end, T pairs of variables and their respective predictors are collected in one dataset, and it is rearranged as ${\{y_t,{\mathit{z}}_t\}}_{t=1}^T$, where t is an index with range from 1 to T.

The conditional probability $p(y_t\mid {\mathit{z}}_t)$, factorized as in Equation (2), is then reorganized in the following form:

$$p(y_t\mid {\mathit{z}}_t)={\int }_{x_{1,t}}\cdots {\int }_{x_{q,t}}p(y_t\mid {\mathit{x}}_t)\prod _{j=1}^{q}p(x_{j,t}\mid z_{j,t})$$

(5)

where ${\mathit{x}}_t\equiv (x_{1,t},\cdots ,x_{q,t})$ denotes the latent class-allocation variables.

For index $j=1,\cdots ,q$ and index $t=1,\cdots ,T$, it then follows that

$$\begin{array}{cc}\hfill x_{j,t}\mid {\mathit{\omega }}^{\left(j\right)},z_{j,t}& \sim \mathbf{Mult}({\mathit{\omega }}^{\left(j\right)}\left(z_{j,t}\right))\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill y_t\mid \tilde{\mathit{\lambda }},{\mathit{x}}_t& \sim \mathbf{Mult}({\tilde{\mathit{\lambda }}}_{{\mathit{x}}_t})\hfill \end{array}$$

(7)

where $\mathbf{Mult}(•)$ is the multinomial distribution [16] and ${\mathit{\omega }}^{\left(j\right)}\equiv {\{{\{{\omega }_s^{\left(j\right)}\left(c\right)\}}_{s=1}^{k_j}\}}_{c=1}^{C_j}$ is the mixture probability matrix. The cth row ${\mathit{\omega }}^{\left(j\right)}\left(c\right)\equiv {\{{\omega }_s^{\left(j\right)}\left(c\right)\}}_{s=1}^{k_j}$ in this mixture probability matrix is a probability vector itself (i.e., it sums to 1). Moreover, $\tilde{\mathit{\lambda }}\equiv {\{{\mathit{\lambda }}_{s_1,\cdots ,s_q}\}}_{(s_1\cdots s_q)}$ is a conditional probability tensor where ${\mathit{\lambda }}_{s_1,\cdots ,s_q}\equiv {\{{\lambda }_{s_1,\cdots ,s_q}\left(c\right)\}}_{c=1}^{C_0}$ is a probability vector for each string $(s_1,\cdots ,s_q)$.

The hierarchical reformulation of HOSVD above illustrates the following features of this model in Equation (5):

• Soft clustering for each one of the predictors ${\mathit{z}}_j\equiv {\{z_{j,t}\}}_{t=1}^T$ is implemented following Equation (6). This allows for inheritance of statistical strengths across different categories.
• The distribution of variable $y_t$ is determined by a probability tensor $\tilde{\lambda }$ of reduced order, following Equation (7).
• In order to capture the interactions among different predictors, class assignment variables ${\mathit{x}}_j\equiv {\{x_{j,t}\}}_{t=1}^T$ are used. They work in an implicit and parsimonious way by allowing the latent populations with the index of $(s_1,\cdots ,s_q)$ to be shared across various state combinations of predictors.

Remark 5.

Here it is very critical to distinguish these two different concepts: (1) the number of clusters ${\tilde{k}}_j$ generated by the latent class variables ${\mathit{x}}_j$ and (2) the dimensions $k_j$ of the probability vector ${\mathit{\omega }}^{\left(j\right)}\left(c\right)$ from the mixture probability matrix. The former one represents the number of groups generated by the data, and is smaller than the latter. It should be noted that ${\tilde{k}}_j$ determines if the predictor ${\mathit{z}}_j$ should be included in the model, because $p(y_t\mid {\mathit{z}}_t)$ will not change with $z_{j,t}$ if ${\mathit{z}}_j$ has just a single latent cluster. Thus the significance of some particular predictor could be tested using on ${\tilde{k}}_j$, which is later elaborated in Section 3.2.

In many real-world applications, the tensor $\tilde{\mathit{\lambda }}$ often has more components than needed, since the product ${\prod }_{j=1}^q{k}_j$ can be large even for modest values of q and $C_j$. To deal with this problem, tensor $\tilde{\mathit{\lambda }}$ is then clustered within different combinations of $(s_1,\cdots ,s_q)$ nonparametrically by imposing a Pitman-Yor process prior [17]. Then, by using the stick-breaking representation of the Pitman-Yor process [18], it follows that

$$\begin{array}{cc}\hfill {\mathit{\lambda }}_l\mid \gamma & \sim \mathbf{Dir}\left(\alpha \right),\mathrm{for}l=1,\cdots ,\infty \hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill V_k\mid a,b& \sim \mathbf{Beta}(1-b,a+kb),\mathrm{for}k=1,\cdots ,\infty \hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\pi }_l& =V_l\prod _{k=1}^{l-1}(1-V_k),\mathrm{for}l=1,\cdots ,\infty \hfill \end{array}$$

(10)

where the bold symbols $\mathbf{Dir}(•)$ and $\mathbf{Beta}(•)$ represents the uniform Dirichlet distributions and Beta distributions [16] respectively, and ${\mathit{\lambda }}_l\equiv ({\lambda }_l\left(1\right),\cdots ,{\lambda }_l\left(C_0\right))$. Moreover, $0\le b<1$ and $a>-b$. For each combination $(s_1,\cdots ,s_q)$, it follows that

$${\varphi }_{s_1,\cdots ,s_q}\mid \mathit{\pi }\sim \mathbf{Mult}\left(\pi \right)$$

(11)

where $\mathit{\pi }\equiv ({\pi }_1,{\pi }_2,\cdots )$. For $t=1,\cdots ,T$,

$$y_t\mid \mathit{\lambda },\mathit{\varphi },{\mathit{x}}_t\sim \mathbf{Mult}\left({\mathit{\lambda }}_{{\varphi }_{{\mathit{x}}_t}}\right)$$

(12)

where $\mathit{\lambda }\equiv {\{{\mathit{\lambda }}_l\}}_{l=1}^{\infty }$ and $\mathit{\varphi }\equiv {\{{\varphi }_{s_1,\cdots ,s_q}\}}_{(s_1,\cdots ,s_q)}$.

The next step assigns priors to the mixture probability matrix ${\mathit{\omega }}^{\left(j\right)}$. Here the dimension of ${\mathit{\omega }}^{\left(j\right)}$ grows linearly as $k_j$ increases (unlike the tensor $\tilde{\mathit{\lambda }}$). Therefore, further clustering of ${\mathit{\omega }}^{\left(j\right)}$ is not necessary. Hence, we assign independent priors to the rows of ${\mathit{\omega }}^{\left(j\right)}$ for $j=1,\cdots ,q$ in the following way:

$${\mathit{\omega }}^{\left(j\right)}\left(c\right)\mid k_j,{\beta }_j\sim \mathbf{Dir}\left({\beta }_j\right),\mathrm{for}c=1,\cdots ,C_j$$

(13)

Lastly, we assign priors to the dimension of the mixture probability vector $k_j$, i.e., for $j=1,\cdots ,q$,

$$p(k_j=k\mid {\mu }_j)\propto \exp (-{\mu }_{j}k),\mathrm{for}k=1,\cdots ,C_j$$

(16)

where ${\mu }_j\ge 0$ and $\mathbf{k}\equiv {\{k_j\}}_{j=1}^q$.

Remark 6.

As the parameter ${\mu }_j$ grows larger, the exponential prior in Equation (14) will assign increasing probabilities to smaller values of $k_j$, and it becomes a uniform prior distribution on $\{1,\cdots ,C_j\}$ when ${\mu }_j$ is zero. Commonly, people have prior beliefs that as time lags increase, they will a have vanishing impact on the distribution of the current response variable. To impose this prior belief, we can assign larger ${\mu }_j$ to time lags further back in the history.

By combining Equations (6)–(14) together, a Bayes network representation of the model is created and Figure 1 illustrates its structures.

3. Estimation and Inference

  This section presents the details of an algorithm for computing posteriors as well as Bayesian hypothesis testing by using Bayes factors.

3.1. Posterior Computation

Despite the fact that the posterior distribution does not have any specific analytical form, we can still perform the inference of the corresponding Bayes network by using Gibbs sampling method. Because the dimension of ${\mathit{\omega }}^{\left(j\right)}$ may vary with $k_j$, constructing a stationary Markov chain by plain Gibbs sampling is difficult. To infer a model with variable dimensions, a common analytical tool—the reversible jump Monte Carlo Markov chain (MCMC) [19], which does trans-dimensional exploration in the model space—is often used.

Product partition modeling [20,21] can help alleviate difficulties occurring in trans-dimensional modeling by constructing a stationary Markov chain on the clustering space. For this proposed method, the dimension ${\mathit{\omega }}^{\left(j\right)}$ is being integrated out for the sampling of $k_j$ directly from $p(k_j\mid {\mathit{x}}_j,{\mathit{z}}_j)$, which will create a partially collapsed Gibbs sampler [22] that alternates between these two spaces: (1) the space with all the variables and (2) the space with all the variables but $\mathit{\omega }={\{{\mathit{\omega }}^{\left(j\right)}\}}_{j=1}^q$.

To compute the posterior probabilities of the Pitman-Yor process, the infinite-dimensional tensors $\mathit{\pi }$ and $\mathit{\lambda }$ after their Lth component are truncated, as performed in [18]. For achieving desired accuracy, an appropriate L needs to be chosen. Other than this, the posterior sampling is rather straightforward. The detailed process is presented in Algorithm 1, in which it is not explicitly mentioned that $\mathit{x}\equiv {\{{\mathit{x}}_t\}}_{t=1}^T$ and $\mathit{\xi }$ collects the variables.

Algorithm 1 Gibbs sampling for the proposed method

Input: Datasets ${\{y_t,{\mathit{z}}_t\}}_{t=1}^T$; hyperparameters a, b, $\alpha$, ${\{{\mu }_j\}}_{j=1}^q$, ${\{{\beta }_j\}}_{j=1}^q$; number of truncating components L; number of all samples N; initial sample ${}_{\left(0\right)}\mathit{\varphi },{}_{\left(0\right)}\mathit{\pi },{}_{\left(0\right)}\mathit{\lambda },{}_{\left(0\right)}\mathit{\omega },{}_{\left(0\right)}\mathit{x},{}_{\left(0\right)}\mathit{k}$. Output: All posterior samples ${\{{}_{\left(n\right)}\mathit{\varphi },{}_{\left(n\right)}\mathit{\pi },{}_{\left(n\right)}\mathit{\lambda },{}_{\left(n\right)}\mathit{\omega },{}_{\left(n\right)}\mathit{x},{}_{\left(n\right)}\mathit{k}\}}_{n=1}^N$ 1: for$n=1$ to N do 2:    For each one string $(s_1,\cdots ,s_q)$, collect a sample ${\varphi }_{s_1,\cdots ,s_q}$ from its multinomial full conditional $$p({\varphi }_{s_1,\cdots ,s_q}=l\mid \mathit{\xi })\propto {\pi }_l\prod _{c=1}^{C_0}{\{{\lambda }_l\left(c\right)\}}^{n_{s_1,\cdots ,s_q}\left(c\right)}$$    where $n_{s_1,\cdots ,s_q}\left(c\right)={\sum }_1^T\mathbf{1}\{x_{1,t}=s_1,\cdots ,x_{q,t}=s_q,y_t=c\}$. 3:    For $l=1,\cdots ,L$, update ${\pi }_l$ by the following rules $$\begin{array}{c}{V}_l\mid \mathit{\xi }\sim \mathbf{Beta}(1-b+n_l,a+lb+\sum _{k>l}{n}_k),l<L\\ V_L=1,{\pi }_l=V_l\prod _{k=1}^{l-1}(1-V_k)\end{array}$$    where $n_l={\sum }_{(s_1,\cdots ,s_q)}\mathbf{1}\{{\varphi }_{s_1,\cdots ,s_q}=l\}$. 4:    For $l=1,\cdots ,L$, collect samples ${\mathit{\lambda }}_l$ from their respective Dirichlet full conditionals $${\mathit{\lambda }}_l\mid \mathit{\xi }\sim \mathbf{Dir}\{\alpha +n_l\left(1\right),\cdots ,\alpha +n_l\left(C_0\right)\}$$    where $n_l\left(c\right)={\sum }_{(s_1,\cdots ,s_q)}\mathbf{1}\{{\varphi }_{s_1,\cdots ,s_q}=l\}{n}_{s_1,\cdots ,s_q}\left(c\right)$. 5:    For $j=1,\cdots ,q$, for $c=1,\cdots ,C_j$, collect samples $${\mathit{\omega }}^{\left(j\right)}\left(c\right)\mid \mathit{\xi }\sim \mathbf{Dir}\{{\beta }_j+n_{j,c}\left(1\right),\cdots ,{\beta }_j+n_{j,c}\left(k_j\right)\}$$    where $n_{j,c}\left(s_j\right)={\sum }_{t=1}^T\mathbf{1}\{x_{j,t}=s_j,z_{j,t}=c\}$. 6:    For $j=1,\cdots ,q$, for $t=1,\cdots ,T$, collect samples $x_{j,t}$ from their corresponding multinomial full conditionals $$p(x_{j,t}=s\mid \mathit{\xi },x_{i,t}=s_i,i\ne j)\propto {\omega }_s^{\left(j\right)}\left(z_{j,t}\right){\mathit{\lambda }}_{{\varphi }_{s_1,\cdots ,s,\cdots ,s_q}}\left(y_t\right)$$ 7:    For $j=1,\cdots ,q$, collect samples $k_j$ from their respective multinomial full conditionals $$p(k_j=k\mid \mathit{\xi })\propto \exp (-{\mu }_{j}k)\prod _{c=1}^{C_j}{n}_{j,c}^{-k{\beta }_j},k_j=\underset{t}{max}\{x_{j,t}\},\cdots ,C_j$$    where $n_{j,c}={\sum }_{t=1}^T\mathbf{1}\{z_{j,t}=c\}$. 8: end for

To successfully run Algorithm 1, certain hyperparameters need to be chosen. The aforementioned determination of ${\mu }_j$ and L have been carefully discussed along with their implications, so we focus on the other hyperparameters. Among those hyperparameters, a and b will determine the clustering ability of the Pitman-Yor process (which are set to be 1 and 0 in this case), rendering it a Dirichlet process; this is sufficient for applications discussed in this paper. It should be noted that $\alpha$ and ${\beta }_j$ are Dirichlet Distribution’s hyperparamters and serve the role of pseudo-counts. The determination of these reflects the users’ prior belief. They are often manually chosen to be some small values without additional information which can justify larger values. In the following sections, they are chosen to be: $\alpha =1$ and ${\beta }_j=1/C_j$ across different applications.

3.2. Bayesian Factor and Hypothesis Testing

This subsection discusses hypothesis testing techniques on the significance of all the predictors to the regression Equation (1). It can be used to make causal inference in order to provide a better understanding of the model and to better allocate computational resources for the sequential classification task by including only the important predictors (and discard the unimportant ones). As previously noted, a particular predictor $z_j$ is considered important if and only if the number of clusters ${\tilde{k}}_j$ formed by their corresponding latent class allocation variables $x_j$ is greater than 1.

Let $\Lambda \subset \{1,\cdots ,q\}$ be the set of predictors under consideration. To perform the Bayesian hypothesis testing, we only need to compute the Bayes factor [23] in favor of $H_1:{\tilde{k}}_j>1$ for some $j\in \Lambda$ against $H_0:{\tilde{k}}_j=1$ for any $j\in \Lambda$, given by

$$B{F}_{10}=\frac{p\left(H_1\right|\mathit{y},\mathit{z})/p\left(H_1\right)}{p\left(H_0\right|\mathit{y},\mathit{z})/p\left(H_0\right)}$$

(15)

where $\mathit{y}\equiv {\{y_t\}}_{t=1}^T$, $\mathit{z}\equiv {\{{\mathit{z}}_t\}}_{t=1}^T$; and $p\left(H_0\right|\mathit{y},\mathit{z})$, $p\left(H_1\right|\mathit{y},\mathit{z})$ are numerically computed as the fraction of samples in which the ${\tilde{k}}_j$’s conform to $H_0$ and $H_1$, respectively; the prior probabilities $p\left(H_0\right)$ and $p\left(H_1\right)$ can be obtained by the following probability equation:

$$\begin{array}{cc}\hfill p({\tilde{k}}_j=1)& =\sum _{k=1}^{C_j}p(k_j=k)\sum _{l=1}^{k}p(x_{j,t}=l\forall t|k_j=k)\hfill \\ & =\left(\prod _{r=1}^{C_j}{\gamma }_j^{\left(n_{j,c}\right)}\right)\left(\sum _{k=1}^{C_j}\frac{p(k_j=k)k}{{\prod }_{k=1}^{C_j}{\left(k{\gamma }_j\right)}^{\left(n_{j,c}\right)}}\right)\hfill \end{array}$$

Specifically, to test whether $\theta$ Granger-causes y, it is only necessary to choose

$$\Lambda =\{D_1+1,\cdots ,q\}.$$

4. Sequential Classification

In Section 3, a Gibbs sampling algorithm is developed to infer the posterior distribution of model parameters given the observed data. In this section, a classification algorithm for dynamical systems based on the posterior predictive distribution, which is derived by marginalizing the likelihood of unobserved data over the posterior distribution of model parameters, is proposed. This algorithm consists of two phases: (1) off-line training phase and (2) online testing phase. Suppose there are M different classes of dynamical systems that are of interest, ${\mathcal{C}}_i,i=1,2,\cdots ,M$, for each of them we collect a training set ${}^{\left(i\right)}\mathcal{D}{}_{T_i}={\{{}^{\left(i\right)}y{}_t,{}^{\left(i\right)}\mathit{z}{}_t\}}_{t=1}^{T_i}$. The requirement for this dataset is that the data are categorical (e.g., quantized categories from continuous data), and for each class they have an identical number of categories of predictors and variables.

During the training phase, training set ${}^{\left(i\right)}\mathcal{D}{}_{T_i}$ is used to compute the posterior of samples

$${\{{}_{\left(n\right)}^{\left(i\right)}\mathit{\varphi },{}_{\left(n\right)}^{\left(i\right)}\mathit{\lambda },{}_{\left(n\right)}^{\left(i\right)}\mathit{\omega }\}}_{n=1}^M$$

for each one of the class ${\mathcal{C}}_i$, as previously described in Algorithm 1. Then, during the test phase, the test set ${\mathcal{D}}_T$ will be classified. Among these M classes, one will be identified as the class to which ${\mathcal{D}}_T$ most likely belongs. In order to do so, the following conditional probability $p({\mathcal{D}}_T\mid {}^{\left(i\right)}\mathcal{D}{}_{T_i})$ will be computed:

$$\begin{array}{c}p({\mathcal{D}}_T\mid {}^{\left(i\right)}\mathcal{D}{}_{T_i})=\prod _{t=1}^{T}p(y_t\mid {\mathit{z}}_t;{}^{\left(i\right)}\mathcal{D}{}_{T_i})\end{array}$$

(16)

$$\begin{array}{c}\begin{array}{cc}\hfill & p(y_t\mid {\mathit{z}}_t;{}^{\left(i\right)}\mathcal{D}{}_{T_i})\approx \frac{1}{N}\sum _{n=1}^N\left(\sum _{s_1=1}^{k_1}\cdots \sum _{s_q=1}^{k_q}{}_{\left(n\right)}^{\left(i\right)}\lambda {}_{{}_{\left(n\right)}^{\left(i\right)}\varphi s_1,\cdots ,s_q}\left(y_t\right)\prod _{j=1}^q{}_{\left(n\right)}^{\left(i\right)}\omega {}_{s_j}^{\left(j\right)}\left(z_{j,t}\right)\right)\hfill \end{array}\end{array}$$

(17)

Following the above calculation of conditional probabilities $p({\mathcal{D}}_T\mid {}^{\left(i\right)}\mathcal{D}{}_{T_i})$, the posterior probability of the test data ${\mathcal{D}}_T$ belonging to class ${\mathcal{C}}_i$ (denoted as $p({\mathcal{C}}_i\mid {\mathcal{D}}_T)$) can be then calculated as:

$$p({\mathcal{C}}_i\mid {\mathcal{D}}_T)=\frac{p({\mathcal{D}}_T\mid {}^{\left(i\right)}\mathcal{D}{}_{T_i})p\left({\mathcal{C}}_i\right)}{{\displaystyle \stackrel{M}{\sum _{r=1}}}p({\mathcal{D}}_T\mid {}^{\left(r\right)}\mathcal{D}{}_{T_r})p\left({\mathcal{C}}_r\right)}$$

(18)

where $p\left({\mathcal{C}}_i\right)$ is the prior probability of the class ${\mathcal{C}}_i$. Next, the classification result is generated by:

$$D_{\mathrm{class}}=\arg \underset{i}{max}p({\mathcal{C}}_i\mid {\mathcal{D}}_T)$$

(19)

The prior probability $p\left({\mathcal{C}}_i\right)$ reflects user’s subjective beliefs and can also be designed to optimize some objective criterion. The reason that the detection algorithm is “sequential” is due to the fact the conditional probability $p({\mathcal{D}}_T\mid {}^{\left(i\right)}\mathcal{D}{}_{T_i})$ is evaluated one by one as shown in Equation (16). In real-world applications, values of $p(y_t\mid {\mathit{z}}_t;{}^{\left(i\right)}\mathcal{D}{}_{T_i})$ in Equation (17) are often precomputed and stored for various values of $(y_t,{\mathit{z}}_t)$, in order to achieve faster computations.

For the binary classification case, we can construct the likelihood ratio test [24] as:

$$\frac{p({\mathcal{D}}_T\mid {}^{\left(1\right)}\mathcal{D}{}_{T_1})}{p({\mathcal{D}}_T\mid {}^{\left(0\right)}\mathcal{D}{}_{T_0})}\stackrel{1}{\underset{0}{\gtrless }}\Theta$$

(20)

where in this equation $\Theta$ is a certain threshold. To choose the threshold $\Theta$, one could rely on the receiver operating characteristic (ROC). ROC curves are often obtained by changing $\Theta$ in order to make a trade-off between the probability of (successful) detection $p_D=Prob(\mathrm{decide}1\mid 1\mathrm{is}\mathrm{true})$ and the false alarm probability $p_F=Prob(\mathrm{decide}1\mid 0\mathrm{is}\mathrm{true})$. Using those ROC curves, an optimal combination of $p_D$ and test set data length for a given $p_F$ can be selected, which would then determine the threshold $\Theta$.

5. Numerical Example

This section presents a numerical example which utilizes the proposed method to infer causal relationships between two categorical time series. In this example, the data generation model is known and thus can be compared with the results from the proposed algorithm for evaluation of performance. The data generation details are given below.

In this particular numerical example, there are two binary sequences of symbols $y_t$ and ${\theta }_t$. Symbol sequences $y_t$ are generated using a known Markov model $p(y_t\mid y_{t-1},y_{t-3},y_{t-4})$, where only the time-lags $y_{t-1},y_{t-2},y_{t-5}$ are important predictors. Symbol sequences ${\theta }_t$ are generated from another Markov model $p({\theta }_t\mid {\theta }_{t-1},{\theta }_{t-2},y_{t-1},y_{t-3})$, where ${\theta }_{t-1},{\theta }_{t-2}$ and $y_{t-1},y_{t-3}$ are the key predictors. In other words, the variable y Granger-causes the variable $\theta$ but not the other way around because y only depends on its own past.

Table 1. Transition Probabilities for $y_t$ in the Numerical Example.

${\mathit{y}}_{\mathit{t}-1}$	${\mathit{y}}_{\mathit{t}-3}$	${\mathit{y}}_{\mathit{t}-4}$	$\mathit{p}({\mathit{y}}_{\mathit{t}}=1)$	$\mathit{p}({\mathit{y}}_{\mathit{t}}=0)$
0	0	0	0.20	0.80
1	0	0	0.75	0.25
0	1	0	0.70	0.30
1	1	0	0.35	0.65
0	0	1	0.40	0.60
1	0	1	0.38	0.62
0	1	1	0.33	0.67
1	1	1	0.71	0.29

lists the transition probabilities for $y_t$, where it is seen that the predictors are $y_{t-1},y_{t-3},y_{t-4}$ only.

Table 2. Transition Probabilities for ${\theta }_t$ in the Numerical Example.

${\mathit{y}}_{\mathit{t}-1}$	${\mathit{y}}_{\mathit{t}-3}$	${\mathit{\theta }}_{\mathit{t}-1}$	${\mathit{\theta }}_{\mathit{t}-2}$	$\mathit{p}({\mathit{\theta }}_{\mathit{t}}=1)$	$\mathit{p}({\mathit{\theta }}_{\mathit{t}}=0)$
0	0	0	0	0.40	0.60
1	0	0	0	0.65	0.35
0	1	0	0	0.70	0.30
1	1	0	0	0.40	0.60
0	0	1	0	0.50	0.50
1	0	1	0	0.47	0.53
0	1	1	0	0.33	0.67
1	1	1	0	0.69	0.31
0	0	0	1	0.45	0.55
1	0	0	1	0.75	0.25
0	1	0	1	0.30	0.70
1	1	0	1	0.50	0.50
0	0	1	1	0.75	0.25
1	0	1	1	0.66	0.34
0	1	1	1	0.65	0.35
1	1	1	1	0.20	0.80

lists the transition probabilities for ${\theta }_t$, where the predictors are $y_{t-1},y_{t-3},{\theta }_{t-1}$, and ${\theta }_{t-2}$ only.

To estimate the regression model in Equation (1) with the parameter $T=1005$, samples of ${\{y_t\}}_{t=1}^{1005}$ and ${\{{\theta }_t\}}_{t=1}^{1005}$ are being collected simultaneously. Based on the prior belief that $y_{t-D}$ and ${\theta }_{t-D}$ are no longer important for making predictions about $y_t$ and ${\theta }_t$ when D is greater than 5, predictors for both $y_t$ and ${\theta }_t$ are set as follows:

$${\mathit{z}}_t\equiv (y_{t-1},y_{t-2},y_{t-3},y_{t-4},y_{t-5},{\theta }_{t-1},{\theta }_{t-2},{\theta }_{t-3},{\theta }_{t-4},{\theta }_{t-5})$$

(21)

From these data sets, 1000 training samples are chosen for testing the proposed algorithm.

To calculate posteriors using Algorithm 1 for $p\left(y_t\right|{\mathit{z}}_t)$, since there is no other prior knowledge, ${\mu }_j$ is set to be 1 across $j=1,\cdots ,10$. Initially, 200,000 samples are used in a burn-in period: they are fed into the algorithm and then discarded. The next 50,000 samples (after burn-in) are downsampled further by taking every 5th sample to reduce their autocorrelation. Figure 2 summarizes the results, in which Figure 2a displays the log-likelihood for 10,000 iterations of this model and Figure 2b illustrates the ability to correctly identify all the important predictors for the proposed method. For this example, the key predictors should be 1, 3 and 4, and the results from the prediction ($y_{t-1},y_{t-3}$ and $y_{t-4}$) are the same as the ground truth. Figure 2c shows the relative frequency of number of predictors that are important. Furthermore, the proposed method also creates parsimonious representations of the model as seen in Figure 2d,e. As previously discussed in Section 2.2, the tensor ${\mathit{\lambda }}_{s_1\cdots s_q}\left(y_t\right)$ has more components than needed but it can be clustered in a nonparametric way to reduce the number of combinations. Referring to [13], Figure 2f shows the Bayes factors calculation as mentioned in Section 3.2 for all of the predictors. Bayes factor $B{F}_{10}$ in Equation (15) can be regarded as the evidence against $H_0$. After setting a commonly-used threshold of $t=20$, it can be concluded that those predictors with higher $B{F}_{10}$ have implications of their evidences being strong. Furthermore, having $B{F}_{10}>150$ indicates even stronger evidence against the hypothesis $H_0$ [13]. It should be noted here that when the inclusion proportions of different lags in Figure 2b are equal to 1, then their corresponding Bayes factors in Figure 2f should tend to infinity (as for predictors $1,3$ and 4 in this example).

Similarly, Figure 3 shows the results using the same set of data as in Figure 2 but instead of estimating $p\left(y_t\right|{\mathit{z}}_t)$, we are estimating $p\left({\theta }_t\right|{\mathit{z}}_t)$ here. Figure 3a–f have the same implications as those previously stated for Figure 2a–f. It can be seen that in this case for $p\left({\theta }_t\right)$, the key predictors should be $1,3,6$ and 7, and the results confirm this in Figure 3.

Besides the ability to correctly identify the structure of the model, the proposed method can also perform transition probability estimation. Figure 4 illustrates two arbitrarily selected cases from Table 1 and Table 2. Setting $y_{t-1}=0$, $y_{t-3}=1$, and $y_{t-4}=0$, from Table 1 we can get the transition probability of the model of $(y_t=1)$ is $0.70$. Similarly, setting $y_{t-1}=1$, $y_{t-3}=0$, ${\theta }_{t-1}=1$, and ${\theta }_{t-2}=0$, from Table 2 we can get the transition probability of $(y_t=1)$ is $0.50$. In Figure 4, the estimated transition probability using the proposed method is displayed along with their running mean as well as their $5\%$ and $95\%$ percentiles. From both subplots of Figure 4, it is observed that the running mean of the transition probability is actually close to the true transition probability as given in the data generation tables. Even with a limited amount of data, the proposed method can not only estimate the transition probabilities, but also give an uncertainty bound in terms of their respective quantiles.

The causal relationship between $\mathit{y}$ and $\theta$ is identified by Bayes factor analysis (see Section 3.2. The results are summarized in

Table 3. Hypothesis Testing of Granger causality in the Numerical Example.

Null Hypothesis	Bayes Factor ${\mathbf{BF}}_{10}$
$\theta$ does not Granger-cause y	0.43
y does not Granger-cause $\theta$	Infinity

, which show that y Granger-causes $\theta$ but not the other way, which is in line with the ground truth.

6. Validation with Experimental Data: The Combustor Apparatus

This section validates the nonparametric regression model with experimental data generated from a swirl-stabilized lean-premixed laboratory-scale combustor apparatus [14].

6.1. Background and Description of the Experimental Procedure

This subsection presents a brief background of thermoacoustic instabilities in the combustor apparatus along with the experimental details for data collection. Thermoacoustic instabilities occur from highly nonlinear coupled phenomena that evolve from mutual interactions among thermofluid dynamics, unsteady heat release, and acoustics of the combustor chamber. The resulting self-sustained high-amplitude pressure oscillations often impose severe negative impacts on the performance and operational life of gas turbine engines [25,26,27].

Technical literature abounds with studies on combustion instabilities and their early detection by time series analysis, especially by using Markov chains [28,29]. However, current methods are largely limited to individual investigations of pressure or chemiluminescence measurements, and have apparently not taken the machine-learning-theoretic approach to information fusion into consideration; consequently, fast detection of thermoacoustic instabilities may not be achieved to the full extent based on the individual information of different sources only. Moreover, parameter estimation is difficult in current methods, even for moderately high-order Markov chains, due to the paucity of data, let alone a more sophisticated information fusion model. As for the detection procedure, empirical thresholds are often used in existing literature, without taking advantage of methods in statistical detection theory (such as sequential testing techniques); therefore, those applications are very limited in real-time detection cases.

Figure 5 presents a schematic diagram of the combustor apparatus [14] that 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.

Experiments have been conducted at 62 different operating conditions by varying the equivalence ratio and percentage of pilot fuel, as listed in

Table 4. Operating conditions.

	Parameters	Values
Variables	Equivalence Ratio	0.525, 0.538, 0.575, 0.625
	Pilot Fuel (percent)	0–9% (0.5% increment)
Fixed Conditions	Inlet Temperature	$250{}^{\circ }\mathrm{C}$
	Inlet Velocity	40 m/s
	Combustor Length	0.625 m

. Under each operating condition, 8 s of pressure and chemiluminescence measurements have been collected at the sampling rate of 8192 Hz, where stable and/or unstable modes are recorded along with each time series data. To alleviate the problem of (possible) oversampling, the pressure and chemiluminescence measurements from combustors are first downsampled, which is obtained from first minimum of the average mutual information [30]. Then, the continuously varying time series data for both stable and unstable modes are quantized using maximum entropy partitioning [31,32] with a ternary alphabet $\Sigma =\{1,2,3\}$. The quantized pressure measurements are denoted as $y_t$ and the chemiluminescence measurements are denoted as ${\theta }_t$ at time instant t.

6.2. Training Phase

This subsection describes details in the nonparametric regression model training, wherein 500 samples have been used after downsampling the quantized pressure time series data under stable and unstable conditions. The maximum memory D of each of $y_t$ and ${\theta }_t$ in this dataset is observed to be generally limited to 5 for both stable and unstable cases. Hence, predictors of $y_t$ or ${\theta }_t$ are set to be ${\mathit{z}}_t\equiv (y_{t-1},y_{t-2},\cdots ,y_{t-5},{\theta }_{t-1},{\theta }_{t-2},\cdots ,{\theta }_{t-5})$ and the corresponding regression model is hereafter referred to as “full order model”. Since $y_t$ and ${\theta }_t$ has three categories, it follows that $C_y=C_{\theta }=3$.

To compute posteriors, as in Algorithm 1, the values

$$[1,1.5,2.0,2.5,3.0,1.0,1.5,2.0,2.5,3.0]$$

are assigned to ${\mu }_j$ for $j=1,\cdots ,10$. After discarding 200,000 data points during the burn-in period, remaining 50,000 samples are then downsampled by taking every 5th data point to reduce their autocorrelation. Gibbs sampling results of pressure data are represented as $p(y_t\mid y_{t-1},\cdots ,y_{t-5},{\theta }_{t-1},\cdots ,{\theta }_{t-5})$ in Figure 6a,b for a stable mode and in Figure 6c,d for an unstable mode. Similarly, Gibbs sampling results of chemiluminescence data are represented as $p({\theta }_t\mid y_{t-1},\cdots ,y_{t-5},{\theta }_{t-1},\cdots ,{\theta }_{t-5})$ in Figure 7a,b for a stable mode and in Figure 7c,d for an unstable mode.

Figure 6a,c and Figure 7a,c show the log likelihood with different iterations for pressure and chemiluminescence data under stable and unstable conditions, respectively. Similarly, Figure 6b,d and Figure 7b,d illustrate the Bayes factors of predictors for pressure and chemiluminescence data under stable and unstable conditions, respectively. Based on the Bayes factor analysis, the important predictors for stable pressure data are identified as:

$$y_{t-1},y_{t-2},y_{t-3},y_{t-4},{\theta }_{t-1},{\theta }_{t-3}\mathrm{and}{\theta }_{t-4}$$

while those for unstable pressure data are identified as

$$y_{t-1},y_{t-3},y_{t-4},y_{t-5},{\theta }_{t-1},{\theta }_{t-3},{\theta }_{t-4}\mathrm{and}{\theta }_{t-5}$$

Using the identical set of hyperparameters and number of iterations, Gibbs sampling has been performed on the same set with pressure data $y_t$ only; this is referred to as the “reduced order model” in the following text. In this case the predictors are set as ${\mathit{z}}_t\equiv (y_{t-1},y_{t-2},\cdots ,y_{t-5})$. The stable and unstable cases are shown in Figure 8a–d respectively. The important predictors for $y_t$ using this reduced order model are: $y_{t-2},y_{t-4}$, and $y_{t-5}$ for the stable mode, and $y_{t-1},y_{t-2},y_{t-4}$, and $y_{t-5}$ for the unstable mode.

Similarly, for chemiluminescence data, the important predictors are identified as:

$$y_{t-1},y_{t-3},y_{t-4},y_{t-5},{\theta }_{t-1},{\theta }_{t-2},{\theta }_{t-3},{\theta }_{t-4}\mathrm{and}{\theta }_{t-5}$$

while those for unstable chemiluminescence data are identified as:

$$y_{t-1},y_{t-2},{\theta }_{t-1},{\theta }_{t-2},{\theta }_{t-3},{\theta }_{t-4}\mathrm{and}{\theta }_{t-5}$$

6.3. Granger Causality

To identify the Granger causal relationship between pressure and chemiluminescence data, Bayes factor analysis has been performed for both stable and unstable cases as described in Section 3.2. The results are summarized in

Table 5. Hypothesis Testing of Granger Causality.

Null Hypothesis	Operating Condition	${\mathbf{BF}}_{10}$
$\theta$ does not Granger-cause y	Stable	Infinity
y does not Granger-cause $\theta$	Stable	Infinity
$\theta$ does not Granger-cause y	Unstable	Infinity
y does not Granger-cause $\theta$	Unstable	Infinity

, which show that pressure and chemiluminescence measurements Granger-cause each other under both stable and unstable conditions; this implies that fusion of these two measurements can enhance the accuracy of prediction. This kind of mutual interaction between pressure and chemiluminescence measurements could be caused by a third unknown physical quantity, the exploration of which is a topic of future research.

6.4. Sequential Classification

For evaluation of the performance of the sequential classification for thermoacoustic instability identification, 100 instances of 50-sample datasets, which are not included in the training set, have been selected (also from their downsampled quantized pressure measurements for both stable and unstable modes). Figure 9 exhibits the profiles of posterior probability of each class as a function of the length of the observed data, where the top plate (i.e., Figure 9a) uses the full order model, and the bottom plate (i.e., Figure 9b) uses the reduced-order model for the same test data sequence. While the test sequences are correctly classified by both models, the reduced-order model is slower than the full-order model that contains more information.

Figure 10 shows the receiver operating characteristic (ROC) curves for the proposed detection algorithm with different lengths of the test data. These ROC curves are plotted for both full-order and reduced-order models to show that, when testing with the same dataset, the full order model achieves better detection performance in terms of the area under the ROC. In other words, the full-order model may achieve the same performance as the reduced order model in a shorter time, which is desirable for active control of thermoacoustic instabilities in real time. It is also observed that the ROC curves tend to improve (i.e., move toward the top left corner) considerably as the length of test data is increased from 5 to 9. This is expected because the information contents monotonically increase with the length of test data and hence better results are obtained.

7. Validation with Economics Data

This section validates the nonparametric regression model with (publicly available) real-world economics data. Specifically, monthly data of the U.S. consumer price index (CPI) and the U.S. Dollar London Interbank Offered Rate (LIBOR) interest rate index with one-month maturity from January 1986 to December 2016 are used. It is noted that: (i) U.S. CPI is a measure of the average change over time in the prices paid by urban consumers for U.S. market of consumer goods and services, and (ii) U.S. Dollar LIBOR is a benchmark for short-term interest rates around the world, which is not a monetary measure associated with any country, and which does not reflect any institutional mandate in contrast to, e.g., when the Federal Reserve sets interest rates. Economics theory [33] indicates that low interest rates can cause high inflation, and empirical research [34] has been conducted to investigate the causal relationship between inflation and nominal or real interest rates for the same country or region.

To avoid spurious regression [35], the raw data of U.S. CPI and U.S. Dollar LIBOR are preprocessed to achieve stationarity. U.S. CPI raw data are used to calculate the monthly percentage increase, and then this percentage increase is converted into a categorical variable by discretizing to quintiles (e.g., 5-quantiles in this study) that are denoted as $y_t$; the rationale for discretization of (noise-contaminated) continuously varying data is to improve the signal-to-noise ratio [36]. Similarly, U.S. LIBOR raw data are used to calculate the monthly difference, and then this difference is convertedin to a categorical variable by discretizing to quintiles, denoted as ${\theta }_t$. The entire dataset is used for training the proposed algorithm.

To estimate the regression model in Equation (1), based on the assertion that $y_{t-D_y}$ and ${\theta }_{t-D_{\theta }}$ are not important for predicting $y_t$ and ${\theta }_t$ if both $D_y$ and $D_{\theta }$ are greater than 6 (i.e, six months for both CPI and LIBOR), the predictor for $y_t$ and ${\theta }_t$ is set as:

$${\mathit{z}}_t\equiv (y_{t-1},y_{t-2},\cdots ,y_{t-6},{\theta }_{t-1},{\theta }_{t-2},\cdots ,{\theta }_{t-6})$$

(22)

To compute the posterior probabilities using the proposed Algorithm 1, ${\mu }_j$ are assigned to be $j/2$ for $j=1,\cdots ,6$ and $(j-6)/6$ for $j=7,\cdots ,12$. After the initial 100,000 samples are discarded during the burn-in period, the remaining 50,000 samples are then downsampled by taking every 5th to reduce their autocorrelation. Figure 11 and Figure 12 respectively summarize the results for $y_t$ and ${\theta }_t$. These figures have similar characteristics to their counterparts in the numerical example in Section 5. The results show that, for $y_t$ or CPI, the important lags are $y_{t-1},y_{t-2},y_{t-3}$ and ${\theta }_{t-1}$. Similarly, for ${\theta }_t$ or LIBOR, the important lags are ${\theta }_{t-1},{\theta }_{t-1},{\theta }_{t-3}$. These results show that LIBOR Granger-cause CPI, but not vice versa. This conclusion is summarized by Bayes factor analysis in

Table 6. Hypothesis test of Granger Causality for economics data.

Null Hypothesis	Bayes Factor ${\mathbf{BF}}_{10}$
US CPI does not Granger-cause LIBOR	7.29
LIBOR does not Granger-cause US CPI	Infinity

.

8. Summary, Conclusions, and Future Work

The proposed Bayesian nonparametric method provides a flexible model for information fusion of heterogeneous, correlated time series data. The proposed method has been validated on a real-world application by using the experimental data collected from a laboratory-scale swirl-stabilized combustor apparatus, as well as on the publicly available economics data. It is demonstrated that the proposed method is capable of enhancing the accuracy for real-time detection of thermoacoustic instabilities and correctly identifying the Granger causal relationship between key economic variables.

There are many promising directions in which the proposed model can be further explored, such as:

• Variational inference algorithm development for the proposed model [37].
• Extension of the present analysis to hidden Markov models (HMM) [38] and information transfer [39].
• Exploration of an unknown physical quantity that may cause the appearance of mutual interactions between pressure and chemiluminescence measurements.
• Investigation of the empirical performance of the proposed approach utilizing extensive simulation studies.