There has been a growing interest in applications where data about individuals are not accessible, instead aggregate population-level observations in the form of counts of the individuals are available [31, 18] . For various reasons including measurement fidelity, privacy preservation, cost of data collection, and scalability, data is often collected as aggregates. For example, in human ensemble flow analysis, individual trajectories may not be readily accessible due to privacy concerns, but the number of individuals in a certain geographical area can typically be counted by cell phone carriers. More examples include voter turnout based on demography from census data [15] and bird migration analysis [33] . One fundamental part in modeling such aggregate data is the training phase for estimating the model parameters. Learning the underlying individual model from aggregate observations is a challenging task since the full trajectory of each individual is not accessible. We are interested in learning hidden Markov models (HMMs) using aggregate data. The HMM is a popular graphical model used in various scenarios involving unobservable (hidden) data sequences arising in ecology, social dynamics, and emergence of an epidemic [28, 6, 8, 32] . Due to their ability to address the nonstationarity in observed data sequences, HMMs are capable of modeling a rich class of problems. In aggregate HMM settings, a large set of homogeneous individuals transit from one state to another according to the underlying HMM and at each time-step, corresponding aggregated observations are recorded. For example, in epidemiology, one can model spread of an infectious disease such as COVID-19 over time in a geographical area using the population level aggregate data generated by an HMM. In this work, we consider the problem of estimating the parameters of a time-homogeneous hidden Markov model, i.e., transition and observation probabilities, from noisy aggregate data.
A traditional method for learning HMM is the Baum-Welch algorithm [1, 2] , which is a special case of the expectation-maximization (EM) algorithm [7] . For the given observations sampled from a model consisting of latent variables (variables that are not observable) with unknown parameters, the EM algorithm aims to find the maximum likelihood estimates of the model parameters. In its first step (E-step), the EM algorithm estimates a function of the expected values of the latent variables and subsequently in the second step (M-step), it finds the maximum likelihood parameter estimates. For the case of learning HMM parameters, inference algorithms such as belief propagation (BP) algorithm [27] is utilized in the E-step of the EM algorithm. The Baum-Welch algorithm for estimating an HMM uses the forward-backward inference algorithm, one type of BP algorithms, in the E-step to complete the data. Unfortunately, traditional HMM learning methods such as Baum-Welch algorithm [2] can not be applied to aggregate setting. Learning the individual model from such population-level observations becomes challenging since great amount of information about individuals is lost due to data aggregation and observation noise.
Recently, the learning and inference problems in aggregate settings have been formalized under the collective graphical model (CGM) framework [31] . Within the CGM framework, for learning the parameters of the individual model, several aggregate inference methods such as non-linear belief propagation (NLBP) [33] and Bethe-RDA [35] algorithms has been utilized in the E-step of the EM algorithm aiming to maximize the complete data likelihood. Both of the inference algorithms work on an explicit observation model. In addition, since NLBP does not exhibit convergence guarantee, it does not lead to stable learning methods.
The primary contribution of our work is a novel algorithm for estimating the HMM parameters with theoretical guarantees from noisy aggregate observations. We utilize a modified EM algorithm for the learning task, where the E-step of the algorithm is solved using recently proposed aggregate inference method, the Sinkhorn belief propagation (SBP) algorithm [32] . We show that our algorithm exhibits a convergence guarantee. Moreover, our algorithm naturally reduces to the standard Baum-Welch algorithm when the observed data is based on a single individual. We further extend our algorithm to learn the model parameters with continuous observation noise model. We evaluate the performance of our algorithm on a variety of scenarios including bird migration and human ensemble flow on real-world dataset.
Related Work: Estimating Markov chains from aggregate data, also referred to as macro data in earlier works, has a long history. It was first studied in [17] where the transition matrices were estimated based on maximum likelihood method. In [34, 21, 14] , the modeling of a single Markov chain was studied by maximizing the aggregate posterior. More recent learning methods from aggregate data include [18, 25] . After the CGM framework proposed by [31] , there have been a few works on learning the underlying individual model from aggregate data. The non-linear belief propagation algorithm [33] , a message passing type algorithm for approximate inference in CGMs, has been utilized in EM for the task of learning a Markov chain. Another existing aggregate inference algorithm utilized in the E-step of the EM algorithm is Bethe-RDA [35] which exhibits convergence guarantees. Finally, [4] proposed a method of moments estimator for learning a Markov chain within the CGM framework. Other works along this line include estimating spatio-temporal population flow [11] and recurrent estimation of HMM [19] from aggregate data, learning stochastic behaviour of aggregate data [20] , and estimating group behavior from ensemble observations [39] .
The rest of the paper is organized as follows. In Section 2, we briefly discuss related background including probabilistic graphical models, collective graphical models, and the Sinkhorn belief propagation algorithm. We present our main results and algorithms in Section 3 for discrete observations. The counterpart with continuous observations is developed in Section 4. This is followed by experimental results in Section 5 and a concluding remark in Section 6.
Our algorithmic framework for learning HMM is based on a modified version of the EM algorithm that utilizes SBP for inference over a graphical model. In this section, we present related background including probabilistic graphical models, their extension to aggregate settings, and the SBP inference algorithm.
A probabilistic graphical model (PGM) [36] represents the dependencies between a collection of random variables using a graph. Let the set of J random variables be X 1 , X 2 , . . . , X J , where each variable takes one of the d possible discrete states from X (|X | = d). Consider a graph G = (V, E) with the set of nodes of the graph V representing random variables and the set of edges E representing dependencies between the variables. Then, the joint probability of the distribution of random variables can be written as
where x 1 , x 2 , . . . , x J are realizations of the corresponding random variables, ψ ij are edge potentials, and Z is a normalization constant. The hidden Markov model is a special PGM.
There are two fundamental problems in PGMs: learning and inference. The learning problem in PGMs is concerned with estimating the parameters and structure of the underlying graphical model using observation data sampled from the model. The inference problem aims to infer the statistics of the node variables given the model parameters. The inference algorithm is a key component of learning algorithms such as EM algorithm [23] .
Belief Propagation: One of the most effective inference algorithms in PGMs is belief propagation (BP) [27] , which estimates the marginal distribution of each node based on the messages from the neighbors of the node. Let m i→j (x j ) be the message from variable node i to variable node j, then
where N (i) is the set of neighboring nodes of i, and N (i)\j denotes the set of neighbors of i except for j. The message in the above equation represents belief of node i on the marginal of node j. The BP algorithm iteratively updates the messages in (2) over the graph.
The BP algorithm is guaranteed to converge globally when the underlying graph is a tree [37] and one can recover the true marginals exactly upon convergence aŝ 
J } contains a value assignment to all the J variables in the model. Then, the maximum-likelihood method estimates the model parameters Ψ (set of all the edge potentials) by maximizing the log-likelihood L(Ψ; x (1) , x (2) , . . . , x (M ) ) = M m=1 log p(x (m) |Ψ). However, in practice due to economic or feasibility reasons, some variables are not observable (known as hidden or latent variables) and observations are made corresponding to only a subset of variables.
Expectation-maximization (EM) [5] is a general technique for maximum-likelihood estimation of the model parameters in presence of hidden (unobservable) variables. The EM algorithm is an iterative method that involves two steps in each iteration: E-step and M-step. In the E-step, the values associated with the hidden variables are estimated to make the data complete and then, in M-step, the parameters of the underlying model are optimized based on the complete data likelihood.
Assume that only a subset of variables Γ ⊂ V are observable in the given graphical model. (1) , o (2) , . . . , o (M ) } represents the observed data and Ψ be the unknown set of parameters to be learned. The E-step of the EM algorithm computes an auxiliary function which is the expected value of L(Ψ; X, o) given the observed data o and the current estimate of parameters Ψ old :
Subsequently, in the M-step the parameters are estimated as
The above two steps are repeated iteratively until convergence. It has been proven that the EM algorithm is guaranteed to converge at least to a local maximum. In case of estimating the parameters of a PGM, the BP algorithm is used in the E-step of the algorithm.
Collective graphical model (CGM) [31] is a framework for learning and inference from noisy aggregate data derived from a graphical model describing the behavior of individuals.
The aggregate data is generated from M independent individuals following a certain individual graphical model as in (1) . Denote the state of the m th individual at node i as random variable X 
The vectors n 1 , . . . , n J constitute the aggregate data and the aggregate edge distributions n ij represent sufficient statistics of the individual model [31] . Let n = {n i , n ij } denote the the aggregate node distributions n i together with the aggregate edge distributions n ij . In its original form, the CGM explicitly models the observation noise as a conditional distribution p(y|n) with y being the aggregate noisy observations. The learning problem in CGMs is concerned with estimating the individual model parameters of the underlying graph 1 from noisy aggregate observations y. For learning the parameters of the model, the EM algorithm is utilized that consists of two operations: Expectationstep (computes the complete data log-likelihood) and Maximization-step (maximizes the complete data loglikelihood). The E-step of an EM algorithm for learning the individual model from aggregate data requires inferring the aggregate node distributions n.
The goal of inference in CGMs is to estimate n from the aggregate noisy observations through the conditional distribution p(n|y) ∝ p(n)p(y|n), where p(n) is known as the CGM distribution [31] which is derived from the individual model (1) . For the tree structured graph, the CGM distribution is given by
where Ψ = {ψ ij (x i , x j )} is the set of parameters of the graphical model and d i is the number of neighbors of node i in the underlying graph G. The first term in the above equation is a count of the number of different ordered samples contributing to the sufficient statistics n; it does not depend on the parameters Ψ. The second term is the joint probability of the entire population. Moreover, the support of the CGM distribution p(n) is such that each entry of n is an integer and satisfies the following constraints
The exact inference of n based on p(n|y) is computationally intractable for large populations as the computational complexity increases very quickly with increase in the population size M [30] . It was first discovered in [30] that, − ln p(n|y) can be approximated by (up to a constant addition and multiplication) the CGM free energy
After relaxing the integer constraints on n i (x i ), n ij (x i , x j ) and under the assumption that the observation model p(y | n) is log-concave, the resulting problem of minimizing F CGM is a convex optimization problem. This is the approximate MAP [30] inference problem in the CGM framework. Note that the problem size of minimizing F CGM is independent of the population size M . Some of the earlier algorithms for aggregate inference problem in CGMs include non-linear belief propagation (NLBP) [33] and Bethe-RDA [35] . These algorithms are based on an explicit noise model that utilize standard BP on a modified set of model parameters. In contrast, recently proposed SBP algorithm [32] consider the observation noise within the graph and exhibits convergence guarantees. Next, we discuss the SBP algorithm in detail that we use to solve the E-step of the EM algorithm for the purpose of learning the underlying model parameters.
Sinkhorn belief propagation [32] is a recently proposed algorithm to solve the inference problem in CGMs efficiently based on Sinkhorn algorithm for multi-marginal optimal transport [24, 26, 3] . As compared to the other existing methods including NLBP and Bethe-RDA for aggregate inference under the CGM framework that explicitly model the observation noise, the observation noise model employed in SBP is incorporated within the underlying graph, ensuring that it reduces to standard BP in case of fixed delta observations. Moreover, SBP enjoys convergence guarantees and is faster than its counterparts. The generative model in the SBP algorithm is such that the noise model is incorporated withing the underlying graph. This can be explained via the graph depicted in Figure 1 , where the aggregate observations are made corresponding to variables X 4 and X 5 , whereas rest of the variables are unobservable. More generally, let G = (V, E) be the underlying CGM with joint aggregate distribution n and let Γ ⊂ {1, . . . , J} be the set of indices representing observed variables such that n i = y i for a given set of observations y i , for i ∈ Γ 2 . Here, without loss of generality, it can be assumed that the observations are recorded corresponding to a subset of leaf nodes of the underlying graph G [32] .
The goal of the aggregate inference algorithm is to find such aggregate distribution n that maximizes the posterior log p(n | y) ∝ log p(n, y) given model parameters Ψ (set of edge potentials). The SBP algorithm employs the approximation of the log-likelihood of the CGM distribution −log p(n, y | Ψ) by F Bethe (n, y | Ψ) and solves the following
where F Bethe (n) is the the Bethe free energy [38] given by
Moreover, (8b) corresponds to the aggregate observation constraints, (8c) are consistency constraints, and (8d) are the normalization constraints. All the steps of the SBP algorithm are listed in Algorithm 1. Note that in the step (ii) of the Algorithm 1, i next is next node followed by i inΓ withΓ representing a sequence of nodes such that each element in Γ appear inΓ infinitely often. The SBP algorithm is guaranteed to converge when the underlying graph is a tree [32] .
One can interpret the expressions (10) and (11) in Algorithm 1 as messages between nodes, similar to the standard BP algorithm. The messages (10) can be understood as a scaling step, which guarantees that the constraints (8b) remain satisfied.
Remark 1 It is worth noting the connection between collective and standard filtering problems. In case of a single individual case (M = 1), the aggregate distributions are concentrated at the fixed observations and subsequently, the SBP algorithm for aggregate inference coincides with the standard BP. We refer the reader to [10] for more details.
In this section, we first propose a learning algorithm called approximate EM algorithm with collective data in general graphs. We then specialize it to learning the parameters of a discrete hidden Markov model (HMM) from aggregate data.
Initialize the messages m i→j (x j ) to 1,
ii) Update all the messages on the path from i to i next
end for end while Compute the required aggregate distributions as
3.1 Aggregate Learning based on Approximate EM Algorithm The EM algorithm [22] increases the log likelihood of the observed data through maximizing expected complete data log likelihood iteratively. However, in the aggregate setting, it becomes intractable [30] to calculate expected complete data log likelihood Q(Ψ, Ψ ( −1) ) = E p(n|y,Ψ ( −1) ) [log p(n, y|Ψ)] at iteration precisely. Fortunately, as pointed out in [30] , the conditional distribution p(n|y, Ψ ( −1) ) concentrates on the minimizer n * of the Bethe free energy in (8), i.e., n * = argmin n F Bethe (n, y | Ψ ( −1) ). This can be explained by the large deviation theory. We refer the reader to [9, 32] for more justifications on this approximation.
Subsequently, the auxiliary function is approximated as
with n * = argmin n F Bethe (n, y | Ψ ( −1) ). The approximation error vanishes as the population size M increases. Based on (12), we propose approximate EM algorithm for parameter learning in CGM as listed in Algorithm 2.
The convergence and effectiveness of Algorithm 2 are characterized by Theorem 1 below. (2)) converges. Moreover, the data likelihood approximately increases after each iterations.
Proof. See Appendix A.
The approximate EM algorithm given by Algorithm 2 works on a general graph. The maximization-step may requires solving another optimization problem. Next we turn our attention to learning discrete HMMs from aggregate data, where the maximization-step has a closedform solution.
X 1 X 2 X T O 1 O 2 O T
An HMM is a Markov chain where the variables are not directly observable, but corresponding noisy variables are observed. Denote the unobserved variables as X 1 , X 2 , . . . and observed variables as O 1 , O 2 , . . .. Assume that each hidden variable takes one of the discrete values from a finite set X h and each observation variable from a finite set X o , where in general |X h | = |X o |.
A time-homogenous HMM is parameterized by the initial distribution π(X 1 ), the state transition probabilities p(X t+1 | X t ), and the observation probabilities Note that the model given by (13) is a directed graphical model but it can be equivalently converted to the undirected graphical model (1) by considering the transition and observation probabilities as edge potentials ψ ij in (1) . Denote the set of parameters to be learned
The aggregate data constitute n 1 , n 2 , . . . , n J (J = 2T ) and let the aggregate observation be y i , ∀i ∈ Γ. For collective HMM, the aggregate observation model for length of T = 3 is depicted in Figure 3 , wherein aggregate data is generated based on M number of individual trajectories 3 . 
In case of collective HMM, the SBP algorithm is nothing but the generalization of traditional forward-backward algorithm [16] to aggregate settings, known as collective forwardbackward algorithm [32] , which we present next.
Collective Forward-Backward Algorithm: The collective forward-backward algorithm (CFB) is a special case of the general SBP algorithm (Algorithm 1) when the underlying graph is characterized by an HMM. Moreover, for the case of single observation trajectory (M = 1), it reduces to standard forward-backward algorithm. The messages in collective HMM are shown in Figure 4 with α t (x t ) being the messages in the forward direction and β t (x t ) are the messages in the backward direction. Moreover, γ t (x t ) denote the messages from observation node to hidden node and ξ t (o t ) are the messages from hidden node to observation node. These messages are characterized by
with boundary conditions α 1 (x 1 ) = π(x 1 ) and β T (x T ) = 1.
The sequence of update steps are listed in Algorithm 3.
Initialize all the messages α t (x t ), β t (x t ), γ t (x t ), ξ t (o t ) while not converged do Forward pass: for t = 2, 3, . . . , T do i) Update γ t−1 (x t−1 ) ii) Update α t (x t ), ξ t (o t ) end for Backward pass: for t = T − 1, . . . , 1 do i) Update γ t+1 (x t+1 ) ii) Update β t (x t ), ξ t (o t ) end for end while Estimate required marginals as n t (x t ) ∝ α t (x t )β t (x t )γ t (x t ), (15a) n t,t+1 (x t , x t+1 ) ∝ p(x t+1 |x t )α t (x t )γ t (x t ) β t (x t+1 )γ t (x t+1 ) (15b) n t,t (x t , o t ) ∝ p(o t |x t )α t (x t )β t (x t ) ξ t (o t ) (15c)
Now we address the problem of parameter learning in an approximate way via EM algorithm with aggregate observation noise model described above. Theorem 1 is applicable to general graph with aggregate observations. The practical implementation requires efficient E-step and M-step. In this subsection, we show there exists an efficient algorithm for aggregate HMM. We propose to implement the E-step via the SBP algorithm and derive a practical method to address the maximization-step characterized via Proposition 1.
The Maximization-step updates in learning aggregate HMM are given by
.
(16c)
Proof. See Appendix B.
Remark 2 If parts of the parameters are known, then we only need to update the other parameters in the maximization step. For instance, if the emission probability are known, then only the steps (16a)-(16b) are needed.
Based on Proposition 1 along with the Approximate EM algorithm (Algorithm 2), we list the steps to learn the individual HMM parameters in Algorithm 4. Moreover, Algorithm 4 converges because of Theorem 1.
Require: Aggregate observations y 1 , y 2 , . . . , y T Initialize π(x 1 ), p(x t+1 | x t ), and p(o t | x t ) repeat Compute the hidden counts n t,t+1 (x t , x t+1 ), n t (x t ), n t,t (x t , o t ) using CFB π(x 1 ) = n 1 (x 1 )
until convergence Algorithm 4 is for a single sequence of aggregate data generated from a certain number of samples. Learning from multiple sequences was initially explored in the Baum-Welch algorithm [29] to improve learning process. Sharing with same merits, Algorithm 5 extends to an ensemble of K number of aggregate observation sequences generated from the same HMM model. Here K number of sequences are generated by aggregating trajectories of N individuals such that each sequence results from M individual trajectories, i.e., K = N/M . Proof. See Appendix C.
Next, we extend our learning framework to HMMs with continuous emission densities, more specifically, to HMMs with Gaussian emission probabilities.
Require: Ensemble of aggregate observations y k 1 , y k 2 , . . . , y k T , for k = 1, 2, . . . , K
until convergence 4 Learning aggregate HMMs with continuous observations Now we turn our attention to the problem of estimating the parameters of continuous HMMs from a set of continuous observation trajectories. A continuous HMM is similar to the discrete HMM except for the continuous emission densities, i.e., instead of taking values from a finite number of discrete symbols, the observations are allowed to take values from continuous s-dimensional observation space R s . In standard HMM case, the continuous observation model has been studied in [12, 13] . In this section, we extend our learning algorithm to continuous emission densities in aggregate observation settings. Although all the observation trajectories are recorded, the assignments to the corresponding state is not known. Now that the observation space is continuous, we represent the assignments of the sample observations as n , has been generated by hidden state x t . Recently, the inference problem in aggregate HMMs with continuous emission densities has been studied in [40] . It was shown that the required marginals can be estimated as (Corollary 2, [40] ) T , ∀m = 1, 2, . . . , M Initialize π(x 1 ), p(x t+1 | x t ), µ(x t ), and Σ(x t ) repeat
Compute the hidden marginals
until convergence that are the fixed points of the following updates
with α 1 (x 1 ) = π(x 1 ) and β T (x T ) = 1.
The inference estimates given by above are applicable to any general continuous emission density. Now, we derive the formulas for parameter estimation of the underlying continuous observation HMM with Gaussian emission density.
Assuming the Gaussian noise model for emission density, it takes the form
i.e., each (discrete) hidden state corresponds to a single Gaussian density parameterized by mean µ(x t ) and variance Σ(x t ). In such a model, an observation o (m) t corresponding to the m th individual at time t is nothing but a sample from one of the Gaussian densities.
Note that with this aggregate Gaussian-HMM settings, the estimation updates of the initial distribution and transition probabilities remain same as described in Algorithm 4. For learning the model parameters, the required marginals in the E-step of the algorithm are computed using Equation (17) and updates in the M-step are characterized via the following. 
where prime denotes transpose operation.
Proof. See Appendix D.
Based on Proposition 3, the parameters of a Gaussian-HMM are estimated using Algorithm 6.
Remark 3 Convergence of Algorithm 6 follows from Theorem 1.
Remark 4 Similar to Algorithm 5, one can extend the learning of aggregate Gaussian-HMM to the case of an ensemble of continuous observation sequences.
To validate the efficacy of the proposed aggregate learning algorithms, we perform multiple sets of experiments on synthetic as well as real-world datasets.
We begin with synthetic data. We perform multiple sets of experiments for performance comparison of fitted time-invariant HMM models with discrete as well as continuous observations. The initial state probability is sampled from the uniform distribution over the probability simplex. To produce the transition matrix, we firstly randomly permute rows of noised Identity matrix I + 0.05 × √ d × exp(U nif orm[−1, 1]). We scale rows of the permuted matrix so that the resulting matrix is a valid conditional distribution. To evaluate the performance on HMM with discrete observation, the emission matrix is generated in a similar way as transition matrix, but with a different random seed. In case of HMMs with continuous observations, we consider the Gaussian emission model. For each state, the corresponding Gaussian distribution is parameterized by a random mean and variance. The mean is sampled from U nif orm[−5d, 5d] and variance is from U nif orm [1, 5] . In continuous observation setting, the algorithm is required to estimate the initial distribution, the transition matrix and the means of Gaussian emission densities. We generate N individual trajectories from an HMM parameterized with Ψ * . Each aggregate sequence consists of collective observations of M independent trajectories of length T . After learning the parameters by [ N M ] sequences, we generate another N individual trajectories for testing purpose.
We use the negative log likelihood (N LL) as a metric for evaluating performance of our learning algorithm. The difference of N LL between the learned model Ψ and
The model with learned parameters is evaluated on test datasets. Figure 5 shows the performance of our algorithm for different values of state dimension d and population size M on HMMs with discrete observations. Curves in the same figure show learning performance with different values of M but fixed d, T , and N . Figure 6 demonstrates the performance of our algorithm for Gaussian observation models. It shows that our algorithm can effectively learn the generative models. Larger aggregate size corresponds to lower convergence rate, which is intuitive; experiments with larger aggregate size lose more information about individuals. We also observe that in both discrete as well as continuous observation models, the performance of our aggregate learning algorithm degrades as model dimension increases. To further demonstrate the scalability of our algorithm, we conduct experiments with various HMM lengths and sample sizes as depicted in Figure 7 and Figure 8 , respectively. In Figure 7 , the curve in different color depicts the learning performance with various HMM lengths. We observe the larger T leads to better performance. This is because larger T is associated with more training data. Figure 8 demonstrates that the overfitting problem can be eased with more data available in aggregate learning.
Now we evaluate our aggregate learning algorithm in estimating the population movement around Japanese city Tokyo. The dataset 4 consists of anonymous individual trajectories. We discretize the whole city area into 14×16 blocks with each block representing a 15km×15km area. We interpolate individual trajectories every 30 minutes. A total of 6,432, 9,166, 6,822, 10,134, 6,646, 10,338 trajectories are collected respectively on July 1, July 7, October 7, October 13, December 16 and December 29 in the year 2013. We assume that the observation model suffer from a small Gaussian Noise. Moreover, with a small chance, a point in the center block can be categorized to eight neighbouring blocks incorrectly, which account for sensor inaccuracy. In Figure 9 , we show the aggregate observation at timestamps 2:00 and 14:00 generated from the noisy observation model. Our task is to estimate the population flow at timestamps 2 : 00, 8 : 00, 14 : 00, and 20 : 00. At each timestamp, we take the observation generated in one and an half hour and construct a time-invariant HMM graph with M = 20 to infer the movement of the population. Figure 10 presents the comparison between our estimation and ground truth movement at the four timestamps, from which we see that our algorithm successfully recovers the underlying movement of population with noisy aggregate observations.
In this paper, we proposed an algorithm for learning the parameters of a time-homogeneous HMMs from aggregate data. Our algorithm is based on a modified version of the EM algorithm, wherein we utilized the Sinkhorn belief propagation algorithm to infer the unobservable states. In contrast to the existing state-of-the-art algorithms that explicitly consider the aggregate observation noise, our algorithm employs the aggregate observation noise within the graphical model and due to which it is consistent with the standard Baum-Welch algorithm when aggregate data consists of only a single individual. Moreover, our algorithm enjoys convergence guarantees. We further extended our algorithm to incorporate continuous observations and presented estimates for Gaussian observation model. In this work, we have assumed that the HMMs are time-homogeneous, which restricts the modeling capability of the data. We plan to explore learning of time-varying HMMs in our future research.
A Proof of Theorem 1 (Convergence and effectiveness of Algorithm 2) Proof. Proof of convergence The Expectation-Step and Maximization-Step in Algorithm 2 are coordinate ascent updates of −F Bethe (n, y | Ψ) with respect to n and Ψ, and thus the objective function is monotonically increasing. This together with the upper boundedness of the objective function guarantees the convergence of Algorithm 2.
Proof of effectiveness We prove that the log-likelihood L(Ψ) := log p(y|Ψ) monotonically increases approximately. The improvement in log-likelihood at iteration is
.
Since p(n|y, Ψ −1 ) approximately concentrates on n * ,
.
This together with the approximation −F Bethe (y, n | Ψ) of log p(y, n | Ψ) points to
The approximate monotonicity of likelihood then follows from the definition of the M-step in Algorithm 2. Let the Lagrange multipliers be λ, ν, and µ corresponding to the constraints given by (B.1b), (B.1c), and (B.1d), respectively. Then, the Lagrangian can be written as
Now differentiating the Lagrangian w.r.t. the variables and equating to zero, we get ∂ L ∂ π(x 1 ) = ∂ F Bethe (n, y | Ψ) ∂ π(x 1 ) − λ = n 1 (x 1 ) 1 π(x 1 ) − λ = 0,
Solving above Equations, in view of the constraints (B.1b)-(B.1c)-(B.1d), we obtain π(x 1 ) = n 1 (x 1 ), (B.3a)
.
(B.3c) C Proof of Proposition 2 Proof. For a single individual case with M = 1, a sequence of observationsô 1 ,ô 2 , . . . ,ô T is recorded. In such a scenario, the aggregate observations take the form
where δ(·) denotes Dirac function. Then the messages in collective forward-backward algorithm coincide with the messages in standard forward-backward algorithm [32] and take the following form
γ t (x t ) = p(ô t |x t ), (C.2c)
Using above messages, the required marginals can be estimated as n t (x t ) ∝ p(ô t |x t )α t (x t )β t (x t ), (C.3a) n t,t+1 (x t , x t+1 ) ∝ α t (x t )p(x t+1 |x t )β t (x t+1 ) p(ô t |x t )p(ô t+1 |x t+1 ), (C.3b) n t,t (x t ,ô t ) = n t (x t ).
(C.3c)
Then, the parameter update equations given in Algorithm 4 reduce to standard Baum-Welch algorithm. Let the Lagrange multipliers be λ and ν corresponding to the constraints given by (D.2b) and (D.2c), respectively. Then, the Lagrangian can be written as
Now differentiating the Lagrangian w.r.t. the variables and equating to zero, we get
Solving above Equations, we get π(x 1 ) = n 1 (x 1 ), (D.3a)
.
(D.3b)