{'fraction_non_alphanumeric': 0.045570752908091525, 'fraction_numerical': 0.015435255017256807, 'mean_word_length': 4.382352941176471, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Disentangled sticky hierarchical Dirichlet process hidden Markov model

The Hierarchical Dirichlet Process Hidden Markov Model (HDP-HMM) has been used widely as a natural Bayesian nonparametric extension of the classical Hidden Markov Model for learning from sequential and time-series data. A sticky extension of the HDP-HMM has been proposed to strengthen the self-persistence probability in the HDP-HMM. However, the sticky HDP-HMM entangles the strength of the self-persistence prior and transition prior together, limiting its expressiveness. Here, we propose a more general model: the disentangled sticky HDP-HMM (DS-HDP-HMM). We develop novel Gibbs sampling algorithms for efficient inference in this model. We show that the disentangled sticky HDP-HMM outperforms the sticky HDP-HMM and HDP-HMM on both synthetic and real data, and apply the new approach to analyze neural data and segment behavioral video data.


Introduction
Hidden Markov models (HMMs) provide a powerful set of tools for modeling time series data. In the HMM we assume that the time series observations are modulated by underlying latent time-varying variables which take a discrete set of states. This model class is useful in its own right and can also be incorporated as a building block for more complicated models. It has been widely used in speech recognition [10], musical audio analysis [26,14], acoustic-phonetic modeling [18], behavior segmentation [9,31,1], sequential text modeling [34,13], financial time series data analysis [9,33], computational biology [17], and many other fields.
Selecting the number of HMM states is an important question for practitioners. Classical model selection techniques can be used, but these methods can be computationally intensive and are sometimes unreliable in practice [3,29]. Also, for real datasets, it is often reasonable to assume that the number of latent states may be unbounded, violating classical assumptions needed to establish consistency results for model selection. Based on previous work in [2], [34] proposed the Hierarchical Dirichlet process HMM (HDP-HMM), a Bayesian arXiv:2004.03019v1 [stat.ML] 6 Apr 2020 nonparametric framework. In the HDP-HMM the transition matrix follows a hierarchical Dirichlet process (HDP) prior. [10] noted that the HDP-HMM tends to rapidly switch among redundant states, and proposed the sticky HDP-HMM (S-HDP-HMM), which strengthens the self-persistence probability. This modification often leads to significant improvements in modeling real data.
In the HMM, it is important to distinguish three features: 1, the similarity of the rows of the transition matrix; 2, the average self-persistence probability of the latent states (controlled by the mean of the diagonal of the transition matrix); and 3, the strength of the self-persistence prior across states (i.e., the inverse prior variance of the diagonal elements of the transition matrix). In the HDP-HMM, there is only one parameter controlling feature 1. The sticky HDP-HMM adds one more parameter to control feature 2, but still entangles features 1 and 3 with only one parameter, thus limiting the expressiveness of the prior.
We show that we can add one additional hyperparameter to generalize the sticky HDP-HMM formulation, obtaining three degrees of freedom to model the three features discussed above. We call this new model the disentangled sticky HDP-HMM (DS-HDP-HMM).
The rest of the paper is organized as follows. In section 2, we provide a brief introduction to Bayesian HMM, HDP-HMM, and sticky HDP-HMM. In section 3, we discuss the limitations of these models. In section 4, we introduce disentangled sticky HDP-HMM, and in section 5 we develop efficient Gibbs sampling inference methods for this new model. Section 6 demonstrates the effectiveness of the disentangled sticky HDP-HMM on both synthetic and real data, including applications to analyzing neural data and segmenting behavior video. A table of notation can be found in the supplementary material.

Background on Bayesian HMM and HDP-HMM
Our goal here is to fit an HMM to time series data. On its face, this would seem to be a solved problem; after all, we can compute the HMM likelihood easily, and the basic expectation-maximization algorithm for HMM fitting is textbook material [30]. Nonetheless, a fully Bayesian solution to this problem has remained elusive. Specifically, we would like to be able to compute a posterior over all of the unknown HMM parameters (including the number of latent states). Quantification of posterior uncertainty is critical in many applications: for example, given short time series data, often we do not have enough data to sufficiently identify the HMM parameters. Even for longer time series data, we might want to fit richer models as we collect more data. Here "richer models" correspond to more latent states, and since the number of parameters in the HMM grows quadratically with the number of states, we may be left again with some irreducible uncertainty about the model parameters.
The HDP-HMM [34] provides a useful starting point for fully Bayesian HMM inference. The basic idea here is to sample a global transition distribution prior from a Dirichlet process (described below), and then for each latent state we sample a transition distribution from this shared (random) global prior distri-bution. To develop the details of this HDP-HMM idea we first need to define some notation for the Dirichlet process (DP). Given a base distribution H on a parameter space Θ and a positive concentration parameter γ, a Dirichlet process G ∼ DP(γ, H) (sometimes also denoted by DP(γH)) can be constructed by the following stick-breaking procedure [32]: let where β ∼ GEM(γ) is a random probability mass function (p.m.f.) defined on a countably infinite set as follows: Then the discrete random measure G = i β i δ θi is a sample from DP(γH), where δ θi denotes the Dirac measure centered on θ i .
The HDP-HMM [34] uses the DP to define a prior on the rows of the HMM transition matrix in a setting where the number of latent states is unbounded. The HDP-HMM is defined as Transition matrix prior : π j iid ∼ DP(αβ), j = 1, 2, · · · Latent states : z t ∼ π zt−1 , t = 1, · · · , T Observations : y t ∼ f (y|θ zt ), t = 1, · · · , T Here, β and {θ j } ∞ j=1 are defined as in the DP described above, and then each transition distribution π j for state j is defined as a random sample of a second DP with base measure β and concentration parameter α. Here α controls how similar π j is to the global transition distribution β. Finally, as usual, z t denotes the state of a Markov chain at time t, and the observation y t is independently distributed given the latent state z t and parameters {θ j } ∞ j=1 , with emission distribution f (·). The sticky HDP-HMM from [10] modifies the transition matrix prior by adding a point mass distribution with stickiness parameter κ to encourage selfpersistence: where δ j denotes the Dirac measure centered on j. Figure 1(a) provides the graphical model for the sticky HDP-HMM.

Limitations of the HDP-HMM and sticky HDP-HMM
The HDP-HMM uses the concentration parameter α to control the following feature of the HMM: feature 1: the strength of the transition matrix prior, or the similarity of the rows of the transition matrix.
In other words, a large value of α here means that the transition probability for each state is close to the global transition distribution β.
A flexible model should have additional parameters to control two additional features: feature 2: the average self-persistence probability, or the mean of the diagonal of the transition matrix. feature 3: the strength of the self-persistence prior, or the similarity of the diagonal elements of the transition matrix.
The sticky HDP-HMM adds just one parameter κ compared to the HDP-HMM. feature 2 is controlled by κ/(α + κ), while both feature 1 and feature 3 are controlled by α + κ. Note that a strong prior (large α + κ) means that both the self-persistence probability and the transition probability are quite similar across states, and it is impossible to separate the strength of these two elements using this parameterization. These three features should occupy three degrees of freedom in total, but the formulation of sticky HDP-HMM is only able to traverse a two-dimensional sub-manifold of this three-dimensional space, limiting the expressiveness of the sticky HDP-HMM prior.
More concretely, consider the speaker diarization example studied in [10]. The task is to distinguish the speakers in an audio recording of a conversation; the current speaker is the hidden state. Suppose that some speakers are likely to speak for a very long time while others are terse (implying a small feature 3, i.e. small α+κ), but the identity of the next speaker is independent of the identity of the previous speaker (implying a big feature 1, i.e. big α + κ). The sticky HDP-HMM would have trouble in this scenario. The opposite case is plausible as well: as an example, consider a group of people sitting in a circle and expressing their opinions in a clock-wise fashion (implying a small feature 1 since the distribution of the next speaker is highly dependent on the previous speaker); if each speaker talks for a very similar amount of time (implying a big feature 3), then the sticky HDP-HMM would have difficulty with this scenario as well.

Disentangled sticky HDP-HMM
Now that we have diagnosed this lack of flexibility in the HDP-HMM and sticky HDP-HMM, we can construct a new more flexible model that separates the strength of the self-persistence from the similarity of the transition probabilities. Specifically, we modify the transition matrix prior as where the transition distribution π j is a mixture distribution. A sample from π j has the self-persistence probability of κ j to come from a point mass distribution at j, and has probability 1 − κ j to come fromπ j , a sample from the DP with base measure β. We call this new model the disentangled sticky HDP-HMM.
Here the beta(ρ 1 , ρ 2 ) prior has the flexibility to control both the expectation of self-persistence (feature 2), and the variability of self-persistence (feature 3). Meanwhile α is free to control the variability of the transition probability around the mean transition β (feature 1). In short, we use 3 parameters (ρ 1 , ρ 2 , α) rather than the 2 parameters (κ, α) to separate the strength of the self-persistence and the transition priors.
When ρ 1 = 0 and ρ 2 > 0, all the κ j = 0, and the disentangled sticky HDP-HMM reduces to HDP-HMM. Importantly, the sticky HDP-HMM is also a special case of the disentangled sticky HDP-HMM, as shown in Theorem 1. See the supplementary material for a proof. An equivalent formulation of z t ∼ π zt−1 in the disentangled sticky HDP-HMM is as follows: where we add binary auxiliary variables w t , which decide whether the next step is self-persistent or switching. See Figure 1(b) for the corresponding graphical model. We use this formulation to facilitate inference in section 5.

Gibbs sampling inference
In this section, we introduce a new direct assignment Gibbs sampler (Algorithm 1) (similar to [34]) and a new weak-limit Gibbs sampler (Algorithm 2) (similar to [10]) for inference in the disentangled sticky HDP-HMM. The direct assignment sampler generates samples from the true posterior of the disentangled sticky HDP-HMM when the Gibbs chains converge. The weak-limit sampler uses finite approximation of the HDP-HMM to accelerate the mixing rate of the Gibbs chains and can be easily adapted to parallel computing. [9] noted that the weak-limit sampler was useful for observation models with dynamics such as auto-regressive HMM (ARHMM) or switching linear dynamic system (SLDS). For detailed derivations of these two algorithms, see the supplementary material.

Direct assignment sampler
The direct assignment sampler for the HDP-HMM marginalizes transition distributions π j and parameters θ j and sequentially samples z t given all the other states z \t , observations {y t } T t=1 , and the global transition distribution β. The main difference between our direct assignment sampler and the corresponding HDP-HMM sampler is that instead of only sampling z t , we sample {z t , w t , w t+1 } in blocks. We sample parameters in the HDP prior only with z t whose corresponding w t = 0.
For step 1, we sequentially compute the probability for each possible case of the posterior p(z t , w t , w t+1 |z \t , w \{t,t+1} , {y t } T t=1 , α, β, {κ j } K+1 j=1 ), and sample {z t , w t , w t+1 } from the corresponding multinomial distribution. If z t = K + 1, i.e. a new state appears, we will increment K, sample self-persistence probability κ K+1 for a new state from the prior, and update β using stick-breaking. For step 2, given w t+1 whose corresponding z t = j, we can sample κ j using beta-binomial conjugacy. For step 3, by introducing auxiliary variables {m jk } K j,k=1 , we sample β using Dirichlet-multinomial conjugacy. For step 4, we compute the empirical transition matrix {n jk } K j,k=1 , where n jk is the number of transitions from state j to k with w t = 0 in {z t } T t=1 , and introduce additional auxiliary variables. Then the posterior of α and γ are gamma-conjugate, given the auxiliary variables. We approximate the posterior of ρ 1 , ρ 2 by finite grids. The complexity for each step in Algorithm 1 is O(T K), O(K), O(K), and O(K) respectively, so the total complexity per iteration is O(T K).
It is worth noting that instead of modeling {κ j } K+1 j=1 as samples from a beta distribution, it is natural to consider any distribution on the [0, 1] interval. The Gibbs algorithm here is easily adaptable to cases where we have extra prior information on the self-persistence probability.

Weak-limit sampler
The weak-limit sampler for the sticky HDP-HMM constructs a finite approximation to the HDP prior based on the fact that converges to the HDP prior when L goes to infinity. Using this approximation, one can jointly sample latent variables {z t } T t=1 with the HMM forward-backward procedure [30], which accelerates the mixing rate of the Gibbs sampler.
The main difference between our weak-limit Gibbs sampler and the corresponding sticky HDP-HMM sampler is that we now have two dimensional latent variables {z t , w t } T t=1 to sample. For step 1, we apply the forward-backward procedure to jointly sample the two dimensional latent variables {z t , w t } T t=1 .
Step 2 is the same as in Algorithm 1. For step 3, we sample β andπ based on Dirichlet-multinomial conjugacy, given auxiliary variables {m jk } L j,k=1 , the empirical transition matrix {n jk } L j,k=1 , and the approximate prior in equation 7. For step 4, we place a conjugate prior on θ j and use conjugacy to sample from the posterior.
Step 5 is the same as in By jointly sampling the full latent sequence {z t , w t } T t=1 , the weak-limit sampler greatly improves the mixing rate. The correlated observations in ARHMM and SLDS further slows the mixing rate of the direct assignment sampler, so jointly sampling is especially important for those models with dynamics [9].

Empirical results
In this section, we apply the disentangled sticky HDP-HMM to both simulated and real data. We compared the performance of our disentangled sticky HDP-HMM with the two baseline models: sticky HDP-HMM and HDP-HMM. We evaluated the model performance according to two metrics: normalized Hamming distance (defined as the element-wise difference between the inferred states and the true underlying states) on training data, and the predictive negative loglikelihood on held-out test data. To compute the Hamming distance, we used the Munkres algorithm [21] to map the indices of the estimated state sequence to the set of indices that maximize the overlap with the true sequence. To compute the predictive negative log-likelihood, we used the set of parameters inferred every 10th Gibbs iteration after it converges, and ran the forward algorithm of [30].
For all the experiments, we adopted a full Bayesian approach, and used the following hyperpriors. For all three models, we placed a Gamma(1, 0.01) prior on the concentration parameters α (and for the sticky HDP-HMM on α+κ) to cover a wide range of α values. We placed a Gamma(2, 1) prior on γ to avoid extreme γ samples, because small γ will cause numerical instability when sampling β, while big γ will generate too many states. For our model and sticky HDP-HMM, we placed non-informative priors on the self-persistence parameters. We placed a Unif([0, 1]) prior on the self-persistence proportion parameter φ = ρ1 ρ1+ρ2 . For our model, we placed a Unif([0, 2]) on the self-persistence scale parameter η = (ρ 1 + ρ 2 ) −1/3 , and cut [0, 1] × [0, 2] (the support of φ, η) into 100 × 100 grids (for simulated data) or 30 × 30 (for real data) to numerically compute the posterior for φ, η.
We used the direct assignment sampler to fit the simulated data, because the simulated data here has a relatively small sample size, and the direct assignment sampler generates samples from the true posterior of the model. We used the weak-limit sampler to fit real data (with a larger sample size) in parallel. We have also applied the direct assignment sampler to a short version of the hippocampal data; the results obtained are qualitatively similar.

Simulated data
In the simulation studies, we focus on two settings that serve to clearly illustrate the differences between the sticky versus disentangled sticky models. In both settings, we experimented with both multinomial and Gaussian emissions. See the Gaussian emission results in the supplementary material.
Different self-persistence, same transition. We simulated data using a transition matrix (Figure 2(a)) with different κ j and the sameπ j across states, which corresponds to a large transition concentration parameter α (big feature 1), and small self-persistence parameters ρ 1 , ρ 2 (small feature 3). We assigned the multinomial observation to be equal to the latent state with probability 0.9 and to other states with equal chance with probability 0.1. Note that we plot κ, α for sticky HDP-HMM in the histogram of ρ1, ρ2, since κ, α would be equal to ρ1, ρ2 respectively if we treat sticky HDP-HMM as a special case of our model. Our model learns a big α and small ρ1, ρ2, which is consistent with the data.
Results of three models fit on multinomial emission are shown in Figure 2. As shown in Figure 2(e), our model learns a big α and small ρ 1 , ρ 2 , which is consistent with the data, while sticky HDP-HMM model entangles these two parameters and learns something in the middle. Better recovery of the hyperparameters leads to better fits on the data: our model outperforms the two baseline models in terms of negative log-likelihood and Hamming distance (Figure 2(b)(c)). The advantages of our model are even more clear under the Gaussian emission (see figure in the supplementary material). The three models generally learn similar γ, which means that they infer similar numbers of states.
For each emission model, we compared these three models on 10 datasets generated from the same model using different random seeds. The conclusions are consistent across the 10 different datasets.
Same self-persistence, different transition. We simulated data using a transition matrix (Figure 3(a)) with the same κ j and differentπ j across states, which corresponds to small α (small feature 1), and large ρ 1 , ρ 2 (large feature 3). We assigned the multinomial observation to be equal to the latent state with probability 0.8 and to other states with equal chance with probability 0.2.
Again, consistently across 10 replications, our model learns hyperparameters consistent with the data (Figure 3(e)) and outperforms the two baseline models (Figure 3(b)(c)), though the advantages are not as big as in the previous scenario. This is likely because in the previous example, our model can learn the hyperparameters based on the similarity among rows of the transition matrix, while in this scenario, it can learn hyperparameters mostly from the similarity among diagonal elements of the transition matrix, which contain much less information.

Inferring rat hippocampal population codes
Next we applied our model to a public hippocampus dataset [27,28] 3 . In the experiment, a rat freely explored in an open square environment (∼ 50cm×40cm), while neural activity in the hippocampal CA1 area was recorded using a multielectrode array. See Figure 4(a) for an example trace illustrating the position of the rat over the course of the experiment.
We selected the 100 most active putative pyramidal neurons and binned the ensemble spike activities with a frame rate of 10Hz. The dataset consists of ∼ 36k frames. We used the first 8k frames, cut it into blocks of 500 frames, and randomly took 8 blocks as training data (4k frames), and the remaining 8 blocks as test data. The spike count at time t for cell c is modeled as Poisson with rate λ zt,c , i.e. y t,c ∼ Poisson(λ zt,c ).  Fig. 3. Results for the simulated data with same self-persistence (ρ1, ρ2 small), different transition (α large) and multinomial emission. We ran 3 chains, each with 30000 iterations, with 20000 iterations as burn-in for this data. Conventions as in Figure 2. Our model learns a small α and big ρ1, ρ2, which is consistent with the data.
We set L = 200 and ran 7 MCMC chains, each with 15000 iterations, with 11000 iterations as burn-in (Figure 4(d)). Our model achieves the smallest negative log-likelihood on test data (Figure 4(c)). A two-sample t-test for the mean difference between negative log-likelihood of DS-HDP-HMM and S-HDP-HMM for 7 chains is significant (p-value < 0.05). The inferred states are correlated with the spatial locations of the rat (Figure 4(b)). Compared to sticky HDP-HMM, our model infers a bigger α and smaller ρ 1 , ρ 2 (Figure 4(e)). Bigger α implies smaller variability of the switching transition, consistent with the rat quickly going to other locations since it has fast running speed. Smaller ρ 1 , ρ 2 imply bigger variability of the self-persistence probability, consistent with the rat spending different durations at different locations. Note that the disentangled sticky HDP-HMM and sticky HDP-HMM perform similarly on the rat hippocampal data; bigger differences are seen in the mouse behavior video data in the next section.

Segmenting mouse behavior video
We also applied our model to a public mouse behavior dataset [4,24]. In this experiment, a head-fixed mouse performed a visual decision task while neural  Figure 2, though with more traces and using the weak-limit sampler in (d). Our model infers a bigger α than ρ1, ρ2 (big feature 1, small feature 3). activity across dorsal cortex was optically recorded using widefield calcium imaging. We only used the behavior video data (128x128 pixels each grayscale video frame), which was recorded using two cameras (one side view and one bottom view). The behavior video is high dimensional, so we directly adopted the dimension reduction result in [1], which output 9-dimensional continuous variables estimated using a convolutional autoencoder. The dataset consists of 1126 trials across two sessions, with 189 frames each trial (30 Hz frame rate). We randomly chose 100 trials as training data (∼ 20k frames), and 30 trials as test data (∼ 6k frames).
We assumed the observation follows an ARHMM, i.e. y t ∼ N (A zt y t−1 , Σ zt ). As in [9], we standardized the observations, and assumed {A j , Σ j } follows a conjugate matrix-normal inverse-Wishart (MNIW) prior. The details of this prior are in the supplementary material. We initialized the states {z t } T t=1 of the three Bayesian nonparametric models using the state assignments result of a 32-states parametric ARHMM in [1]. We set L = 40 and ran 7 MCMC chains, each with 10000 iterations, with 9000 iterations as burn-in (Figure 5(c)). From Figure 5(a), we can see that our model has less rapid switches among states than the parametric ARHMM even if we have more states. The diagonal elements are strongly variable across states (small feature 3) (Figure 5(d)). For this dataset, α, ρ 1 , and ρ 2 are all small, but they still have different values. Our model infers a bigger α than ρ 1 , ρ 2 (big feature 1, small feature 3), as shown in Figure 5(e). Again, our model achieves much smaller negative predictive log-likelihood on test data than the other two Bayesian nonparametric models ( Figure 5(b)), due to disentanglement of hyperparameters, improving the expressiveness of the HDP prior.

Discussion and conclusion
In this paper, we propose an extension of the sticky HDP-HMM, to decouple the strength of the self-persistence prior (or the similarity of the diagonal elements of the transition matrix) and transition prior (or the similarity of the rows of the transition matrix). We develop two novel Gibbs samplers for performing efficient inference. We also show in simulated and real data that our extension outperforms existing approaches.
The work [16] proposed a Bayesian nonparametric prior to Hidden Semi-Markov Model [22] which also extends the sticky HDP-HMM. The idea is to break the Markovian assumption and instead model the distribution of duration of self-persistence for each state explicitly. The flexibility of the choice of duration distribution provides a rich set of possible priors, though the choice of duration distribution can be challenging for practitioners. Our formulation can be thought as a special case of their model, but with a few important advantages. First, we clarify that the limitations of the sticky HDP-HMM might not necessarily come from the Markovian assumption but could instead be due to the entanglement of the prior. Second, our extension adds a single hyperparameter, making it easy for practitioners to use. Third, our formulation still maintains the Markovian structure, which is useful in many applications. For example, there are applications of HDP-HMM prior to Markov decision processes [5,7,6], in which the Markovian structure is critical. Also, by harnessing the Markovian structure, our model enjoys scalable O(T K 2 ) message passing algorithms, while the message passing algorithm of HSMM has complexity of O(T 2 K + T K 2 ) in general.
In the future we hope to explore alternative computational approaches (specifically, stochastic variational [15,36] and amortized inference [25] methods), which may further improve the scalability of the model introduced here. Another important direction is to adapt the mixture of finite mixtures (MFM) model [20] to the HMM setting, as a replacement for the HDP prior; as in the mixture modeling setting discussed in [20], we expect that the MFM prior may lead to better estimates of the number of latent states.

Supplementary material
A Notation Table 1 summarizes the notation used in this paper.
General notation x1:T the sequence {x1, · · · , xT } x \t the sequence {x1, · · · , xt−1, xt+1, · · · , xT } x·j i xij xi· j xij γ concentration parameter in DP prior H prior distribution of parameters DP(γH) Dirichlet process distribution with concentration parameter γ and base measure H GEM(γ) stick-breaking distribution with parameter γ β global transition distribution α concentration parameter in transition matrix prior θj parameter of distribution of observations associated with state j πj transition distribution for state j zt latent state at time t yt observation at time t Sticky HDP-HMM κ scalar stickiness parameter Disentangled sticky HDP-HMM κj self-persistence probability for state j ρ1, ρ2 parameters in prior distribution of self-persistence probabilitȳ πj transition distribution if not self-persistent for state j wt binary variable indicating next step is self-persistent or switching n jk number of transitions from state j to k with wt = 0 in {zt} T t=1 n −t jk number of transitions from state j to k with wt = 0 in {zt} T t=1 , not counting zt−1 to zt or zt to zt+1 m jk number of tables in restaurant j serving dish k Table 1. Notation

C Derivation of Gibbs samplers
C.1 Derivation of direct assignment sampler Step 1: sequentially sample z t , w t , w t+1 The joint posterior of z t , w t , and w t+1 has the following form where y \t are all the observations except for y t , w \{t,t+1} are all the w t except for w t , w t+1 . The predictive observation likelihood p(y t |y \t , z t = k, z \t ) can be easily computed if we use a conjugate prior on the parameter in observation likelihood. For all the emissions in this paper, we used conjugate priors. See the normal-inverse-Wishart conjugate model in [23], Dirichlet-multinomial conjugate model in [35], Poisson-gamma conjugate model in [19], and matrix-normal inverse-Wishart conjugate model in [9]. p(z t = k, w t , w t+1 |z \t , w \{t,t+1} , α, β, {κ j } K+1 j=1 ) ∝ p(z t = k, w t , w t+1 , z t+1 |z \{t,t+1} , w \{t,t+1} , α, β, Let z t−1 = j, z t+1 = l, then equation 8 has the following cases Similar as [10], we have the following equations with a close-form integration. , πj p(z t = k|w t , π j )p(π j |{τ |z τ −1 = j, w τ = 0, τ = t}, α, β)dπ j · π k p(z t+1 = l|w t+1 , π k )p(π k |{τ |z τ −1 = k, w τ = 0, τ = t + 1}, α, β)dπ k where n zt−1zt is the number of transitions from state z t−1 to state z t with w t = 0 in z 1:T . n −t zt−1zt is the number of transitions from state z t−1 to z t with w t = 0, without counting the transitions z t−1 to z t or z t to z t+1 . n j· is the sum of number of transitions from j to other states.
If z t = K + 1, i.e. a new state appears, we will increment K, sample selfpersistence probability κ K+1 for a new state from prior, and update β in the following stick-breaking way. Sample b ∼ beta(1, γ), and assign β K = bβ knew , β knew = Step 2: sample {κ j } K+1 j=1 The posterior of κ j is derived as follows using the beta-binomial conjugate property where κ K+1 serves for the self-persistence probability of a new state.
Step 3: sample β To sample the global transition distribution β, we introduce auxiliary variables m jk is the auxiliary variable according to the Chinese restaurant franchise (CRF) formulation of HDP prior [34]. (Note that it can be think as number of tables in restaurant j which serve dish k.) We first update m jk , then sample β.
Step 4: sample hyperparameters α, γ, ρ 1 , ρ 2 Sampling α, γ are the same as in [34,8]. Basically, if we apply gamma prior on α and γ, then by introducing some auxiliary variables, we can have the posterior of α and γ to be gammaconjugate.

E Simulation results for Gaussian emission
We simulated data using a transition matrix (Figure 6(a)) with different κ j and the sameπ j across states, which corresponds to a large transition concentration parameter α (big feature 1), and small self-persistence parameters ρ 1 , ρ 2 (small feature 3). To generate Gaussian observations, we assumed that y t ∼ N (θ zt , 0.5 2 ), θ j iid ∼ N (3.5, 6 2 ), j = 1, · · · , 8. We placed a Gaussian prior on the Gaussian mean parameters, with mean and variance equal to the empirical mean and variance, and fixed the observation noise as N (0, 0.5 2 ). We ran 3 MCMC chains, each with 30000 iterations, with 26000 iterations as burn-in. See Figure 6 for the results.