Variational Beta Process Hidden Markov Models with Shared Hidden States for Trajectory Recognition

Abstract

Hidden Markov model (HMM) is a vital model for trajectory recognition. As the number of hidden states in HMM is important and hard to be determined, many nonparametric methods like hierarchical Dirichlet process HMMs and Beta process HMMs (BP-HMMs) have been proposed to determine it automatically. Among these methods, the sampled BP-HMM models the shared information among different classes, which has been proved to be effective in several trajectory recognition scenes. However, the existing BP-HMM maintains a state transition probability matrix for each trajectory, which is inconvenient for classification. Furthermore, the approximate inference of the BP-HMM is based on sampling methods, which usually takes a long time to converge. To develop an efficient nonparametric sequential model that can capture cross-class shared information for trajectory recognition, we propose a novel variational BP-HMM model, in which the hidden states can be shared among different classes and each class chooses its own hidden states and maintains a unified transition probability matrix. In addition, we derive a variational inference method for the proposed model, which is more efficient than sampling-based methods. Experimental results on a synthetic dataset and two real-world datasets show that compared with the sampled BP-HMM and other related models, the variational BP-HMM has better performance in trajectory recognition.

1. Introduction

Trajectory recognition is important and meaningful in many practical applications, such as human activities recognition [1], speech recognition [2], handwritten character recognition [3] and navigation task with mobile robot [4]. In most practical applications, the trajectory is affected by the hidden features corresponding to each point. The hidden Markov model (HMM) [2], hierarchical conditional random field (HCRF) [5,6] and the HMM-based models, such as the hierarchical Dirichlet process hidden Markov model (HDP-HMM) [7], the Beta process hidden Markov model (BP-HMM) [8,9,10] and the Gaussian mixture model hidden Markov model (GMM-HMM) [2] are used to model sequential data and identify their classes [11,12,13,14].

The HMM is a popular model which has been applied widely in human activity recognition [1,15], speech recognition [2,16] and remote target tracking [2,17]. Besides, the HMM is becoming a more significant part as a building block of smart cities and Industry 4.0 [18,19] and implemented in extensive applications such as driving behaviors prediction [20] and the inernet of thing (IoT) signature anomalies [21]. One drawback of the HMM is having to ensure in advance the number of hidden states that need to be selected or cross-validated. To address this problem, several methods based on model selection are employed, such as BIC [22] or some Bayesian non-parameter prior like the BP [23] and the HDP [24]. Besides, directly using the original HMM for classification has another disadvantage, in which each HMM is trained for one class separately and thus information from different classes cannot be shared. It is worth mentioning that the sampled BP-HMM proposed by Fox et al. [9] can not only learn the number of hidden features automatically but also obtain the sharing features between different classes, which has been proved to be meaningful for human activity trajectory recognition. The sampled BP-HMM learns the shared states among different classes by jointly modeling all trajectories together, in which a hidden state indicator for one trajectory with a BP prior is introduced and thus a state transition matrix for each trajectory is maintained. When used for classification, the sampled BP-HMM calculates the class-specific transition matrix by averaging the transition matrices of the trajectories from the corresponding class. However, from the perspective of performance or efficiency, if the sampled BP-HMM [1,7] is used for classification, there is still a lot of room for improvement.

From the perspective of performance, the classification procedure in the sampled BP-HMM [1] is too rough to make full use of the trained model, in which the state transition matrix for each class is calculated by averaging the transition matrixes of all the trajectories. Obviously, this will lead to the loss of information, especially when the training set has some ambiguous trajectories. For instance, a “running” class has some “jogging” trajectories. One naive method to solve it is to select the K best HMMs for each class. However, it will cost plenty of time to select representatives for each class. In order to take account of both performance and efficiency, we change the way of modeling data in BP-HMMs. Differently from those versions of BP-HMMs [1,8,9,10,25], in variational BP-HMMs, an HMM is created for each class instead of for each trajectory.

From the perspective of efficiency, the existing approximate inference for the BP-HMM is based on sampling methods [1,9] which often converge slowly. This drawback of the sampled BP-HMM [1] is inconvenient to practical applications. To provide a faster convergence rate than sampling methods, we develop variational inference for the BP-HMM. If the variational lower bound is unchanged or almost unchanged, the iteration will stop. To be amenable to the variational method, we use the stick-breaking construction of the BP [26] instead of the Indian buffet process (IBP) construction [27] in the sampled BP-HMM.

In this paper, we propose a variational BP-HMM for trajectory recognition, in which the way of the data modeling and the inference method are novel compared with the previous sampled BP-HMM. On the one hand, the new method of modeling trajectories enables the model to obtain better classification performance. Specifically, the hidden state can be optionally shared, and the class-specific state indicator is more suitable for classification than the trajectory-specific state indicator in the sampled BP-HMM. The transition matrix is actually learned from the data instead of averaging all the trajectory-specific transitions. On the other hand, the derived variational inference of the BP-HMM makes the model more efficient. In particular, we use the two-parameter BP as the prior of the class-specific state indicator, which is more flexible than the one-parameter Indian buffet process in the sampled BP-HMM. We apply our model to the navigation task of mobile robots and human activity trajectory recognition. Experimental results on the synthetic and real-world data show that the proposed variational BP-HMM with sharing hidden states has advantages to trajectory recognition.

The remainder of this paper is organized as follows. Section 2 gives an overview of the BP and HMM. In Section 3, we review the model assumption of the sampled BP-HMM. In Section 4, we present the proposed variational BP-HMM including the model setting and its variational inference procedure. Experimental results on both synthetic and real-world datasets are reported in Section 5. Finally, Section 6 gives the conclusion and future research directions.

2. Preliminary Knowledge

In order to explain the variational BP-HMM more clearly, the key related backgrounds including BP and HMM will be introduced in the following sub-sections.

2.1. Beta Process

The BP is defined by Hjort [28] for applications in survival analysis. It is a significant application as a non-parametric prior for latent factor models [23,26], and used as a non-parameter prior for selecting the hidden state set of the HMM [8,9,25]. At the beginning, the BP is defined on the positive real line $\left({\mathbb{R}}^{+}\right)$ then extended to more general spaces $\mathsf{\Omega }$ (e.g., $\mathbb{R}$).

A BP, $B\sim \mathrm{BP}\left(\alpha ,B_0\right)$, is a positive Lévy process. Here, $\alpha$ is the concentration parameter and $B_0$ is a fixed measure on $\mathsf{\Omega }$. Let $\gamma =B_0\left(\mathsf{\Omega }\right)$. The $\mathrm{BP}\left(\alpha ,B_0\right)$ is formulated as

$$\begin{array}{c}\hfill \begin{array}{cc}& B_K=\sum _{k=1}^{\infty }{\pi }_k{\delta }_{{\omega }_k},\hfill \\ & {\omega }_{i_j}\stackrel{i.i.d.}{\sim }\frac{1}{\gamma }{B}_0,\hfill \end{array}\end{array}$$

(1)

where $\left\{\omega \right\}$ are atoms in B. If $B_0$ is continuous, the Lévy measure of the BP is expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}& \nu \left(d\omega ,d\pi \right)=\alpha \left(\omega \right){\pi }^{-1}{\left(1-\pi \right)}^{c\left(\omega \right)-1}d\pi B_0\left(d\omega \right).\hfill \end{array}\end{array}$$

(2)

If $B_0$ is discrete, in the form of $B_0={\sum }_k{q}_k{\omega }_k$, the atoms in B and $B_0$ have the same location. It can be represented as follows

$$\begin{array}{c}\hfill \begin{array}{cc}& B_K=\sum _{k=1}^K{\pi }_k{\delta }_{{\omega }_k},\hfill \\ & {\pi }_k\stackrel{i.i.d.}{\sim }\mathrm{Beta}\left(\frac{\alpha \gamma }{K},\alpha \left(1-\frac{\gamma }{K}\right)\right),\hfill \\ & {\omega }_k\stackrel{i.i.d.}{\sim }\frac{1}{\gamma }{B}_0.\hfill \end{array}\end{array}$$

(3)

As $K\to \infty$ and $H_K\to \infty$, B represents a BP [29].

The BP is conjugate to a class of Bernoulli process, denoted by $\mathrm{BeP}\left(B\right)$. For example, we define a Bernoulli process $F\sim \mathrm{BeP}\left(B\right)$. In this article, we focus on the discrete Bernoulli process in the form of $B={\sum }_k{\pi }_k{\delta }_{{\omega }_k}$, and then the Bernoulli process can be expressed as $F={\sum }_k{b}_k{\delta }_{{\omega }_k}$, where $B\in [0,1]$, $b_k$ is the independent Bernoulli variable with the probability ${\pi }_k$. If B is a BP, then

$$\begin{array}{c}\hfill \begin{array}{cc}& B\sim \mathrm{BP}\left(\alpha ,B_0\right),\hfill \\ & F\sim \mathrm{BeP}\left(B\right),\hfill \end{array}\end{array}$$

(4)

is called the Beta-Bernoulli process.

Similarly to Dirichlet process which has two principle methods for drawing samples, (1) the Chinese restaurant process [30], (2) the stick-breaking process [31], the BP generates samples using the Indian buffet process (IBP) [23] and the stick-breaking process [29].

The original IBP can be seen as a special case of the general BP, i.e., an IBP is a one-parameter BP. Similarly to the Chinese restaurant process, the IBP is described in the view of customers choosing dishes. It is also employed to construct two-parameter BPs but with some details changed. Specifically, the procedure for constructing BP$(\alpha ,B_0),\gamma =B_0\left(\mathsf{\Omega }\right)$ is as follows:

• The first customer takes the first $\mathrm{Poisson}\left(\gamma \right)$ dishes.
• The nth customer then takes dishes that have been previously sampled with probability $\frac{m_k}{\alpha +n-1}$, where $m_k$ is the number of people who have already sampled the dish k. He also takes Poisson $\left(\frac{\alpha \gamma }{\alpha +n-1}\right)$ new dishes.

The BP has been shown as a de Finetti mixing distribution underlying the Indian buffet process, and an algorithm has been presented to generate the BP [23].

The stick-breaking process of the BP, $B\sim \mathrm{BP}\left(\alpha ,B_0\right)$, is provided by Paisley et al. [29]. It is formulated as follows.

$$\begin{array}{c}\hfill \begin{array}{cc}& B=\sum _{i=1}^{\infty }\sum _{j=1}^{C_i}{V}_{i_j}^{\left(i\right)}\prod _l^{i-1}\left(1-V_{i_j}^{\left(l\right)}\right){\delta }_{{\omega }_{i_j}},\hfill \\ & C_i\stackrel{i.i.d.}{\sim }\mathrm{Poisson}\left(\gamma \right),\hfill \\ & V_{i_j}^{\left(l\right)}\stackrel{i.i.d.}{\sim }\mathrm{Beta}\left(1,\alpha \right),\hfill \\ & {\omega }_{i_j}\stackrel{i.i.d.}{\sim }\frac{1}{\gamma }{B}_0.\hfill \end{array}\end{array}$$

(5)

It is clearly shown from the above equations that in every round (indexed by i), $C_i$ atoms have been selected, the weights of them follow an i-times stick-breaking process in which each breaking has the $\mathrm{Beta}\left(1,a\right)$ probability and $C_i$ is drawn from $\mathrm{Poisson}\left(\gamma \right)$.

2.2. Hidden Markov Models

The HMM [2] is a state space model where each sequence uses a Markov chain of discrete latent variables, with each observation conditioned on the state of the corresponding latent variables. Obviously, they are appropriate to model the data varying over time, and the data can be considered to be generated by the process that switches between different phases or states at different time-points. The HMM has been proved as a valuable tool in human activity recognition, speech recognition and many other popular areas [32].

Suppose that the trajectory observation $X=\{x_1,\dots ,x_N\}$ is an $N\times d$ matrix and $Z=\{z_1,\dots ,z_N\}$ is a N dimensional latent variable vector which has a value set ${\mathsf{\Omega }}_1$ with size K. The joint distribution of X and Z is expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p\left(X,Z\mid \theta \right)=& p\left(z_1|\mathbf{\pi }\right)\left\{\prod _{t=2}^{T}p\left(z_t|z_{t-1},\mathsf{\Pi }\right)\right\}\hfill \\ & \prod _{t=1}^{T}p\left(x_t|z_t,\varphi \right),\hfill \end{array}\end{array}$$

(6)

where $\theta =\{{\mathbf{\pi }}_0,{\mathbf{\pi }}_k,\varphi \}$, and A is a $K\times K$ matrix with ${\pi }_{jk}=p\left(z_{t+1}=k|z_t=j\right)$, $t=\{1,\dots ,T-1\},i,j\in {\mathsf{\Omega }}_1$ with ${\sum }_k{\pi }_{jk}=1$, and ${\pi }_{0k}$ is a K dimensional vector with ${\pi }_{0k}=p\left(z_1=k\right),k\in {\mathsf{\Omega }}_1$ with ${\sum }_k{\pi }_{0k}=1$. Furthermore, in the nonparametric version of HMM, the matrix $\mathsf{\Pi }$ can be assumed to obey a Dirichlet distribution, i.e., ${\mathbf{\pi }}_j\sim \mathrm{Dir}\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K\right)$ where $\sum {\alpha }_k=1.$ The probabilistic graphical model is represented in Figure 1.

If ${\mathbf{x}}_t$ is discrete with value set ${\mathsf{\Omega }}_2$ in the size of D, $\varphi$ is a $K\times D$ matrix with element ${\varphi }_{i_j}=p\left({\mathbf{x}}_t=j|z_t=i\right),i\in {\mathsf{\Omega }}_1,j\in {\mathsf{\Omega }}_2$. $\mathsf{\Pi }$ and $\varphi$ are named respectively as the transition matrix and emission matrix. If ${\mathbf{x}}_t$ is continuous, the emission matrix will be replaced by the emission distribution, where $\varphi$ is often defined as a distribution like Gaussian distribution $p\left({\mathbf{x}}_t|z_t=k\right)=\mathcal{N}\left({\mathbf{x}}_t|{\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k\right)$, $k\in {\mathsf{\Omega }}_1$. In the fully Bayesian framework, ${\mu }_k,{\mathsf{\Sigma }}_k$ can be regarded as random variables with distribution like normal inverse Wishart or Gaussian with Gamma distribution.

Marginal likelihood is often used to evaluate how an HMM is fit for the trajectories. Therefore, the HMM is usually trained by maximizing the marginal likelihood over the training trajectories. Baum–Welch (BW) algorithm, as an EM method is a famous algorithm for learning parameters of HMMs. The parameters include the transition matrix, the initial state distribution and the emission matrix (distribution’s parameters). In the BW algorithm, the forward-backward algorithm is employed to calculate the marginal probability. It should be noted that since the BW algorithm can only find the local optimum, multiple initializations are usually used to obtain better solutions. Given the learned parameters, the most likely state sequence corresponding to a trajectory is required in many practical applications. Viterbi algorithm is an effective method to obtain the most probable state sequence.

HMMs are a kind of generative model; they model the distribution of all the observed data. In trajectory classification tasks, such as activity trajectory recognition, different HMMs are used to model different classes of trajectories separately. After training these HMMs, the parameters in different HMMs are used to evaluate the newly come trajectory to find the most probable class. Specifically, to model large multiple trajectories from different classes, a separate HMM is defined for each class of trajectories, where ${\theta }^c$ represents its parameters. Given the trained HMMs, the class label $y^{*}$ of a new test trajectory ${\mathbf{x}}^{*}$ is determined according to

$$y^{*}=\underset{c}{argmax}lnp\left({\mathbf{x}}^{*}\right|{\theta }^c),$$

(7)

where $p\left({\mathbf{x}}^{*}|{\theta }^c\right)$ can be calculated using the forward-backward algorithm.

3. The Sampled BP-HMM

The sampled BP-HMM [9] is proposed to discover the available hidden states and the sharing patterns among different classes. It jointly models multiple trajectories and learns a state transition matrix for each trajectory. The sampled BP-HMM is successfully applied to trajectory recognition tasks, such as human activity trajectory recognition [1,10].

The sampled BP-HMM uses HMMs to model all the trajectories from all the classes and uses the BP as the prior of the indicator variables with each one corresponding to one trajectory. Suppose $X=\{X^{\left(1\right)},X^{\left(2\right)},\dots ,X^{\left(n\right)},\dots ,X^{\left(N\right)}\},N\in {\mathbb{N}}^{+}$ where $X^{\left(n\right)}$ is the $n\mathrm{th}$ trajectory. Each trajectory is modeled by an HMM. These HMMs share a global hidden feature set $\mathsf{\Omega }$ with the size of ∞. The sampled BP-HMM uses a hidden state selection matrix F with the size of $N\times \infty$ to indicate the available states for each trajectory, i.e., $f_{nk}=\{0,1\}$ indicators whether the $n\mathrm{th}$ HMM owns the $k\mathrm{th}$ state. The prior of the transition matrix ${\mathsf{\Pi }}^{\left(n\right)}$ for each trajectory is related to F, The transition matrix of the $n\mathrm{th}$ HMM is

$${\mathbf{\pi }}_j^{\left(n\right)}\sim Dir\left([r,r,\dots ,r+\kappa ,r,\cdots ]\odot {\mathbf{f}}_n\right),j>1,$$

(8)

and the initial state probability vector ${\mathbf{\pi }}_0^{\left(n\right)}$ is also related to F,

$${\mathbf{\pi }}_0^{\left(n\right)}\sim Dir\left([r,r,\dots ,r,r,\cdots ]\odot {\mathbf{f}}_n\right).$$

(9)

Similarly to the standard HMM, the latent variable ${\mathbf{z}}^{\left(n\right)}$ is a discrete sequence with

$$z_1^{\left(n\right)}\sim {\mathbf{\pi }}_0^{\left(n\right)},z_{t+1}^{\left(n\right)}|z_t^{\left(n\right)}\sim {\mathbf{\pi }}_{z_t^{\left(n\right)}}^{\left(n\right)},t=1\dots T,$$

(10)

and the emission distribution of the $n\mathrm{th}$ HMM is

$$\begin{array}{c}\hfill \begin{array}{cc}& X_t^{\left(n\right)}|z_t^{\left(n\right)}\sim \mathcal{N}\left({\mathbf{\mu }}_{z_t^n},{\mathsf{\Sigma }}_{z_t^n}\right),\hfill \\ & \left({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k\right)\sim \mathrm{NIW}\left(u_0,{\lambda }_0,{\mathsf{\Phi }}_0,{\upsilon }_0\right).\hfill \end{array}\end{array}$$

(11)

In order to build a non-parameter model, the hidden states selection matrix F is constructed by a BP−BeP.

$$\begin{array}{c}\hfill \begin{array}{cc}& B\sim \mathrm{BP}\left(\alpha ,B_0\right),\hfill \\ & f_i|B\sim \mathrm{BeP}\left(B\right).\hfill \end{array}\end{array}$$

(12)

From the perspective of the characteristic of BPs, we can find that the greater the concentration parameter $\alpha$, the sparser the hidden state selection matrix F, and greater $\gamma$ will lead to more hidden features.

Given the above model assumptions, the sampled BP-HMM uses the Gibbs sampling method to train the model and uses the gradient based method to learn the parameters. With the state transition matrix for each trajectory being learned, the average state transition matrix for each class can be calculated by the mean operation. The new test trajectories are classified according to their likelihood probabilities conditional on each class.

4. The Proposed Variational BP-HMM

In this section, we will introduce the proposed variational BP-HMM which has more reasonable assumptions and more efficient inference procedure than the sampled BP-HMM. We first describe key points of our model and present our stick-breaking representation for the BP which allows for variational inference. Then we give the joint distribution of the proposed BP-HMM and the variational inference for the BP-HMM.

4.1. BP-HMM with the Shared Hidden State Space and Class Specific Indicators

As introduced above, the existing sampled BP-HMM can jointly learn the trajectories from different classes by sharing a same hidden state space. It can also automatically determine the available states and the corresponding transition matrices for one trajectory by the introducing state selection matrix F. However, in the sampled BP-HMM, the state transition matrix and initial probabilities are trajectory-specific, and it is not appropriate to perform mean operation on these transition matrices and probabilities to obtain a average matrix and probabilities for each class.

In order to model trajectories from different classes more reasonably, we introduce a shared hidden state space and class-specific indicators. We define a state selection vector ${\mathbf{f}}_c$ for each class which are used to distinguish the differences between classes and define state initial probabilities ${\mathbf{\pi }}_0$ and transition matrix ${\mathbf{\pi }}_j$ for each class which are used to capture the commonness with one class. The transition matrix of the cth class from state j is

$${\mathbf{\pi }}_j^{\left(c\right)}\sim Dir\left([r,r,\dots ,r+\kappa ,r,\cdots ]\odot {\mathbf{f}}_c\right),j>0,$$

(13)

and the initial state probability vector ${\mathbf{\pi }}_0^{\left(c\right)}$ is also related to F,

$${\mathbf{\pi }}_0^{\left(c\right)}\sim Dir\left([r,r,\dots ,r,r,\cdots ]\odot {\mathbf{f}}_c\right).$$

(14)

Similarly to the standard HMM, the latent variable ${\mathbf{z}}^{\left(n\right)}$ for the nth trajectory is a discrete sequence with

$$z_1^{\left(n\right)}\sim {\mathbf{\pi }}_0^{\left(y_n\right)},z_{t+1}^{\left(n\right)}|z_t^{\left(n\right)}\sim {\mathbf{\pi }}_{z_t^{\left(n\right)}}^{\left(y_n\right)},t=1,\dots ,T.$$

(15)

where $y_n$ denotes the class of nth trajectory.

From the way of modeling, the proposed new version of the BP-HMM is different from the sampled BP-HMM [1] which learns an HMM for each trajectory, and it is also different from the traditional HMMs which learn an HMM for each class separately. The proposed BP-HMM can use all the sequences from different classes to jointly train a whole BP-HMM with each HMM corresponding to one class. Therefore, the proposed BP-HMM can better model the trajectories from multiple classes and can further make better classification.

4.2. A Simpler Representation for Beta Process

Besides the model assumption, the proposed variational BP-HMM has different representation of the BP. As introduced in Section 2, the IBP construction of the BP describes the process by conditional distributions. This kind of representation is only suitable for sampling methods which are similar to the Chinese restaurant construction of DPs. Therefore, different from the sampled BP-HMM which uses the IBP construction for the BP to lend it to a Gibbs sampler, we use the stick-breaking construction for the BP to adapt to variational inference. There is some work in constructing stick-breaking representation of BPs for variational inference. The stick-breaking construction is used for the IBP which is closely related to the BP and can be seen as a one-parameter BP [26]. The two-parameter BP is also constructed through stick-breaking processes to server for variational inference [29]. Recently, a simpler representation of the two-parameter BP based on stick-breaking construction is developed to make simpler variational inference [33]. In order to approximate posterior inference to the BP with variational Bayesian method more easily, we refer to the simpler representation of the BP [33]. Let $d_k$ mark the round in which the $k\mathrm{th}$ atom appears. That is,

$$\begin{array}{c}\hfill \begin{array}{cc}& d_k=1+\sum _{i=1}^{\infty }\delta \left(\sum _{j=1}^i{C}_j<k\right).\hfill \end{array}\end{array}$$

(16)

Note $\delta \left(·\right)$ is a binary indicator and it equals to 1 if the formula is true. Using the latent indicators, the representation of B in (6) is simplified as

$$\begin{array}{c}\hfill \begin{array}{cc}& B=\sum _{k=1}^{\infty }{V}_{k_{d_k}}\prod _{l=1}^{d_k-1}\left(1-V_{k_l}\right){\delta }_{{\omega }_k},\hfill \end{array}\end{array}$$

(17)

with $\omega$ and V drawn as before.

Let $T_k=-{\sum }_{l<d_k}ln\left(1-V_{k_l}\right)$. Since each individual term $-\ln \left(1-V_{k_l}\right)\stackrel{iid}{\sim }\mathrm{Exponential}\left(\alpha \right)$, it follows that $T_k\stackrel{iid}{\sim }\mathrm{Gamma}\left(d_k-1,\alpha \right)$. This gives the following representations of the BP,

$$\begin{array}{c}\hfill \begin{array}{cc}& B=\sum _{k=1}^{\infty }{V}_k{e}^{-T_k}{\delta }_{{\omega }_k},\hfill \\ & V_k\stackrel{i.i.d.}{\sim }\mathrm{Beta}\left(1,\alpha \right),\hfill \\ & T_k\sim \mathrm{Gamma}\left(d_k-1,\alpha \right),\hfill \\ & \sum _{k=1}^{\infty }1\left(d_k=r\right)\stackrel{i.i.d.}{\sim }\mathrm{Poisson}\left(\gamma \right),r\in {\mathbb{N}}^{+},\hfill \\ & {\omega }_k\stackrel{i.i.d.}{\sim }\frac{B_0}{\gamma }.\hfill \end{array}\end{array}$$

(18)

Here we should notice that each $d_k$ does not have a distribution, but the cardinality of $\{d_k=r\}$ is drawn by $\mathrm{Poisson}\left(\gamma \right)$. In addition, $T_k=0$ with probability one when $d_k=1$. In this BP, the atom ${\omega }_k=\{{\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k\}$ and Gamma priors with hyper-parameters $\{a_1,a_2\}$, $\{b_1,b_2\}$ are given to $\alpha$ and $\gamma$:

$$\begin{array}{c}\hfill \begin{array}{cc}& \alpha \sim \mathrm{Gamma}\left(a_1,a_2\right),\hfill \\ & \gamma \sim \mathrm{Gamma}\left(b_1,b_2\right).\hfill \end{array}\end{array}$$

(19)

4.3. Joint Distribution of the Proposed BP-HMM

Assume that the total class number is C and the trajectory number is N. Let X represent the data, $W=\{\alpha$, $\gamma ,$$\{{\mathbf{\mu }}_k$, ${\mathsf{\Sigma }}_k\}$, $\left\{d_k\right\}$, $\left\{V_k\right\}$, $\left\{T_k\right\}$, $\left\{f_{ck}\right\}$, $\left\{{\mathbf{\pi }}_k^{\left(c\right)}\right\}$, $Z\}$ represents the set of all latent variables in the model, including $\theta$ which is the set of all the hyper-parameters, and Y is the set of all the class labels.

The probabilistic graphical model is shown in Figure 2, where its joint likelihood is

$$\begin{array}{c}\hfill \begin{array}{cc}& p\left(X,W|\theta \right)=p\left(X|W,\theta \right)\times p\left(W|\theta \right).\hfill \end{array}\end{array}$$

(20)

The likelihood $p\left(X|W,\theta \right)$ is defined as a multi-normal distribution by

$$p\left(X|W,\theta \right)=\prod _{n=1}^N\prod _{t=1}^T\prod _{k=1}^K\mathcal{N}{\left({\mathbf{x}}_t^{\left(n\right)}\right|{\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k)}^{\delta (z_t=k)}.$$

(21)

The prior distribution of the parameter W and detailed setup are expressed in Appendix A.

4.4. Variational Inference for the Proposed BP-HMM

We use a factorized variational distribution over all the latent variables to approximate the intractable posterior $p\left(W|X,\theta \right)$. Two truncations are set in the inference: one is truncation of the number of hidden states at K and the other is the truncation of the round number at R. Specifically, we assume the variational distribution as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Q=& q\left(\alpha \right)q\left(\gamma \right)\prod _{k=1}^K\{q\left({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k\right)q\left(d_k\right)q\left(V_k\right)q\left(T_k\right)\hfill \\ & \left.\times \prod _{c=1}^{C}q\left(f_{ck}\right)q\left({\mathbf{\pi }}_k^{\left(c\right)}\right)\right\}\prod _{c=1}^{C}q\left({\mathbf{\pi }}_0^{\left(c\right)}\right)\prod _{n=1}^{N}q\left({\mathbf{z}}^{\left(n\right)}\right),\hfill \end{array}\end{array}$$

(22)

where

$$\begin{array}{ccc}& & q\left(\alpha \right)=\mathrm{Gamma}\left(\alpha |k_1,k_2\right),\hfill \\ & & q\left(\gamma \right)=\mathrm{Gamma}\left(\gamma |{\tau }_1,{\tau }_2\right),\hfill \\ & & q\left({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k\right)=\mathrm{NIW}\left({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k|u_k,{\lambda }_k,{\mathsf{\Phi }}_k,{\upsilon }_k\right),\hfill \\ & & q\left(d_k\right)=\mathrm{MultiNomial}\left(d_k|{\mathbf{\phi }}_k\right),\hfill \\ & & q\left(V_k\right)=\mathrm{Beta}\left(V_k|{\tau }_{k_1},{\tau }_{k_2}\right),\hfill \\ & & q\left(T_k\right)=\mathrm{Gamma}\left(T_k|u_k^{\prime },v_k^{\prime }\right),\hfill \\ & & q\left(f_{ck}\right)=\mathrm{Bernoulli}\left(f_{ck}|{\upsilon }_{ck}\right),\hfill \\ & & q\left({\mathbf{\pi }}_k^{\left(c\right)}\right)=\mathrm{Dir}\left({\mathbf{\pi }}_k^{\left(c\right)}|{r_{k1}^{\prime }}^{\left(c\right)},{r_{k2}^{\prime }}^{\left(c\right)},\cdots ,{r_{kK}^{\prime }}^{\left(c\right)}\right),\hfill \\ & & q\left({\mathbf{z}}^{\left(n\right)}|y_n\right)=\prod _{t=1}^T\prod _{k_1=1}^K\prod _{k_2=1}^K{a}_{k_1{k}_2}^{{\left(y_n\right)}^{\delta \left(z_t^{\left(n\right)}=k_1,z_{t+1}^{\left(n\right)}=k_2\right)}}\hfill \\ & & \times \prod _{k=1}^K{a_{0k}^{\left(y_n\right)}}^{\delta \left(z_0^{\left(n\right)}=k\right)}\prod _{t=1}^T\prod _{k=1}^K{b}_{tk}^{{\left(y_n\right)}^{\delta \left(Z_t^{\left(n\right)}=k\right)}}.\hfill \end{array}$$

It is obvious that $V_k$ and $T_k$ do not have conjugate posterior. Thus the distributions are selected for better accuracy and more convenience. Here $a_{0k}^{*}$ is an estimation of the probability of the initial state distribution, $a_{j_{1j_2}}^{*}$, where $j_1>0$ and $j_2>0$ is an estimation of the probability of transition from state $j_1$ to $j_2$ and $b^{*}{t}_j$ is an estimation of the emission probability density given the system in state j at time point t. In order to simplify our representation, we do not use sub-index. Here $a_i=\left\{a_{ij}\right\}$, $j=1,\dots ,K$. Let $\varphi$ be the set of variational parameters. We expand the lower bound as $\mathcal{L}\left(X,\varphi \right)={\mathbb{E}}_Q\left(\ln P\left(X,W|\theta \right)\right)-{\mathbb{E}}_Q\left[\ln Q\right]$ which is expressed in detail in Appendix B.

4.5. Parameter Update

In the framework of variational mean field approximation, the parameters of some variational distributions can be analytically solved using

$$\ln q\left({\mathbf{w}}_j\right)={\mathbb{E}}_{q(W\ne {\mathbf{w}}_j)}\left[\ln p\left(X,W|\theta \right)\right]+\mathrm{const}.$$

(23)

However, in some cases that the prior distribution and posterior distribution over one latent variable are not conjugate, the variational distribution over this variable cannot have an analytical solution. The parameters of this variational distribution should be optimized through gradient based methods with the variational lower bound being the objective.

In our model, the variational distributions $q\left(\alpha \right)$, $q\left(\gamma \right)$, $q({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k)$, $q\left(d_k\right)$, $q\left({\mathbf{\pi }}_k\right)$, $q\left(Z\right)$ have a closed form solution, and we can get their parameter update formulas according to (23). While the variational distributions $q\left(V_k\right)$, $q\left(T_k\right)$, $q\left(f_{ck}\right)$ cannot be analytically solved, we can update their parameters by corresponding gradients. Next, we give the way of calculating variational distributions and show the procedure for training the variational BP-HMM in Algorithm 1. The detailed parameters update formulas or the gradients with respect to the parameters are presented in Appendix C.

Algorithm 1 Variational Inference for the Proposed BP-HMM

1: Initialize $\theta$ and $\varphi$.   2: Given R and threshold and Initialize RunTime = 0;   3: while$|\mathcal{L}-{\mathcal{L}}_{\mathbf{old}}|<\mathrm{threshold}$ or RunTime < R do   4:  ${\mathcal{L}}_{\mathbf{old}}=\mathcal{L}$   5:  for each trajectory n do   6:   Update $q\left({\mathbf{z}}^{\left(n\right)}\right)$   7:   Calculate $q(z_t^{\left(n\right)}=k)$ and $q(z_t^{\left(n\right)}=k_1,z_{t+1}^{\left(n\right)}=k_2)$   8:  end for   9:  for each class c do 10:   Update each $q\left({\mathbf{\pi }}_k^{\left(c\right)}\right)$, $k=0,\dots ,K$ 11:   Update each $q\left(f_{ck}\right)$, $k=1,\dots ,K$ 12:  end for 13:  for each $k=1,\dots ,K$ do 14:   Update $q({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k)$, $q\left(d_k\right)$, $q\left(T_k\right)$, $q\left(V_k\right)$ 15:  end for 16:  Update $q\left(\alpha \right)$, $q\left(\gamma \right)$  Calculate $\mathcal{L}$ 17: end while

4.5.1. Calculation for $q\left(\alpha \right)$, $q\left(\gamma \right)$, $q({\mathit{\mu }}_k,{\mathsf{\Sigma }}_k)$, $q\left(d_k\right)$, $q\left({\mathit{\pi }}_k^{\left(c\right)}\right)$, $q\left(Z\right)$

$$\begin{array}{ccc}\hfill & & lnq\left(\alpha \right)={\mathbb{E}}_q[lnp\left(\alpha \right)+\sum _{k=1}^{K}lnp\left(V_k\right|\alpha )+lnp\left(T_k\right|d_k,\alpha )],\hfill \\ & & lnq\left(\gamma \right)={\mathbb{E}}_q[lnp\left(\gamma \right)+\sum _{k=1}^{K}lnp\left(d_k\right|\gamma )],\hfill \\ & & lnq({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k)={\mathbb{E}}_q[lnp({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k|\theta )\hfill \\ & & +\sum _{n=1}^N\sum _{t=1}^{T}p({\mathbf{x}}^{\left(n\right)}|{\mathbf{z}}^{\left(n\right)},{\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k)],\hfill \\ & & lnq\left(d_k\right)={\mathbb{E}}_q[lnp\left(d\right|\gamma )+lnp\left(T_k\right|d_k,\alpha )\hfill \\ & & +\sum _{c=1}^{C}lnp\left(f_{ck}\right|V_k,T_k,d_k)],\hfill \\ & & lnq({\mathbf{\pi }}_k^{\left(c\right)})={\mathbb{E}}_q[lnp({\mathbf{\pi }}_k^{\left(c\right)}|f_{ck},r,\kappa )\hfill \\ & & +\sum _{n=1}^N\sum _{t=1}^{T-1}\delta (y_n=c)lnp(z_{t+1}^{\left(n\right)}|{\mathbf{\pi }}_k^{\left(c\right)},z_t^{\left(n\right)}=k)],\hfill \\ & & lnq({\mathbf{z}}^{\left(n\right)})={\mathbb{E}}_q[lnp({\mathbf{x}}^{\left(n\right)}|{\mathbf{z}}^{\left(n\right)},{\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k)\hfill \\ & & +\sum _{n=1}^{N}lnp({\mathbf{z}}^{\left(n\right)}|{\mathsf{\Pi }}^{\left(y_n\right)})],\hfill \end{array}$$

4.5.2. Optimization for $q\left(V_k\right)$, $q\left(T_k\right)$, $q\left(f_{ck}\right)$

The variational parameters of $q\left(V_k\right)$, $q\left(T_k\right)$, $q\left(f_{ck}\right)$ include $\{{\tau }_{k_1},{\tau }_{k_2}\}$, $\{u_k^{\prime },v_k^{\prime }\}$, $\left\{{\upsilon }_{ck}\right\}$. They are updated by the gradient based method where the gradients of the lower bound $\mathcal{L}$ with respect to these parameters should be calculated.

4.5.3. Remarks

Note that when updating the BP parameters, we should calculate the expectation as

$$\begin{array}{ccc}\hfill {\mathbb{E}}_q[lnp\left(f_{ck}\right|V_k,T_k)]& =& v_{ck}{\mathbb{E}}_q[ln{V}_k{e}^{-T_k}]\hfill \\ & & +(1-v_{ck}){\mathbb{E}}_q[1-ln{V}_k{e}^{-T_k}],\hfill \end{array}$$

of which the second term is intractable. We refer the work in [33] to use a Taylor expansion to ${\mathbb{E}}_q[ln(1-V_k{e}^{-T_k})]$ about the point one,

$${\mathbb{E}}_q[ln(1-V_k{e}^{-T_k})]=-\sum _{m=1}^M\frac{1}{m}{\left(V_k{e}^{-T_k}\right)}^m.$$

(24)

For clarity, we define each term $\frac{1}{m}\mathbb{E}\left[{\left(V_k{e}^{-T_k}\right)}^m\right]$ in the Taylor expansion using the notation ${\Delta }_k\left(m\right)$ as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\Delta }_k\left(m\right)& =\frac{1}{m}\frac{\mathsf{\Gamma }\left({\tau }_{k_1}+{\tau }_{k_2}\right)}{\mathsf{\Gamma }\left({\tau }_{k_1}+{\tau }_{k_2}+m\right)}\frac{\mathsf{\Gamma }\left({\tau }_{k_1}+m\right)}{\mathsf{\Gamma }\left({\tau }_{k_1}\right)}{\left(\frac{v_k^{\prime }}{v_k^{\prime }+m}\right)}^{u_k^{\prime }}\hfill \\ & =\prod _{i=1}^m\frac{{\tau }_{k_1}+i-1}{{\tau }_{k_1}+{\tau }_{k_2}+i-1}{\left(\frac{v_k^{\prime }}{v_k^{\prime }+m}\right)}^{u_k^{\prime }},\hfill \end{array}\end{array}$$

and define ${\Delta }_k\left(·\right)={\sum }_{m=1}^M{\Delta }_k\left(m\right)$. Therefore, ${\mathbb{E}}_q[ln(1-V_k{e}^{-T_k})]=-{\Delta }_k\left(·\right)$.

4.6. Classification

Our model is applicable to trajectory recognition like human activity trajectory recognition. We use the proposed variational BP-HMM to model all the training data from different classes, with each HMM corresponding to a class. Given the learned model with the hyperparameters and variational parameters $\{\theta ,\varphi \}$, a new test trajectory ${\mathbf{x}}^{*}$ can be classified according to its marginal likelihood $p\left({\mathbf{x}}^{*}\right|\theta ,\varphi )$. Denote $y^{*}$ as the label of the test trajectory; the classification criteria can be expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y^{*}& =\underset{c}{argmax}lnp\left({\mathbf{x}}^{*}\right|\{{{\mathbf{a}}_k^{*}}^{\left(c\right)},u_k,{\lambda }_k,{\upsilon }_k,{\mathsf{\Phi }}_k\},{{\mathbf{a}}_0^{*}}^{\left(c\right)}),\hfill \\ & =\underset{c}{argmax}ln(\int p\left({\mathbf{x}}^{*}\right|\{{\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k\},\mathbf{z})p\left(\mathbf{z}\right|\left\{{{\mathbf{a}}_k^{*}}^{\left(c\right)}\right\})\hfill \\ & \prod _{k=1}^{K}p({\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k|u_k,{\lambda }_k,{\upsilon }_k,{\mathsf{\Phi }}_k)d\mathbf{z}d{\mathbf{\mu }}_{k}d{\mathsf{\Sigma }}_k),\hfill \end{array}\end{array}$$

(25)

where ${a_{jk}^{*}}^{\left(c\right)}$ is an estimate of the probability of transition from state j to k in the cth class. The likelihood can be calculated through the forward-backward algorithm.

This classification mechanism is more reasonable than the method in [1], as the transition matrix is actually learned.

5. Experiment

To demonstrate the effectiveness of our model on trajectory recognition, we conduct experiments on one synthetic dataset and two real-world datasets; the detailed data statistics are illustrated in

Table 1. Data statistics for the CCP, HATR and WFNT datasets and corresponding classes.

Datasets	#Train Trajectories	#Classes (Descriptions)
CCP	20 (5/class)	4 (Normal, Cyclic, IT, DT)
HATR	300 (50/class)	6 (PTSS, PTES, GA, CPH, WFB, WTS)
WFNT	40 (10/class)	4 (F, L, R, B)

and the following subsections. We compare our model with HCRF, LSTM [34], HMM-BIC and the sampled BP-HMM. In particular, in HCRF, the number of hidden states is set to 15 and the size of the window is set to 0. In LSTM, we use a recurrent neural network with one hidden layer as its architecture. In HMM-BIC, the state number is selected from the range $[1,20]$. In the sampled BP-HMM, the hyperparameters are set according to Sun et al. [1]. In the variational BP-HMM, the hyperparameters $\{a_1,a_2,b_1,b_2,r,\kappa \}$ are randomly initialized and selected by maximizing the variational lower bound, and the emission hyperparameters are initialized with k-means. Particularly, the state truncation parameters in variational BP-HMM are set according to specific datasets, e.g., $K=7$ for the synthetic data and $K=20$ for the two real-world data. All experiments are repeated ten times with different training and test division methods, and the average classification accuracy with the standard deviation is reported.

5.1. Synthetic Data

The synthetic data called control chart patterns (CCP) have some quantifiable similarities. They contain four pattern types which can be downloaded from the UCI machine learning repository. The CCP are trajectories that show the level of a machine parameter plotted against time. 400 trajectories are artificially generated by the following four equations [35]:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & 1.\mathrm{Normal}\mathrm{pattern}\left(\mathrm{Normal}\right):y\left(t\right)=m+rs,\hfill \\ & \mathrm{where}m=3,s=2\mathrm{and}0<r<1.\hfill \\ & 2.\mathrm{Cyclic}\mathrm{pattern}\left(\mathrm{Cyclic}\right):y\left(t\right)=m+rs+a\mathrm{Sin}\left(2\pi t/T\right),\hfill \\ & \mathrm{where}0<a,T<15.\hfill \\ & 3.\mathrm{Increasing}\mathrm{trend}\left(\mathrm{IT}\right):y\left(t\right)=m+rs+gt,\hfill \\ & \mathrm{where}0.2<g<0.5.\hfill \\ & 4.\mathrm{Decreasing}\mathrm{trend}\left(\mathrm{DT}\right):y\left(t\right)=m+rs-gt,\hfill \\ & \mathrm{where}0.2<g<0.5.\hfill \end{array}\end{array}$$

Figure 3 shows the generated synthetic data. In this experiment, 20 trajectories are used for training with 5 trajectories for each class. The classification results are presented in

Table 2. Comparisons of the classification accuracy for the proposed method VBP-HMM versus HCRF, LSTM, HMM-BIC and the sampled BP-HMM in CCP.

	Classification Accuracy
Approach	HCRF	LSTM	HMM-BIC	SBP-HMM	VBP-HMM
CCP	0.88 ± 0.03	0.95 ± 0.02	0.97 ± 0.01	0.96 ± 0.02	1.00 ± 0.00

. The results are obtained through five-fold cross-validation. In order to illustrate that the sharing patterns have been learned by our method, the Hinton diagrams of the variational parameter V are given in Figure 4, where the occurrence probabilities of the hidden states are presented by the sizes of the blocks. For example, we can find that IT and DT share the 4th, 5th, 6th features.

We compare our method with HCRF, LSTM, HMM-BIC and the sampled BP-HMM. As we can see from Table 2, our method outperforms all the other methods.

In this experiment, the sharing patterns contribute to improving the performance. Since an HMM is created for each class of trajectories in our proposed method instead of each trajectory in the sampled BP-HMM, our method has better performance than the sampled BP-HMM.

5.2. Human Activity Trajectory Recognition

Human activity trajectory recognition (HATR) [36] is important in many applications such as health care. In our human activity trajectory recognition experiment, parking lot data are collected from the video [1]. We use the data tagged manually [1], which has 300 trajectories with 50 trajectories for each class. Six classes are defined, which are “passing through south street” (PTSS), “passing through east street” (PTES), “going around” (GA), “crossing park horizontally” (CPH), “wandering in front of building” (WFB) and “walking in top street” (WTS). As seen from [1], the sampled BP-HMM is the best method among the methods including HCRF, LSTM, HMM-BIC and the sampled BP-HMM in HATR. Here we use the same training and test data to compare the variational BP-HMM with the sampled BP-HMM.

Table 3. Comparisons of the classification accuracy for the proposed method VBP-HMM versus HCRF, LSTM, HMM-BIC and the sampled BP-HMM in HATR.

	Classification Accuracy
Approach	HCRF	LSTM	HMM-BIC	SBP-HMM	VBP-HMM
HATR	0.68 ± 0.03	0.75 ± 0.03	0.95 ± 0.02	0.91 ± 0.02	0.96 ± 0.02

shows the comparisons of the classification accuracy for the proposed method VBP-HMM versus HCRF, LSTM, HMM-BIC and the sampled-BP-HMMs in HATR. The results are obtained through five-fold cross-validation. As can be seen from Table 3, the accuracy of our method is $0.96$, while the accuracy of the sampled BP-HMM is $0.91$ [1]. The detailed confusion matrix for our method is given in

Table 4. Classification accuracy for human activity trajectory recognition.

	Classification Accuracy
Predicted Class	PTSS	PTES	GA	CPH	WFB	WTS
PTSS	1.00	0.00	0.00	0.00	0.00	0.00
PTES	0.00	1.00	0.00	0.00	0.00	0.00
GA	0.00	0.00	0.97	0.00	0.00	0.00
CPH	0.00	0.00	0.03	0.96	0.00	0.00
WFB	0.00	0.00	0.00	0.03	1.00	0.23
WTS	0.00	0.00	0.00	0.01	0.00	0.77

. The state sharing patterns learned by variational BP-HMM are displayed with the Hinton diagrams in Figure 5, in which GA and CPH, as well as GA and WTS, are more likely to share states. The good performance verifies the superiority of modeling an HMM for each class. Moreover, we take some examples of the correct classification and misclassification results for visualization as in Figure 6 and Figure 7. As illustrated in Figure 7, the misclassified trajectories often contain some deceptive subpatterns such as the trajectory of CPH in subfigure (d) containing a back turn and a left turn like the GA class.

5.3. Wall-Following Navigation Task

We perform the Wall-Following navigation task (WFNT) in which data are collected from the sensors on the mobile robot SCITOS-G5 [4]. We think that this task is a trajectory with historical data, and two ultrasound sensors datasets are selected, because the cost is as low as possible in civil applications with acceptable accuracy. There are 187 trajectories in the data and four classes need to be recognized, which are “front distance” (F), “left distance” (L), “right distance” (R) and “back distance” (B). We randomly select 40 training trajectories with 10 for each class. The confusion matrix of classification is shown in

Table 5. The confusion matrix of the Wall-Following navigation recognition.

	Classification Accuracy
Predicted Class	F	L	R	B
F	0.95	0.18	0.00	0.00
L	0.02	0.73	0.00	0.00
R	0.03	0.09	0.95	0.0
B	0.00	0.00	0.05	1.00

and the state sharing patterns learned by variational BP-HMM are displayed with the Hinton diagrams in Figure 8, where R and F, as well as R and B, have a small number of shared states.

The comparison of the classification accuracy for our method VBP-HMM versus HCRF, LSTM, HMM-BIC and the sampled BP-HMM is shown in

Table 6. Comparisons of the classification accuracy for the proposed method VBP-HMM versus HCRF, LSTM, HMM-BIC and the sampled BP-HMM in WFNT.

	Classification Accuracy
Approach	HCRF	LSTM	HMM-BIC	SBP-HMM	VBP-HMM
WFNT	0.80 ± 0.03	0.73 ± 0.08	0.86 ± 0.04	0.85 ± 0.02	0.89 ± 0.01

. The results are obtained by five-fold cross-validation. It is obvious that our method is much better than the sampled BP-HMM, because we create an HMM for each class of trajectories rather than create an HMM for each trajectory. Although the sharing patterns are not obvious in this experiment, our method has better performance than the other methods. As we have analyzed, sharing patterns among different classes will be learned automatically by our model, which helps to localize precisely the difference of different classes. When there is no sharing pattern among classes, the advantage will be weakened.

5.4. Performance Analysis

In our experiments, the results show that the proposed variational BP-HMM has a great improvement compared to the sampled BP-HMM which uses average transition over trajectories from each class. We analyze the advantages of variational BP-HMM for the following reasons. Due to the small amount of training data in our experiment, the performance of LSTM is not satisfactory. HMM-BIC finds an optimal state number through model selection but it cannot make use of the shared information among classes, and its performance is the second-best overall. Although the sample BP-HMM can share hidden states among classes, it does not make correct use of the shared information in classification and thus does not gain better results. Our proposed variational BP-HMM constructs a mechanism to learn shared hidden states by introducing state indicator variables and maintains class-specific state transition matrices which are very helpful for classification tasks.

Moreover, we give the total cost time of the variational BP-HMM, HMM-BIC, LSTM, HCRF and the sampled BP-HMM in

Table 7. Comparisons of total time cost for the proposed method VBP-HMM versus HCRF, LSTM, HMM-BIC and the sampled-BP-HMM in experiments.

	Total Time Cost (s)
Approach	HCRF	LSTM	HMM-BIC	SBP-HMM	VBP-HMM
CCP	196	6	54	2151	117
HATR	118	13	93	2521	205
WFNT	1005	8	115	1819	312

, where we can see the variational BP-HMM performs much more efficiently than the sampled BP-HMM. This is attributed to the efficiency of the variational methods. Although the sampled BP-HMM and the variational BP-HMM have similar time complexity, due to the sampling operation, the cost time of the sampled BP-HMM is usually several times that of the variational BP-HMM. In other words, the variational BP-HMM converges much faster than the sampled BP-HMM. Besides, compared with HMM-BIC, it only takes about twice the time to achieve significant performance improvements. Above all, we can conclude that the proposed variational BP-HMM is an effective and efficient method for trajectory recognition.

6. Conclusions

In this paper, we have proposed a novel variational BP-HMM for modeling and recognizing trajectories. The proposed variational BP-HMM has shared hidden state space which is used to capture the commonality of the cross-category data and class-specific indicators which are used to distinguish the data from different classes. As a result, in the variational BP-HMM, multiple HMMs are used to model multiple classes of trajectories among which a hidden state space is shared.

The more reasonable assumptions of the proposed model make it more suitable for jointly modeling trajectories over all classes and further making trajectory recognition. Experimental results both on synthetic and real-world data have verified that the proposed variational BP-HMM can find the feature sharing patterns among different classes, which helps to better model trajectories and further improve the classification performance. Moreover, compared with the sampled BP-HMM, the derived variational inference for the proposed BP-HMM can reduce the time cost of the training procedure. The experimental time records also show the efficiency of the proposed variational BP-HMM.