Bayesian Nonparametric Dynamical Clustering of Time Series

Adrián Pérez-Herrero; Paulo Félix; Jesús Presedo; Carl Henrik Ek

stat.MLOct 8, 2025

Bayesian Nonparametric Dynamical Clustering of Time Series

Authors: Adrián Pérez-Herrero, Paulo Félix, Jesús Presedo, Carl Henrik Ek

Organizations: CiTIUS (Centro Singular de Investigación en Tecnoloxías Intelixentes), Universidade de Santiago de Compostela, Santiago de Compostela, 15782 Spain. · Department of Computer Science and Technology, University of Cambridge, 15 JJ Thomson Avenue, Cambridge, CB3 0FD, United Kingdom.

Abstract

We present a method that models the evolution of an unbounded number of time series clusters by switching among an unknown number of regimes with linear dynamics. We develop a Bayesian non-parametric approach using a hierarchical Dirichlet process as a prior on the parameters of a Switching Linear Dynamical System and a Gaussian process prior to model the statistical variations in amplitude and temporal alignment within each cluster. By modeling the evolution of time series patterns, the method avoids unnecessary proliferation of clusters in a principled manner. We perform inference by formulating a variational lower bound for off-line and on-line scenarios, enabling efficient learning through optimization. We illustrate the versatility and effectiveness of the approach through several case studies of electrocardiogram analysis using publicly available databases.

Figures & tables

Figure 1. Several heartbeat examples from record e1302 of European ST-T Database ( Taddei et al., 1992 ) . All of them were manually annotated as Normal except the fifth one, which was annotated as Ventricular ( Goldberger et al., 2000 ) . ECG examples with one distinct beat Line plot showing several ECG heartbeat segments from the same record. All beats except the fifth one have similar QRS-like morphology but with a considerable variation in amplitude, with a narrow positive peak followed by a negative deflection, while the fifth segment has a visibly different shape corresponding to a ventricular annotation.

Figure 2. Graphical model. Several linear dynamics are generated to represent the different clusters that evolve across time. Conditioned on a specific setting of the switch variable

s_{n}=m

, an observation

\mathbf{y}_{n}

is generated by LDS

m

through a linear projection given by equations ( 20 ) and ( 21 ) from the last observation generated by this same LDS. HDP-GPC graphical model Diagram of the proposed HDP-GPC probabilistic model. A hierarchical Dirichlet process defines global and transition probabilities for the switching states; each state selects cluster-specific Gaussian-process and linear-dynamical components, which generate an observed time-series segment after latent evolution and temporal warping.

Factor	Latent variables	Variational family	Role in inference
$q(\beta)$	Stick-breaking variables	$\displaystyle\prod_{k=1}^{K}\operatorname{Beta}\!\left(v_{k}\|\lambda_{k}\eta_{k},(1-\lambda_{k})\eta_{k}\right)$	Controls the global HDP weights $\boldsymbol{\beta}$ and penalises unused clusters.
$q(\pi)$	Transition distributions	$\displaystyle\prod_{j=0}^{K}\operatorname{Dir}\left(\pi^{(j)}\|\kappa_{j1},\ldots,\kappa_{j,K+1}\right)$	Updates Markov transition probabilities from expected transition counts.
$q(\mathcal{S})$	Switching states	$\displaystyle\prod_{n=1}^{N}\prod_{j=1}^{K}\prod_{k=1}^{K}\xi_{njk}^{s_{n-1,j}s_{nk}}$	Provides cluster responsibilities and expected transition counts.
$q(\mathcal{G})$	Warping functions	$\displaystyle\prod_{n=1}^{N}\prod_{k=1}^{K}\delta_{\widehat{\mathbf{g}}^{(k)}_{n}}\!\left(\mathbf{g}^{(k)}_{n}\right)$	Approximates temporal alignment by a MAP monotone warp.
$q(\phi)$	LDS parameters	$\displaystyle\prod_{k=1}^{K}q\left(\mathbf{A}^{(k)},\boldsymbol{\Sigma}^{(k)}_{\omega}\right)q\left(\mathbf{C}^{(k)},\boldsymbol{\Sigma}^{(k)}_{\epsilon}\right)$	Learns cluster-specific evolution and emission parameters.
$q(\mathbf{A}^{(k)},\boldsymbol{\Sigma}^{(k)}_{\omega})$	Latent transition parameters	$\mathcal{MNIW}\left(\widetilde{\mathbf{M}}^{(k)}_{A},\widetilde{\mathbf{W}}^{(k)}_{A},\widetilde{\mathbf{V}}^{(k)}_{A},\widetilde{\nu}^{(k)}_{A},\widetilde{\mathbf{S}}^{(k)}_{\omega}\right)$	Conjugate update for the LDS transition model.

Table 1. Variational factors used to approximate the posterior distribution of HDP-GPC.

Panel of Figure 3 — Figure 3. Model enhancing plasticity. A total of 50 heartbeats are represented, manually annotated as Normal, from record 114 of the MIT-BIH Arrhythmia Database [11:30.000,12:20.000]. The resulting GP $f(t)$ (left) and $x(t)$ (right), after evolving along with these sequences, were sampled over a dense time axis. The mean of the posterior is shown as a black line. A 95% confidence region for the posterior is shown in cyan blue. Kernel hyperparameters for $k_{\theta}(t,t^{\prime})$ are $\sigma^{2}_{f}=16.0^{2}$ , $\ell=2.5$ , $\sigma^{2}_{n}=5.0^{2}$ . LDS priors: $S_{\omega}=10.0^{2}I,S_{\epsilon}=5.0^{2}I$ , where $I$ is the identity matrix. For clarity purposes, no warping function has been used. This model embraces the totality of the heartbeats in a single cluster. Plastic GP model for ECG beats Two-panel ECG plot showing fifty normal heartbeat segments together with posterior Gaussian-process representatives under a plastic configuration. The posterior mean adapts to the evolving morphology, and the uncertainty band of the emission step remains small relative to the variance of the observations indicating a plastic behaviour of the model adapting to new evidence.

Panel of Figure 4 — Figure 4. Model enhancing stability. The same 50 heartbeats as in Figure 3 are represented, from record 114 of the MIT-BIH Arrhythmia Database. The resulting GP $f(t)$ (left) and $x(t)$ (right), after evolving along with these sequences, were sampled over a dense time axis. Kernel hyperparameters for $k_{\theta}(t,t^{\prime})$ are $\sigma^{2}_{f}=16.0^{2}$ , $\ell=2.7$ , $\sigma^{2}_{n}=5.0^{2}$ . LDS priors: $S_{\omega}=4.0^{2}I,S_{\epsilon}=15.0^{2}I$ . For clarity purposes, no warping function has been used. Again, the model embraces the totality of the heartbeats in a single cluster. Stable GP model for ECG beats Two-panel ECG plot showing the same fifty normal heartbeat segments under a stable configuration. The posterior mean represents an average beat morphology, and the uncertainty bands become wider enough to describe variability around this representative.

Figure 5. The GP model for each one of the beats depicted in Figure 1 is shown. Our method correctly includes all the beats in the same evolving cluster, except for the fifth one, which gives rise to a new cluster. Evolving ECG cluster models Grid of Gaussian-process models associated with the ECG beats shown earlier. Beats previously identified as normal are represented by the same evolving morphology even considering the amplitude variations, while the disruptive ventricular beat is assigned to a separate cluster.

Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$	Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$
100	4	2	3	2	201	15	10	6	10
101	4	2	3	2	202	9	7	8	6
102	10	10	9	8	203	33	7	12	38
103	10	4	5	2	205	14	10	7	9
104	16	15	10	27	207	61	9	15	36
105	10	4	10	50	208	28	17	11	30

Table 2. Number of clusters per record in the MIT-BIH Arrhythmia Database.

Figure 6. Clusters identified from record 210 of the MIT-BIH Database. The black line denotes the posterior mean, the pink region represents the 95% confidence interval for the latent process $f(\cdot)$ , and the red region corresponds to the emission process $x(\cdot)$ . Cluster 1 captures nearly all normal heartbeats despite gradual morphological evolution. Cluster 7 represents an aberrated atrial premature morphology, while the remaining clusters correspond to multiform premature ventricular contractions. ECG clusters from MIT-BIH record 210 Grid of ECG morphology clusters inferred from MIT-BIH record 210. Each panel overlays multiple aligned beat traces with a posterior mean curve and uncertainty bands. One large cluster contains most normal beats, while smaller clusters contain ventricular or aberrated morphologies.

Panel of Figure 7 — Figure 7. Respiration signal (in blue) and prediction (in red) for record f1y01, [6:30.000,8:30.000] of the Fantasia Database is shown. Training set involves a total of 500 beats from [0:00.000,6:30.000]. Respiration reconstruction for Fantasia record f1y01 Time-series plot comparing the recorded respiration signal with the respiration signal reconstructed from ECG-derived warping functions for Fantasia record f1y01. The two curves show similar oscillatory behaviour over the displayed interval. Respiration reconstruction for Fantasia record f1o03 Time-series plot comparing the recorded respiration signal with the respiration signal reconstructed from ECG-derived warping functions for Fantasia record f1o03. The reconstructed curve follows the dominant respiratory rhythm, although some local amplitude and phase differences remain. Respiration reconstruction during apnea Time-series plot comparing recorded and reconstructed respiration for Apnea-ECG record a03. The figure includes an annotated apnea interval, during which the respiratory pattern changes; the reconstruction mainly captures the respiratory periodicity with a visible change on the apnea interval.

Figure 10. Directional wave-spectra regimes inferred by the multi-output HDP-GPC model at NDBC station 41013. For each identified regime, the number of wave spectra is included in parentheses. The polar plot shows the mean directional distribution of spectral energy measured in m 2 /Hz, the right column shows the direction-integrated frequency spectra of the cluster and the GP representative. Directional wave-spectrum clusters Summary of inferred directional wave-spectrum regimes. Each cluster is represented by a polar plot of mean directional energy and a frequency-spectrum plot with the assigned spectra and Gaussian-process representative, showing the dominant direction and spectral morphology of each sea-state cluster.

Figure 11. Overview of pedestrian trajectory clustering for PB-CAM6, the 9th of January 2023, using the static HDP-GPC model. Coloured curves show the observed pedestrian trajectories, and bold curves show the GP representatives of the inferred route families. The walking direction is illustrated by an empty point indicating the initial position and a filled one indicating the ending position. The top-right diagram shows a map of the tracked zone from PB-CAM6, with green arrows indicating exits and entrances, the red box indicating the position of a shop, solid lines indicating obstacles and dashed lines indicating the limits of the camera-tracked zone. Pedestrian trajectory clustering overview Overview of pedestrian trajectories recorded by camera PB-CAM6. Many coloured paths are overlaid with thicker Gaussian-process route representatives. A small layout diagram marks entrances, exits, obstacles, a shop area, and the camera-tracked region.

Appendix figures & tables8 assets

Supplementary material from the paper’s appendix.

Appendix

Method	Clusters	Acc./Purity	ARI	AMI	NMI	Homo.	Comp.	V-meas.
Lagerholm et al. ( Lagerholm et al., 2000 )	1200	0.985	–	–	–	–	–	–
Castro et al. ( Castro et al., 2015 )	629	0.972	–	–	–	–	–	–
DTW+k-medoids ( Sakoe and Chiba, 1978b ; Kaufman and Rousseeuw, 1990 )	119	$0.943\pm 0.063$	$0.224\pm 0.340$	$0.221\pm 0.306$	$0.222\pm 0.305$	$0.288\pm 0.322$	$0.232\pm 0.330$	$0.222\pm 0.305$
DTW+hierarchical-average ( Sakoe and Chiba, 1978b ; Murtagh and Contreras, 2012 )	107	$0.934\pm 0.087$	$0.367\pm 0.399$	$0.340\pm 0.372$	$0.341\pm 0.372$	$0.355\pm 0.373$	$0.506\pm 0.413$	$0.341\pm 0.372$
CAE+DBSCAN ( Masci et al., 2011 ; Ester et al., 1996 )	168	$0.946\pm 0.064$	$0.373\pm 0.308$	$0.323\pm 0.270$	$0.325\pm 0.269$	$0.490\pm 0.302$	$0.287\pm 0.243$	$0.325\pm 0.269$
CAE+DPGMM ( Masci et al., 2011 ; Antoniak, 1974 )	906	$0.979\pm 0.041$	$0.271\pm 0.326$	$0.298\pm 0.281$	$0.301\pm 0.281$	$0.754\pm 0.278$	$0.225\pm 0.237$	$0.301\pm 0.281$

Appendix

Table 3. Global clustering results on the MIT-BIH Arrhythmia Database.

Record	Beats	Clusters	Purity	NMI	ARI	Homo.	Comp.	M-F1.
100	2272	2	0.985	0.092	0.056	0.048	1.000	0.664
101	1865	2	0.998	0.584	0.571	0.413	1.000	0.666
102	2187	10	0.991	0.452	0.328	0.872	0.305	0.932
103	2084	4	0.999	0.000	-0.002	0.002	0.000	0.500
104	2229	15	0.918	0.496	0.335	0.785	0.363	0.497
105	2572	4	0.986	0.496	0.661	0.619	0.413	0.566

Appendix

Table 4. Per-record clustering metrics for the multi-output offline HDP-GPC method.

	N	S	V	F	Q
N	90029	1101	152	67	199
S	406	1672	28	0	38
V	87	6	6934	74	46
F	68	1	46	662	0
Q	35	0	73	0	8232

Appendix

Table 5. Global confusion matrix after AAMI label aggregation.

Figure 12. Distribution of execution times across MIT-BIH records for the five inference configurations. Violin shapes summarise empirical runtime density, boxplots show the interquartile range and median, triangles indicate the mean, and points correspond to individual records. Runtime distributions by inference configuration Violin and box plots comparing execution times across five inference configurations. The online single-output configuration without warping has the lowest runtimes, while the offline multi-output configuration with warping has the largest runtimes and widest variability.

Figure 13. Empirical runtime scaling with the inferred number of clusters on the MIT-BIH Arrhythmia Database. Points correspond to individual records, solid lines show least-squares fits within each configuration, and shaded regions indicate bootstrap confidence envelopes. Runtime versus number of clusters Scatter plot of runtime against the inferred number of clusters for five inference configurations. Fitted trend lines increase with the number of clusters, and the warped and multi-output configurations occupy higher runtime ranges.

Figure 14. Empirical runtime scaling with the number of beats on the MIT-BIH Arrhythmia Database. Points correspond to individual records, solid lines show least-squares fits within each configuration, and shaded regions indicate bootstrap confidence envelopes. Runtime versus number of beats Scatter plot of runtime against the number of ECG beats for five inference configurations. Fitted trend lines show increasing computational cost for longer records, with offline warped configurations requiring the most time.

Figure 15. Memory usage by inference configuration. The left panel reports peak CPU resident set size, and the right panel reports the increase in resident memory during execution. Bars show means across records and error bars indicate one standard deviation. Memory usage by inference configuration Pair of bar charts summarising memory usage across inference configurations. The left chart reports peak resident memory, and the right chart reports the increase in resident memory during execution. Offline, multi-output, and warped configurations show higher memory usage.

Figure 16. Complete pedestrian trajectory clustering obtained with the static HDP-GPC model. The top panel shows all trajectories and the inferred GP representatives. The remaining panels show the individual route families separately. The black curve denotes the GP representative of each cluster, and coloured curves denote the assigned trajectories. Direction is represented by fading each trajectory towards its endpoint; for the representatives, the open marker indicates the initial point and the filled marker indicates the final point. Detailed pedestrian trajectory clusters Large grid of pedestrian trajectory clusters. The top panel shows all trajectories and their representatives together, while the remaining panels show individual route families. In each cluster panel, coloured curves represent assigned pedestrian paths and a black curve represents the Gaussian-process route representative.

Time Series Clustering

Time series clustering groups sequences by shared temporal patterns, regimes, trajectories, or behaviors without requiring labeled examples.

See all 27 papers

Poisson-Gamma Dynamical Systems with Time-varying Transition Dynamics
Sep 1, 2026Jiahao Wang, Yijun Wang, Nan Fang +1Dynamical Systems Probabilistic Time Series Forecasting
Department of Electrical and Computer Engineering, University of Arizona, US · School of Computing and Information Technology, Great Bay University, China · Tsinghua University, China +1
Model-free Bayesian Clustering and Multimodal Posterior Summarisation
Oct 9, 2026Qing Li, Sara WadeClustering Unsupervised Clustering
School of Mathematics University of Edinburgh
Amortized Neural Clustering of Time Series based on Statistical Features
May 13, 2026Ángel López-Oriona, Ying SunCluster Number Estimation Amortized Inference
CEMSE Division, Statistics Program, King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia. · Statistics Program, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia

Sep 1, 2026J/K move · Enter open · S save

Poisson-Gamma Dynamical Systems with Time-varying Transition Dynamics

Jiahao Wang, Yijun Wang, Nan Fang, Sikun Yang

Department of Electrical and Computer Engineering, University of Arizona, US · School of Computing and Information Technology, Great Bay University, China · Tsinghua University, China+1

Dynamical Systems Probabilistic Time Series Forecasting

Open paper AskPDF arXiv

Bayesian methodologies for handling count-valued time series have gained prominence due to their ability to infer interpretable latent structures and to estimate uncertainties. Among these Bayesian models, Poisson-Gamma Dynamical Systems (PGDSs) are proven to be effective in capturing the evolving dynamics underlying observed count sequences. However, the state-of-the-art PGDS still falls short in capturing the transition dynamics that are commonly observed in real-world count time series. To mitigate this limitation, a PGDS with time-varying transition kernel (TV-PGDS), is proposed to allow the underlying transition matrices to evolve over time. Three specifically-designed Dirichlet Markov chains (Dir-Dir, Dir-Gam-Dir, PR-Gam-Dir) are constructed to accommodate heterogeneous structural mutations within these dependencies. Leveraging Dirichlet-Multinomial-Beta data augmentation techniques, a fully-conjugate and efficient Gibbs sampler is developed to perform posterior simulation. Experiments show that, in comparison with related models, the proposed PGDS achieves improved predictive performance due to its capacity to learn time-varying dependency structure captured by the time-evolving transition matrices.

Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$	Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$
100	4	2	3	2	201	15	10	6	10
101	4	2	3	2	202	9	7	8	6
102	10	10	9	8	203	33	7	12	38
103	10	4	5	2	205	14	10	7	9
104	16	15	10	27	207	61	9	15	36
105	10	4	10	50	208	28	17	11	30

Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$	Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$
100	4	2	3	2	201	15	10	6	10
101	4	2	3	2	202	9	7	8	6
102	10	10	9	8	203	33	7	12	38
103	10	4	5	2	205	14	10	7	9
104	16	15	10	27	207	61	9	15	36
105	10	4	10	50	208	28	17	11	30

Abstract

Time Series Clustering

Explore similar work

Poisson-Gamma Dynamical Systems with Time-varying Transition Dynamics

Model-free Bayesian Clustering and Multimodal Posterior Summarisation

Amortized Neural Clustering of Time Series based on Statistical Features

Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$	Record	( Castro et al., 2015 )	$N_{MO}$	$N_{OFF}$	$N_{ON}$
100	4	2	3	2	201	15	10	6	10
101	4	2	3	2	202	9	7	8	6
102	10	10	9	8	203	33	7	12	38
103	10	4	5	2	205	14	10	7	9
104	16	15	10	27	207	61	9	15	36
105	10	4	10	50	208	28	17	11	30