Multinomial Sampling of Latent Variables for Hierarchical Change-Point Detection

Bayesian change-point detection, with latent variable models, allows to perform segmentation of high-dimensional time-series with heterogeneous statistical nature. We assume that change-points lie on a lower-dimensional manifold where we aim to infer a discrete representation via subsets of latent variables. For this particular model, full inference is computationally unfeasible and pseudo-observations based on point-estimates of latent variables are used instead. However, if their estimation is not certain enough, change-point detection gets affected. To circumvent this problem, we propose a multinomial sampling methodology that improves the detection rate and reduces the delay while keeping complexity stable and inference analytically tractable. Our experiments show results that outperform the baseline method and we also provide an example oriented to a human behavioral study.


Introduction
Change-point detection (CPD) methods aim to identify abrupt transitions in sequences of observations, for both univariate and multivariate cases.Typically, a change-point (CP) is only considered if there is a noticeable difference between the generative parameters of the data before and after the change-point event.Two classical families of approaches can be found in signal processing and machine learning for this task.First, the main focus of early literature has been on batch settings [7,14], where the entire dataset is available for analysis.Second, online CPD methods [1] avoid the previous assumption to fulfill two intertwined tasks: i) estimation of the generative parameters of the model as new observations come in and ii) segmentation of the data sequence into non-overlapping partitions based on the parameters obtained.
The identifiability of change-points (CP) is directly related to the discrepancy between the distributions governing each partition.In this context, the Bayesian framework provides a reliable solution to obtain uncertainty measures over both the parameters and the CP locations.The Bayesian online CPD algorithm (BOCPD) introduced in [1] uses this idea to derive a recursive exact inference method.However, when observations become high-dimensional and the number of parameters in the model grows exponentially, there is not enough evidence in the sequential data to obtain reliable estimates of the true generative parameters.
Latent variable models are particularly amenable to overcome the high-dimensionality issue.Under the assumption that change-points lie on a lower-dimensional manifold, one can extend the BOCPD algorithm to accept subsets of surrogate discrete latent variables.Each data point is therefore linked to a single assignment, as it is done in mixture models.The main drawback is that true latent class assignments are never observed but inferred, leading to introduce pseudoobservations [11].For this purpose, there are two main strategies: i) use the posterior probability vector as a continuous multivariate datum, i.e. as a Dirichlet distributed variable or ii) observe single point-estimates of the discrete latent Lorena Romero-Medrano and Pablo Moreno-Muñoz have equally contributed.
variable.Despite that the first idea was explored in previous works out of the CPD problem [12], it still requires expensive approximate methods due to non-tractability issues.The second idea allows reliable detection instead, particularly when posterior densities over the latent variables are certain enough.
In this paper, we consider the case of having poor inference of point-estimates over the latent variables that lead to catastrophic results on the CPD.Our contribution is to provide a novel extension for the hierarchical CP model that improves the detection rate and reduces delay even under extremely flat posterior distributions with high variance.The proposed solution considers latent variable samples as multivariate observations, that we model as multinomial distributed.This keeps the original analytic simplicity of the Bayesian CPD inference as well as the complexity cost remains significantly low.In the experiments, we prove the utility of the new inference method on synthetic data and we also provide insights to be applicable in real-world scenarios, such as change-point detection in a human behavior study.

Bayesian Change-Point Detection
Based on the work presented in [1], we assume that a sequence of observations x 1 , x 2 , … , x t may be partitioned into non-overlapping segments.We consider that each segment or partition with = {1, 2, … } has an associated generative distribution p(x| ) where the parameters are unknown and observations are assumed to be independent and identically distributed (i.i.d).The maximum number of partitions is also unknown and bounded by the total number of data points t at each time-step.Therefore, it may increase as new observations x t come in.
We are concerned with discovering the true generative distributions p(x| t ) and hence, their parameters t at each time-step.To alleviate the combinatorial problem of estimating parameters based on every partition hypothesis and time-step t, we introduce an auxiliary random variable (r.v.) r t , also called the run-length in the original version of [1].The discrete variable counts the number of time-steps since the last change-point, that is The main idea behind the run-length , r t , is that it converts the partition hypothesis problem into a Bayesian inference task as well as introduces a relatively simple CP indicator.This strategy augments the model, leading to a double inference mission: i) estimating the posterior distribution over r t given the data, p(r t |x 1∶t ) and ii) obtaining reliable values (1) of the generative parameters t conditioned to the partition hypothesis of r t .
The general method is based on the marginalization of the model parameters t for the generative distribution given the data at each time-step, and the factorization of the joint density p(r t , x 1∶t , t ) .The discrete nature of the r t counting r.v. also makes it appropri- ate for integration, being feasible to obtain the joint distribution in a recursive manner by marginalizing over r t−1 , assuming that the prior probability of r t depends only on the previous value r t−1 .Integrating in (3) as proposed in (2) leads to obtain, where we have previously defined When the computation of Ψ (r)  t is not possible, for instance, due to the underlying generative model p(x| t ) is too expres- sive or complex, other ways for approximate inference must be considered [13,15].
As we can see in equation (7), the learning process of t is conditioned to the run-length r t and hence, the partition hypothesis, carrying out a multiple-thread inference mechanism.Importantly, at each time-step t we can estimate the posterior p(r t |x 1∶t ) , obtaining a probability measure of the last CP location, given by the value of r t .Equivalently, we obtain probability measures for the location of the starting observation of the current partition.For example, having observed x 1∶5 at some time-step t = 5 , we would compute posterior estimates t |r t , x 1∶t , one per each r t value.As a consequence, given the hypothesis r t = 2 , the estimation would be analogous to see t |{x 4 , x 5 } , e.g. a CP is located two observations back, under the previous notation.This example is also depicted in the graphical scheme of Fig. 1.
1 3 However, a key inconvenient of this model appears as the size of observations x t rises, and the Bayesian method works in a potentially high-dimensional setting.In such cases, the complexity of the generative model increases accordingly to the dimensionality of x t .This leads to an extremely large number of parameters t to estimate.In fact, it makes almost impossible to perform CPD in a reliable manner as there is not sufficient statistical evidence given x 1∶t , to update our posterior distribution.In such case, CPs are typically confounded with noise drifts in the underlying parameters.

CPD and Latent Variable Models
Latent variable models are a powerful tool in unsupervised learning, with significant connections with Bayesian statistics.This family of approaches typically assumes that there exists a finite low-dimensional representation of the data that characterizes the generative properties of the observed objects.Particularly, it has become a popular solution in probabilistic modelling when the high-dimensionality problems rises.It allows to easily take decisions about the dimensionality of the latent manifold, its nature (i.e.continuous or discrete) and its conditioning to the rest of r.v.implied in the generative model of the data.
In our particular scenario for Bayesian CPD, we may assume that the sequential observations x 1∶t belong to a lower-dimensional manifold, where the true CPs lie.The generative model is then expressed as where the conditional distribution p(x t |z t ) is now assumed to be fixed and p(z t | t ) is the new likelihood distribution over the latent variable z t , that can be either continuous or discrete.With this approach in mind, we can obtain the posterior distribution p(z t |x t ) at each time-step t, allowing us to work over the latent space instead of performing parameter estimation over the initial observational space.Similar ideas were previously explored in [2,11] as extensions of the BOCPD method, where only discrete z t variables were considered.

Hierarchical CPD
Based on the previous idea, we introduce the hierarchical model presented in [11], where the latent variables at instant t, z t , are considered as categorical r.v. or classes [9], such that z t ∈ {1, 2, … , K} .Hence, they work as the assignments of each observation object x t .In the CPD scenario, it can be understood as a segmentation problem of different latent class models.Moreover, we would like to perform CP detection over the true sequence of assignments z 1∶t , but we cannot observe them.Instead, we only know the sequence of posterior distributions p(z 1∶t |x 1∶t ) that has been previously inferred via, for instance, the online expectation-maximization (EM) algorithm [5] or other continual learning strategies [10].As the chosen approach, we consider to use maximum a-posteriori (MAP) estimates of the latent variables z t as our pseudo- observations.Thus, the point-estimates are obtained from Using the strategy presented in [1], we can obtain p(r t , z ⋆ 1∶t ) as in (2).This also translates the problem of CP detection over the sequence of observations x 1∶t to perform CPD directly over the sequence of MAP estimates z ⋆ 1∶t .Importantly, to compute the posterior distribution of r t , we are also based on the marginalization over the parameters of the joint distribution.We can build the same recursiveness in (3) over the r.h.s.term by marginalizing over r t−1 inside the integral, that is where we also consider, The term p( t |r t−1 , z ⋆ 1∶t−1 ) that appears whitin Ψ (r) t is obtained by the multiple posterior updates depicted in the diagram of Fig. 1 and p(z ⋆ t | t ) is a categorical likelihood density of the new observed assignments z ⋆ 1∶t .The details of the developments and definitions can be found in [11].
Fig. 1 Illustration of the parallel inference mechanism for the estimation of t conditioned on the run-length r t given x 1∶t 1 3

The Problem of Flat Posterior Distributions
The hierarchical CP detector solves the problem of poor estimates of the likelihood parameters for high-dimensional observations.However, working with the sequence of MAP point-estimates over the latent space, may also lead to falsealarm or missing detection problems when the inferred posterior distribution p(z 1∶t |x 1∶t ) is extremely flat, e.g. it is of high variance and there is still uncertainty over the true value of z t .We have recursively seen this behavior when dealing for instance with a considerable amount of missing data over the observed sequence, or when discrimination of observations is hard.In these particular cases, the MAP estimation may not coincide with the true latent class assignment, introducing extra noise in the CPD with undesired results.

Multinomial Sampling
As presented in Section 3, we assume that the observations x 1∶t belong to a lower-dimensional manifold where the true CPs lie.We consider that the generative model is fully expressed as in ( 8) and the latent space is of discrete nature, where we are concerned to perform the CP detection.Our goal now is to obtain a better characterization of the underlying posterior distribution, and hence z t at each time step t, when the density is not well fitted.We have presented one option based on working with the sequence of MAP pointestimates obtained from the posterior p(z 1∶t |x 1∶t ) .However, if this distribution is extremely flat, we are introducing noise in the CP detection process.
We propose to generate a new type of pseudo observations of the latent variable by drawing S i.i.d.samples of the posterior density, such that z (1)  t , z (2)  t , … , z (S) t ∼ p(z t |x t ) ∀t , rather than handling a single point-estimate z ⋆ t .This allows us to work with more information about the unknown true assignment z t of the observation at each time-step.
The new approach addresses the question of how to deal with multiple subsets of S samples instead of just one variable at each time step.A potential idea would be to introduce Monte-Carlo (MC) sampling methods, but it would lead to draw S ⋅ t samples at each time step, becoming unfeasible in the long term.Alternatively, we propose to assume that samples are multinomial distributed, which has the advantage of preserving the prior-conjugacy and is still consistent with the original recursiveness of the BOCPD algorithm [1].
We recall that a multinomial distribution with natural parameters t ∈ S K and N, measures the probability that each class k ∈ {1, ..., K} has been observed n k times over N categorical independent trials with same probabilities t .This model allows us to deal with an augmented number of pseudo-observations from z t at each time t with just the cost of introducing one more parameter in the model: N = S , that is, the total number of samples drawn from the posterior.
To treat the samples as multinomial distributed, we have to transform them.Given the sampled vector z ⋆ t = (z (1)  t , z (2)  t , … , z (S) t ) ∈ {1, ..., K} S , we can define its asso- ciated counting vector c t ∈ ℤ K + where c k t ∶= ∑ S s=1 {z (s) t = k} ∀k .Each component c k t counts the times that class k has been drawn from the S trials, so that we have ∑ K k=1 c k t = S .Thus, at each time t, we can consider the counting vector c t as an i.i.d.observation of a multinomial distribution with natural parameters t ∈ S K and S ∈ ℕ.
With the previous notation in mind, and also assuming the following prior distributions for the parameters where ∈ ℝ K + , the likelihood function expression of a new observation c t is given by Using (12) and the conjugacy property, the posterior update rule of parameters has the following closed form � = + c t , which introduces a direct recursion when a new sample is observed.
Notice from the first term of ( 14) and the definition of c t that by taking the proposed multinomial model, we are not explicitly working with distributions over the S-dimensional sampled vectors themselves anymore, but over their equivalence classes.Here, two sampled vectors are considered equivalent z ⋆ is a permutation of the vector z ⋆ S 1 .We now aim to obtain the posterior distribution p(r t |c 1∶t−1 ) by building up the recursiveness suggested by [1] for the joint distribution.This would allow us to consider a measure of uncertainty of the CP locations and 1 3 estimate the generative parameters of the multinomial distribution at each partition.Following the development of ( 6), we have where again, we have defined As presented in [1,11], we have assumed that there exists a term p(r t |r t−1 ) that only depends on the value of the run- length in the previous instant, r t−1 .This term acts like a con- ditional prior probability that modulates how likely is to detect a new CP, that is, r t = 0 would be more or less likely conditioned to the previous value of run-length hypothesis r t−1 .
We also wish to infer the parameter vector (r) t related to the current run-length r t and its associated data samples, presented in (17), through the posterior p( t |r t−1 , c (r)  1∶t−1 ) .To carry out the inference method depicted in Fig. 1 we need to find Ψ (r)  t , that in fact can be seen as the posterior predictive density of the new pseudo-observation c t con- ditioned to the run length r t−1 and the previous data in the referred partition, The predictive term p(c t |r t−1 , c (r)  1∶t−1 ) has not closed form but it is still a function of the statistics of the model.Thus, its computation is straightforward where we have defined . Additionally, using both the definition of the binomial coefficient and the properties of the Gamma function, Γ(n + 1) = n! for n ∈ ℕ , we transform the previous expression to the following one: The term S grows by S at each time-step, leading to numeri- cal instabilities in the l.h.s term of (19) for high values of t.Therefore, we have considered the following expression that is numerically more stable and it is a result of the manipulations in the terms of (19).Then, it is (18) Ψ (r) t = p(c t |r t−1 , c (r) 1∶t−1 ).
This predictive prob- ability is the one that we introduce in (15) for the estimation of r t , and hence the detection of CPs.

Change-Point Prior Distribution
The prior distribution of a change-point p(r t |r t−1 ) needs to be defined to carry out the recursiveness of equation ( 15) at each time-step t.
In particular, we have considered the prior term to be time independent, as proposed in the original BOCPD version.The conditional distribution takes the form, where H( ) is the hazard function for the geometric distribu- tion with timescale parameter , that results in a memoryless process where H( ) = 1∕ is constant, as detailed in [6].We consider as a model hyperparameter that we fix, but there are also existing works [16] where is learned in an online manner.However, this usually leads to extra complexity in the model, and in our case, falls out of the scope of this work.

Definition of Change-Points
The presented recursive method a llows us to obtain p(r t | 1∶t ) and therefore, the probability that the last CP occurred a number of r t time-steps ago.This fact is given from the normalization of the joint distribution MAP estimates-based CPs.Given the posterior density we can define the sequence of likely CP locations {r ⋆ 1∶t } through the MAP estimates of the posterior distribution of the run-length, This estimation r ⋆ t of the CP hypothesis variables is the most used in the literature and the one that we use in our experiments of section 5, since it defines the most likely CPs along the time sequence.
Cumulative probability-based CPs.We propose a new alternative strategy to characterize the variable r t from p(r t | 1∶t ) .For fixed t, we consider the probability that a CP occurred in the previous n days so, in this approach, we define the sequence {r ⋆ 1∶t } by (20) .
1 3 Note that r ⋆ t = 1 for t ≤ n , leading to need at least n days of data to use this strategy for CP location definition.Example results are shown in the experiments section, considering different number of days for the cumulative probability computation.

Computational Cost
Algorithm 1 presents all steps that must be followed to obtain the sequence of CP location estimates {r ⋆ 1∶t } .As the original method of [1], the time-complexity of the general model equation ( 6) per time-step is linear in the number of data-points so far observed.In the other hand, we see that the complexity related to the number of samples, S, needed for the Multinomial CPD performance, grows linearly with S per time-step.This follows from expression (20) and the fact that ∑ K k=1 c k t = S .The low computational cost is one of the most important contributions of this work.

Experimental Results
In this section we evaluate the performance of the proposed multinomial sampling extension for hierarchical CPD.First, we study the improvements of the method (named Multinomial CPD) over synthetic data, where we may increase or decrease the quality of inference over the latent variables to prove that detection is still reliable.In the second experiment, we evaluate the method using real-world data of a monitored user from an authorized human behavior study, analyzing how we are able to reduce the delay in the whole detection process.In the third experiment, we study the performance of the method for very large number of latent classes, that is, different values of K, the dimension of the latent variables.In the experiments, we consider that a change point is detected at time-step t = t � if there is an abrupt decrease from r ⋆ t � −1 to r ⋆ t ′ , which means that the CP occurred at instant t = t � − r ⋆ t � .We set r ⋆ t < r ⋆ t−1 − 20 as the condition for detection.This can be also adapted if more precision is required.

Synthetic Data Evaluation
In our first experiment, the Multinomial CPD model has been applied to sequences of synthetic data and the results have been summarized in Fig. 2 and Table 1.Particularly, we want to evaluate the performance of the method for several sampling sizes S, drawn at each time-step and for different levels of flatness (uncertainty) of the generative posterior distribution.
We have fixed the number of CPs on the latent sequence to five, that is, six partitions, each one occurring every 100 time steps.Moreover, we have run the algorithm for T = 600 .In the experiment, the posterior distributions p(z t |x t ) of the latent variables are also simulated.This guarantees that the posterior densities are as ill-fitted as we want in every example, avoiding the intervention of inference's stochastic conditions in the results.For each partition , we have generated a set of 100 K-dim vectors ,t from a Dirichlet distribution with parameters .At the same time, these 6 K-dim vectors have been sampled from a Uniform distribution in the interval (0, ) .This parameter defines the flatness of the synthetic posterior distribution, where a lower implies a flatter generative distribution.The hyperparameter K has been fixed to 20 classes for the whole experiment.In the proposed model, each S-vector has been sampled from a Multinomial( ,t , S ) with the vector ,t previously presented.
The prior probability of the run length r t is a function of the hyperparameter −1 , which controls the prior probability of a change: the higher is , the less probable is a change.For the Multinomial CPD method (MCPD), we have defined it as a function of the number of samples = 10 S to make both comparable the terms involved in (6) and also the results in the experiment for different number of samples.The intuition behind this choice is that, for high values of S, we want the prior probability of a change to be almost zero, so that the change-point occurrence is determined from the data.However, more accurate results may be found by tuning the parameter at each particular case.For the comparison with the Hierarchical CPD method (HCPD) we considered the same values except for the hyperparameter lambda, that has been fixed to 10 20 independently of the flatness level of the simulated distributions.
In Fig. 2 we compare the MCPD (left column) for different number of samples S = 10, 50, 100, 150, 200 with the HCPD (right column) and different levels of flatness = 3.0, 10.0, 50.0 (each row).In the upper figures we can see the distributions of the latent variables or the MAP assignments at each time step, respectively.In the bottom figures the MAP estimates of the run-length r t are jointly shown with dashed lines indicating the true change-points.We have also summarized the results of running the method five times for each pair of values (S, ) in Table 1.There, we also show the total precision rate, defined as the ratio of change points detected for each pair, and the mean and standard deviation of the delay, defined as the time points between the instant of the detection t = t � and the real instant t = t � − r ⋆ t � in which the CP occurred.For example, if a CP is detected at t = 150 and r ⋆ 150 = 30 , this means that a change occurred at t = 120 , and the delay of the detection would be 30 steps.Looking at the results included in Fig. 2, with the MCPD we detect the five change points for every value of and many of the values of S considered.In the table we confirm that the precision increase as S grows, detecting less change points when the distribution is highly flat for lower values of S and in particular for the HCPD, that is similar to the limit case in which S = 1 .For = 2.0 no change points are detected using the HCPD method.However, with the MCPD, even if the distribution is that flat we are able to find the change points by increasing the number of samples, obtaining a precision of 88% for 50 samples at = 3.0 ver- sus the 20% in the HCPD case, or even a 100% of precision already at = 4.0 when S = 100.
For higher values of , we can see both in the Fig. 2 and the Table 1 that the performance is good enough for both methods in terms of CPD precision.However, the delay of the detections is always notably lower in the proposed MCPD.In comparison to the HCPD, we can see in the table that the average of the delay in the detections is reduced by more than a half when 100 samples are considered, independently of the flatness of the distribution, with just 23.08 time steps of average delay when = 4.0 or 13.1 when = 10.0.

Computational Cost Analysis
We analyze the precision rate and average delay in the detection for several sampling sizes and their associated computational cost.In Fig. 4 we show these measurements from S = 10 to S = 200 evaluated each 10 samples using the MCPD proposed method.The first point of every line ( ≈ S = 1 ) has been obtained using the HCPD method for every row plot.We consider the same synthetic data generated for subsection 5.1 with K = 20 latent classes, = 10 S for the MCPD experiments and = 10 4 for the HCPD results.In the first row of Fig. 4, we show the computational cost evolution in seconds (red line) as the number of samples increases, jointly with an adjusted linear regression (blue line).As commented in subsection 4.3, we see that the computational cost increases linearly with the number of samples.We have fixed = 3.0 for this experiment, but notice that the computational cost is constant for any level of flatness.The method has been run 3 times and, in the second and third row of the figure, we show the total precision rate (range 0.0-1.0)and average delay (range 0-100 time-steps) in the detection for the same range of sample sizes but different levels of flatness = 2.0 , = 3.0 , = 4.0 , = 6.0 and = 10.0 .Recall that higher corresponds to lower flatness.In this experiment, the delay has been computed just for the detected CPs, leading to some cases of lower delay values for higher flatness.
We see that both the precision rate and the delay tend to the maximum (minimum for the delay) the method can get for each flatness level as the number of samples increases.The precision rate reaches this value around 50 samples for ≥ 3.0 , that is 1.0 for ≥ 4.0 .The delay presents faster stabilization, due to the computation just over the detected CPs.With respect to the associated computational cost, it increases just 0.3 seconds from the use of the HCPD method (1.1 sec) to the use of MCPD (1.4 sec) for S = 50 samples.Looking at these results and taking into account that the maximum accuracy of the method is superseded to the quality of the original data, we can conclude that a size of S = 50 samples could give enough good results in most of the cases.Anyway, choosing a higher sampling size like 100 or even 200 would increase the computational cost in just 1.0 second, ensuring to achieve the maximum precision rate and lower delay independently of the data quality level.

CPD for Large Number of Latent Classes
In this work we have discussed an approach to lead with high-dimensionality data by working over the latent representation of the observations.In this experiment we study the performance of the method in the sense of precision and delay in the detection when the number of latent classes, K, is large, for different levels of flatness, , of the posterior distribution over the latent variables.The number of CPs has been fixed to 5 on the latent sequence, that is, six partitions, each one occurring every 100 time steps.Moreover, we have run the algorithm for T = 600 and S = 100 for every pair ( , K) .The hyperparameter has been fixed to 10 S as before.The details in the generation of the data are the same as explained in subsection 5.1.For metrics computation, we compare with the ground truth, but assuming that a CP is considered as not detected if the delay is higher than 100 time-steps.The delay of a not detected CP computes as 100 for the total metric of a trial.
In Tables 2 and 3 we compare the precision and delay respectively for the Multinomial CPD and Hierarchical CPD considering different number of latent classes K = 10, 20, 40, 50, 100, 200 (columns) and different levels of flatness = 3.0, 4.0, 5.0, 10.0, 20.0 (rows).The precision is presented as the total ratio of detected points over the trials and the delay is shown as the mean and standard deviation of the delay of every CP and every trial.
In the results we see that, in general, the higher is K, the lower is the detection rate and the higher is the average delay.In terms of precision, the Multinomial CPD is able to detect more than 92% of the CPs when the flatness is higher than 4.0.Even though, the detection rate for = 3.0 is always over 84% when K is lower than 100.For these values of flatness, the HCPD precision is lower than 64% when we have more than 40 classes, as expected.
The delay of the MCPD is less than a half compared to the HCPD results of Table 3 for every pair ( , K) .In fact, the average delay in the detection when the flatness is higher enough is not larger than 18.2 for the MCPD even for 100 and 200 classes.In the case of the HCPD, the average delay is always over 80 for K = 200 as expected due to the low detection rate and the value considered in the not-detected cases to compute this metric.

Cumulative Probability Metric
In subsection 4.2 we have proposed an alternative characterization of the variable r t for CP definition from the posterior p(r t | 1∶t ) : to consider the cumulative probability that a CP happened in the previous n days as a CP location indicator.We present an example of this metric using the synthetic data generated in subsection 5.1 for = 4.0.
In Fig. 5 we show the output of the detection over the mentioned dataset for the MCPD method with S = 50 sam- ples, K = 20 latent classes and = 10 50 .In the first row we have the posterior distribution of the latent class indicator z t along time.In second, third and fourth rows we show the sequence of MAP estimates of the run length (green line) jointly with the negative logarithm of the cumulative probability,− log p(r t ≤ n | 1∶t ) , for the 5 (blue line), 10 (orange line) and 15 (brown line) previous days, respectively.The true change points are indicated with dashed lines and the cumulative probability for t ≤ n is not plotted because it is 1 by definition and therefore, not informative for our goal.
We see that the cumulative probability based-approach is consistent with the use of MAP estimates for CP definition.However, it could brings to reduce even more the delay in the detection by considering that a CP occurred if there is a sufficient growth of the cumulative probability at a particular time-step.These increases are presented in the images as peaks, that coincide or occur some time-steps before the MAP estimate peaks.Clear examples happen at t = 400 and t = 500 where the delay could be reduced to 0 for the first CP or to more than a third for the second CP with respect to the MAP strategy.distribution over the latent variable.We can see that the method finds three change points around day 100, day 230 and day 290, clearly partitioning the time in four behavioral periods between the first and last day of monitoring.These changes have not been contrasted with external information of the user yet, but the results are consistent in terms of number of detections for every value of S considered, and seem to be coherent with the overview of the distributions in the third raw of the figure.Moreover, we can see that increasing the number of samples at each time step, we can reduce the delay in the detection almost 50 days w.r.t. the hierarchical CPD method.

Comparative CPD Results
In this section, we show comparative CPD results for four different Bayesian approaches and two additional methods from the classical literature.Particularly, we use data from the magnetic response obtained during the drilling of a well.This dataset has been previously used in the context of CPD methods by [1].We remark that the true location of CPs is not provided.
The methods considered for the comparison are i) the Bayesian CPD algorithm [1], ii) Hierarchical CPD [11], iii) the infinite-dimensional method of [10] and iv) the Multinomial-based approach proposed in this work.The detection curves are shown in Fig. 6, where we observe that the better performance comes from the proposed approach using Multinomial samples.The Bayesian CPD algorithm does not include the latent variable hierarchy, which in cases with corrupted, missing or heterogeneous observations is necessary.Its performance is a bit more noisy than ours around the time-step t = 1000 .The hierar- chical CPD method is similar to ours but only uses one single MAP estimate from the underlying distribution p(z t |x t ) and in this case, it fails when the characterization of CPs requires a higher precision, e.g.around t = 1200 .This is understandable as the sampling methodology allows us to better characterize the latent variable distribution in that sections of the signal.Moreover, the infinite-dimensional approach of [10] which does not consider a fixed dimensionality for the latent space performs similar to our method.However, it does not detect small transitions in the shortlength scale of the time-series, as we do.Examples of these CPs can be found in t < 1000 and t > 3500 from Fig. 6.
Additionally, we compared our method with other non-Bayesian CPD approaches, which are included in the ruptures library for offline CPD methods. 1 In particular, we considered optimal partitioning [8] and the binary segmentation [14] method.The detection results were similar for both methods, and they only marked CPs in the t range [1000,2500].This means that the short-length transitions at the beginning and end of the signal were not considered as the other Bayesian methods do.In this context, it is worthy to mention the work of [4], where a thorough evaluation of CPD methods is performed w.r.t. a large benchmark of datasets.The final results shed light on the advantage of considering Bayesian CPD methods in certain cases, similarly as we see in our experiments.

Conclusion
In this paper, we have presented a novel methodology for improving the Bayesian CPD algorithm of [1] with latent variable models.Under the assumption that CPs lie in a lower-dimensional manifold, inference is carried out with pseudo-observations based on posterior point-estimates of the latent variables given the data.We introduce a multinomial-sampling method that improves the detection rate and reduces the delay when we treat with high-dimensional sequences of observations.The analytical tractability in the inference is maintained as well as a low computational cost.The experimental results show significant improvements in the CPD as posterior estimates become less certain.Interestingly, even under a good inference performance, the multinomial sampling method reduces the delay of detection, what in practice is a key point for its application to real-world problems.We illustrate an example on a human behavioral study, that detects changes in the circadian patterns of a user.In future work, this could be integrated with other CPD methods that consider the dimensionality of the latent variables unbounded [10].Prior to this, he held different teaching positions at Universidad de Vigo, Universidad Politécnica de Madrid, and Universidad de Alcalá, all of them in Spain.He has participated in more than 70 projects and contracts and has coauthored more that 90 journal papers and more than 130 international conference papers.His research interests include signal processing, machine learning, and information theory methods, and its application to health and sensor networks.

Fig. 2
Fig. 2 Comparison between the multinomial CPD, based on sampling from the latent class posterior, and the baseline CPD method.The resulting CPs (bottom figures) are considered as jumps over the MAP estimates (solid lines) of the run-length r t ∀t .Dashed lines indicate the true change points.Left Column: Each row represents an exam-

Fig. 3 Fig. 4
Fig. 3 Human behavior CPD with heterogeneous daily mobility metrics from a user.Three upper rows.Respectively, 310 days of distance wandered, presence at home and number of steps every 30 minutes.Fourth row.Posterior expectations over the K = 20 latent class indicator z t .Fifth row.Hierarchical CPD for several multinomialsampling cases

Fig. 5
Fig. 5 Cumulative probability measurements.First row.Posterior distributions over the K = 20 latent class indicator z t .Second, third and fourth row.Results of the MCPD method for S = 50 samples: MAP estimates (green line) of the run-length r t ∀t and negative log of the cumulative probability for 5 (blue line), 10 (orange line) and 15 (brown line) days, respectively.Dashed lines indicate the true change points

Fig. 6
Fig.6 Comparative CPD results for four different methods based on Bayesian inference.First row: Univariate data corresponds to the nuclear magnetic response obtained during the drilling of a well.Sec-

1 3 at
Universidad Carlos III de Madrid and Evidence-Based Behavior (eB2).Her research interests include mathematical modeling, machine learning, probabilistic methods, and its applications to human behavior and biomedical sciences.Pablo Moreno-Muñoz obtained his B.Sc. and M.Sc in Telecommunication Engineering from Universidad Carlos III de Madrid, Spain in 2014 and 2016, respectively.In 2015, he held a 6-month traineeship at European Space Agency for the investigation of probabilistic methods applied to distance calculations in astronomy.He is currently PhD student at the Dept. of Signal Theory and Communications also at the Universidad Carlos III de Madrid.During the last years, he has been research visitor at the University of Sheffield, UK and the Max Planck Institute for Intelligent Systems in Tübingen, Germany.His research interests include probabilistic machine learning methods for heterogeneous data, Gaussian processes, change-point detection, continual Bayesian inference, and its application to human behaviour modelling in medicine.Antonio Artés-Rodríguez was born in Alhama de Almería, Spain, in 1963.He received the Ingeniero de Telecomunicación and Doctor Ingeniero de Telecomunicación degrees, both from the Universidad Politécnica de Madrid, Madrid, Spain, in 1988 and 1992, respectively.He is a Professor at the Department of Signal Theory and Communications, Universidad Carlos III de Madrid, Madrid.