Learning with Labeled and Unlabeled Multi-Step Transition Data for Recovering Markov Chain from Incomplete Transition Data

1 Introduction

The Markov chain (MC) has been applied to various dynamic systems such as queuing systems, marketing and finance [Neuts, 1981; Pfeifer and Carraway, 2000; Frydman and Schuermann, 2008]. Due to the development and widespread use of sensor devices, MC variants are now used for analyzing the dynamics of urban cities such as traffic, and people flows [Crisostomi et al., 2011; Fan et al., 2015; Iwata et al., 2017].

Since the transition probabilities of MC are unknown in practice, we need to estimate them from observed transition data. In the ideal case, for example, if the information of every visit of each person to each state (location) is available for analyzing people flow in a city 1 , the transition probability can be estimated directly from the number of transitions between the states [Billingsley, 1961]. However, comprehensive data collection is rare because of, for example, constraints imposed by privacy protection, coverage area deficiencies, and limited sensor device precision. Thus in practice the actual data to be analyzed differs from the ideal complete transition data; the transition data available is incomplete as some information is missing.

1 Hereinafter we consider that people's movements are represented by the transition between discretized areas or meshed cells similar to [Fan et al., 2015] for simplicity.

We found that it is essential to be able to analyze two common types of incomplete transition data: hidden transition data where visited state information is partially hidden and dropped transition data where visits to states are partially dropped. See Figure 1. The difference between the two is whether we are aware of the existence of "missing states" or not in the sequence of transitions.

An example of dropped transition data is the transition data of subscribers provided by mobile phone companies. The data contain the transition between states (locations) only where the user stays for more than a certain time to ensure privacy protection and to save memory cost. Therefore, as shown in Fig. 1, although the (true) complete transition consists of two steps, a → b → c, the visit of state b is dropped and only the transition a → c is recorded so the user appears to have visited only state a and c. An example of hidden transition data is the data collected by ourselves using e.g., GPS logger or mobile phone apps. Since the data often contains missing values, the current state (location) information may not be recorded in some logs although the log exists. For example, if the visited state information of state b is missing from a → b → c, only the transition a → ? → c is recorded.

In this paper, we propose a new method that can, in a unified manner, estimate MC parameters from an incomplete transition data set that includes dropped and/or hidden transition data. A key to developing the method is to regard hidden and dropped transition data as labeled and unlabeled multistep transition data, where the label represents the number of steps required for each transition. See the bottom of Fig. 1. The example of the hidden transition data, a → ? → c, and that of dropped transition data, a → c, can be seen as the labeled multi-step transition a → c where the number of steps is 2 and the unlabeled multi-step transition a → c where the number of steps is unknown, respectively. This yields the following description of the generative process of multi-step transition data: the labels indicating the number of steps, k, is first determined and the next state is determined by multistep transition probability, which, according to MC theory, is given by raising the transition probability to the power of k. Figure 2 shows the input and output of the proposed method. By estimating MC parameters, we can recover the (1-step) transition probability that underlies multi-step transitions. It is also shown that the proposed method can be seen as a semisupervised learning method of (possibly infinite) mixture of exponentiated transition probabilities.

We also develop a Majorization Minimization (MM) algorithm [Hunter and Lange, 2004; De Leeuw, 1994] that can guarantee that the objective function improves with each iteration. Note that our formulation also allows us to use complete transition data since it can be represented by labeled multi-step transition data whose labels are all 1. Moreover, our formulation, model, and algorithms are not limited to the case of analyzing people flow and indeed is valid in various scenarios where the available data are represented by labeled and unlabeled multi-step transition data.

The contributions of this paper are summarized below:

• We consider the problem of recovering MC parameters from labeled and/or unlabeled multi-step transition data.
• We develop a new model that can handle labeled and/or unlabeled multi-step transition data; it can be seen as a semi-supervised learning method of (possibly infinite) mixture of exponentiated transition probability.
• We develop an MM algorithm and give proof of its convergence. A method of handling infinite mixtures with theoretical support is also provided.
• We confirm the effectiveness of the proposal by experiments on synthetic and real data sets.

The rest of this paper is organized as follows. §2 details related works and §3 provides a definition and theory of the Markov chain. The proposed method is presented in §4. §5 details the experiments conducted and §6 concludes the paper.

2 RelatedWorks

Recently, new problem formulations for estimating MC parameters have been investigated. For example, [Morimura et al., 2013; Kumar et al., 2015] tackle the problem of recovering MC from a steady state distribution. [Iwata et al., 2017; Akagi et al., 2018] tackle the estimation problem by using snapshots of population data. However, the problem of MC

recovery from multi-step transition data, especially unlabeled multi-step transition data made from dropped transition data, has yet to be studied. Moreover, although handling just hidden transition data is simple since the EM algorithm can be adopted by treating missing state information as latent variables, no unified approach that can handle both dropped and hidden transition data exists. Hidden Markov Model (HMM) cannot be applied to our problem since the time-step of (true) complete transition data and that of dropped transition data do not have one to one correspondence and an infinitely large number of steps may be dropped. This study provides a unified approach to handling the two types of data as well as complete transition data.

Our problem formulation can be seen as semi-supervised learning since the relation between the labeled and unlabeled multi-step data are similar to the way labeled and unlabeled data are treated in semi-supervised clustering [Ghahramani and Jordan, 1994; Nigam et al., 2000; Basu et al., 2002]. However, we emphasize our motivation comes from recovering MC from incomplete transition data.

Our model for MC recovery from multi-step transition data can be seen as a mixture model. However, our model differs from the "standard" mixture of Markov chains that have different component distributions [Goodman, 1961; Frydman, 1984; Gupta et al., 2016]; the k-th component distribution of our model is the transition probability raised to the power of k. This property allows us to approximate the component with large k by a steady state distribution under a mild assumption; it enables us to construct a new type of infinite mixture model that differs from Dirichlet process mixture models [Ferguson, 1973; Antoniak, 1974; Neal, 2000; Blei and Jordan, 2006; Teh, 2010].

3 Preliminaries

Let X = {1, 2, · · · , |X|} be a state space. A discrete time Markov chain (MC) on X is a stochastic process {Xt; t =

0, 1, 2, · · · } that satisfies the following Markov property:

MC M is thus defined by M := {X,P}, where P : X × X → [0, 1] is the transition probability, P(xnext|x) = Pr(Xt+1 = xnext|Xt = x). We also employ adjacency information of Γ = {Γi}i∈X , where Γi is the set of reachable states from state i by one-step transitions.

Later we will use the following theoretical results.

Theorem1. (e.g. Theorem (2.1) [Durrett, 1999]) k-step transition probability is given by P to the power of k, P k .

For example, 2-step and 3-step transition probabilities are given by P 2 and P 3 , respectively. Note that 0-step transition probability is self-transition since a matrix to the power of 0 is an identity matrix. If MCM is irreducible and aperiodic, the following theorem holds.

Theorem 2. (e.g. Theorem (4.5) [Durrett, 1999]) Let π = {πi}|X | i=1 be the steady state probability of Markov chain M. Then, (P k )ij ----→ k→∞ πj.

This states that each row of P k converges to the steady state probability at the limit.

4 ProposedMethod

This section describes the proposed method for estimating Markov chain parameters from multi-step transition data. Figure 2 shows our problem formulation.

4.1 Data

We assume that incomplete transition data are available in the form of unlabeled and labeled multi-step transition data, as shown in Fig 1. We denote the unlabeled multi-step transition data made from dropped transition data (where the number of steps is unknown) as Dum = {Nij}ij∈X where Nij denotes the number of transitions from state i to state j. We also denote the labeled multi-step transition data made from hidden transition data as Dℓm = {Mijk}ij∈X,k∈Z≥0 where Mijk denotes the number of transitions from state i to state j by k-step transitions. If complete transition data is available, we consider that the data are represented by part of the labeled multi-step transition data, {Mij1}ij . We denote the set of the two data types as D = Dum∪Dℓm and refer to it as multi-step transition data. We also use adjacency information Γ as an optional input. In the case of, for example, analysis of people flows in a city, Γ can be derived from maps or street network information. If such adjacency information is not available, we consider Γi to be a set of all states, X. Although Γ is optional, its use is recommended if the labeled multi-step data are unavailable or limited since Γ strongly helps in estimating the number of steps present in the unlabeled multi-step data.

4.2 Model Component

Here we show the model components used by the proposed method.

Label Probability. We decide to use the exponential family for modeling generated labels that indicate the number of steps of each transition since it can express various types of distributions. The density of exponential family f is given by

where η is the natural parameter, T(k) is sufficient statistics and A(η) is the log-normalizer. Examples of distributions belonging to the exponential family are categorical distribution (Cat), geometric distribution (Geo), Poisson distribution (Poi), and zero-truncated Poisson distribution (ZTP):

By assigning specific values to T(x) and A(η), the above distributions are represented by Eq. (1) (See Table 1). We denote the cumulative density function of f as F, where F(C|η) := ∫ C -∞ f(k|η)dk. In a later section, we use meanvalue parameter µ, which is defined as µ := Ef(k|η)[T (k)] = ∂A(η)/∂η. Note that natural parameter η and mean-value parameter µ have a one-to-one correspondence [Amari and Nagaoka, 2007]. Converting µ into η may require numerical computation; we use Newton method for ZTP. We will also use the property that the Hessian of the log-normalizer is positive-definite, since it corresponds to the variance.

Transition Probability. We use parameter ν to model the (1-step) transition probability of a MC. To emphasize parameter dependency, we denote the model of the transition probability as Pν; the (i, j)-th element of matrix Pν represents the (1-step) transition probability from state i to j, i.e., Pr(Xt+1 = j|Xt = i, ν) = (Pν)ij. We denote a Markov chain constructed using the model with parameter ν as Mν = {X, Pν}. Examples of the model are the following tabular model and log-linear model.

Model 1. (Tabular Model) Let us define ν as ν = {{pij}j∈Γi }i∈X where ∑j pij = 1 for all i. The tabular model can be defined as (Pν)ij = pij if j ∈ Γi, and (Pν)ij = 0 otherwise.

Model 2. (Log-LinearModel) Let us define ν as ν = {v,w}. Then, (Pν)ij = exp{vij+φ(i,j)T w} ∑j′∈Γi exp{vij′+φ(i,j ′ )T w} if j ∈ Γi and (Pν)ij = 0 otherwise, where φ(i, j) is a feature vector, such as the (inverse) distance between states i and j 2 .

4.3 Generative Process

The generative process of multi-step transition data can be described using the model component explained in the previous section. We consider that parameters θ = {η, ν} are determined following prior distribution P(θ) = P(η)P(ν). P(η) is a conjugate prior of the exponential family,

where the normalizer Z(ξ0, µ0) is determined by the chosen distribution (See Table 1). As the prior of ν, we can use any distribution. For the case of the log-linear model, a Gaussian distribution is frequently used:

Note that ξ0, µ0 and α0 are hyperparameters 3 .

The stochastic process describing a multi-step transition is determined in the following manner. We denote the process as {X′s ; s = 0, 1, 2, · · · } where s are the time steps of data generation 4 . (i) The initial state, X′0 , is determined by some initial state probability. At every data generation time step s, (ii) the number of steps Ks when moving from current state X′s = i is determined following P(Ks = k|X′s = i) defined using exponential family f,

Then, (iii) next state visited by k step transition is determined following P(X′s+1 |X′s = i,Ks = k) which is defined as, from Theorem 1, Pν to the power of k,

Note that k is a latent variable for unlabeled data and an observed variable for labeled data. The following equations are derived from Eq. (2)(3):

We can interpret this model as being a mixture model because the transition probability defined by Eq. (5) where k is marginalized out is given by the weighted sum of exponentiated transition probability. However, unlike the "standard" mixture model of the Markov chain e.g. [Goodman, 1961;

2 If no such information is available, the term related to the feature and parameters w can be dropped from the model. Note that the tabular model is a special case of the log-linear model since the model without feature vector and parameter w is equivalent to the tabular model if we define pij = exp{vij}/ ∑j′∈Γi exp{vij′}.

3 We set ξ0 = µ0 = α0 = 1.0 in the experiments.

4 Time step s does not correspond to t, the time step of the underlying Markov chain, since the number of (1-step) transitions conducted up to step s′ is given by ∑s′ s=0 Ks.

Frydman, 1984; Gupta et al., 2016], the mixing ratio of our model is given by an exponential family and the k-th component is the transition probability to the power of k. Moreover, when we adopt an exponential family that has unbounded support such as a Poisson distribution, the model becomes an infinite mixture model; a way of handling infinite summation is detailed later.

From Eq. (4)(5), the log-likelihood of the labeled and unlabeled multi-step transition data are given by

The logarithm of joint probability of data D and parameter θ is then given by

Parameter θ is estimated by optimizing this objective.

Note that we can interpret term -logP(θ) as the regularization term; choosing the other prior distribution yields the other type of regularization.

4.4 MajorizationMinimization Algorithm

For minimizing objective function Eq. (6), we use a majorization-minimization (MM) algorithm [Hunter and Lange, 2004; De Leeuw, 1994] that can guarantee that the objective function is improved in each iteration.

The MM scheme indirectly minimizes objective L by using the majorizing function G which is defined as

where Z = {zijk} is an auxiliary variable satisfying zijk ≥ 0 (∀(i, j, k)), ∑k zijk = 1 (∀(i, j)). Majorizing function G has the following two properties:

Note that the equality holds if and only if

MM uses the following 2 step procedures to minimize majorizing function G:

1. Minimize G w.r.t. θ, 2. Minimize G w.r.t. Z. Parameter θ is updated in the following manner.

Update of η. Let ∇θ be the partial derivative operator w.r.t. θ. Setting the partial derivative w.r.t. ηi equal to zero, ∇ηi G(θ, Z) = 0, yields

where µi is the mean value parameter of f(k|ηi). Thus, parameter ηi is obtained by converting µi.

Algorithm 1 Majorization Minimization Algorithm for Markov Chain Recovery from Multi-step Transition Data

Input: D: multi-step transition data, Γ: adjacency information (optional), α0, ξ0, µ0: hyperparameters

Output: ˆ θ = {ˆ η, ˆλ}: estimated parameters

1: Initialize θ.

• 2: repeat
• 3: Update Z following Eq. (8) or (12).
• 4: Update ν by numerical optimization.
• 5: Update η following Eq. (9) or (13).
• 6: until A stopping condition is met

Update of ν. To update ν we need some numerical optimization techniques such as gradient descent. Gradient descent based algorithms update the parameter at optimization step ℓ as follows, νℓ+1 ← νℓ - γℓG-1 ℓ ∇νG(θℓ,Z), where γℓ is the learning rate, Gℓ is the identity matrix for gradient descent, and Hessian is used as the Newton method. ∇ν is the partial derivative operator w.r.t. ν. The partial derivative can be computed as

The experiments described later use the L-BFGS method [Liu and Nocedal, 1989]. The optimization process is summarized in Alg. 1.

We provide here the convergence property of the proposed MM algorithm. To ensure the generality of the analysis to cover the use of various numerical optimization techniques for updating ν, we make the following assumption, similar to the generalized EM algorithm [Dempster et al., 1977].

Assumption 1. Let νold and νnew be parameter ν before and after the update by some (numerical) optimization, respectively. Then, G((η, νnew), Z) ≤ G((η, νold), Z).

Under assumption 1, the following theorem holds.

Theorem 3. The proposed MM algorithm makes objective L(θ) monotonically decreasing. The value is invariant if and only if θ is at a stationary point.

Proof. (Theorem 3) Let us denote the parameters and the auxiliary variables that satisfy L(θ) = G(θ,Z) as θold = (ηold, νold) and Zold, respectively. We also denote η after the first step of the optimization given by Eq. (9) as ηnew and Z after the second step given by Eq. (8) as Znew. Since A(η) is convex, G is also convex w.r.t. η; G((ηnew, ν),Z) ≤ G((η, ν),Z) holds (∀η). From the property of Z shown in Eq. (7), G(θ,Znew) ≤ G(θ,Z) holds (∀Z). Using Assumption 1, L(θold) = G(θold,Zold) ≥ G((ηnew, νold),Zold) ≥ G((ηnew, νnew), Zold)≥G((ηnew, νnew), Znew)=L(θnew).

4.5 Handling Infinite Mixtures

This section shows how to handle infinite mixtures. The key is the use of approximation for computing the infinite summation in P(X′s+1 |X′s ) (Eq. (5)) by adopting the following mild assumption.

Assumption 2. Markov chain Mν is irreducible and aperiodic for any parameter ν.

This assumption allows us to use Theorem 2. By using (sufficiently large) truncation level Ktr, we can approximate P(X′s+1 |X′s ) by ˜ P(X′s+1 |X′s ) as follows:

where πν is the steady state probability of Markov chain Mν. Thus infinite mixture models are handled by using K = Ktr + 1 components, where the final component has mixing ratio 1 - F(Ktr|ηi) and steady state probability πν . Under Assumption 2, the following theorem holds.

Theorem 4. Given that Assumption 2 holds, there exist constants α ∈ (0, 1) and C > 0 such that maxi∈X ‖ ˜ P(X′s+1 |X′s = i) - P(X′s+1 |X′s = i)‖TV ≤ CαKtr , where ‖ · ‖TV is the total variation distance 5 .

Then, the infinite mixture model of Eq. (5) is well approximated by the approximated model of Eq. (11) with sufficiently large Ktr. Proof is shown by using the following theorem.

Theorem 5. [Levin and Peres, 2017] If P is irreducible and aperiodic, with stationary distribution π, there exist constants β ∈ (0, 1) and D > 0 such that maxx∈X ‖P k (·|x)-π‖TV ≤ Dβk .

Proof. (Theorem 4) Since f is a discrete distribution, there exists a finite value fmax := maxk,i f(k|ηi).

The algorithm for (approximate) infinite mixtures is derived in an analogous manner to the finite model.

Update of Z. Similar to Eq. (8),

5 The total variation (TV) distance between probability distribution µ and ν on X is formally defined as ‖µ - ν‖TV = maxA⊂X |µ(A)-ν(A)|. From Proposition 4.2 in [Levin and Peres, 2017], TV is also computed as ‖µ - ν‖TV = 1/2 ∑x∈X |µ(x) - ν(x)|.

Update of η. Similar to Eq. (9),

where ftr(x|η, a, b) is the truncated distribution whose support lies in the range (a, b] 6 and Eftr is the expectation of the distribution

Update of ν. Similar to Eq. (10),

The partial derivative of steady state πν can be computed as 7 ∇νπν = πν∇νPν(I - Pν + 1T 1)-1 , where 1 is a row vector whose elements are all ones.

5 Experiment

This section confirms that the proposed method well handles both unlabeled and labeled multi-step transition data. To evaluate how well the proposed method can recover 1-step transition probability, we use test data given in the form of complete transition data (1-step transition data), see Fig. 1.

5.1 Setting

Synthetic data. In the synthetic data experiment, we set the number of states to 15 and randomly generated chain edges following [Morimura et al., 2009]. The true transition probability was set using the log-linear model where parameters ν ∗ and feature φ were generated using a standard normal distribution. We also added symmetric Dirichlet noise with parameter of 0.3 to the transition probability by taking the weight sum with β = 0.1 for the noise term and 1 - β = 0.9 for the transition probability. The true label probability was set using a Poisson distribution whose parameters {λ∗i } were generated by a gamma distribution with shape = 2.0, scale = 1.0 8 . We generated training and test data sets by generating episodes (sequences of states from initial state) with 20 steps in common. The number of episodes used in generating the test data were 1000. We prepared 5 sets of training and test data.

6 The truncated distribution is often used in survival analysis and is formally defined as ftr(x|η, a, b)={F(b|η)-F(a|η)}-1 f(x|η) if x in (a, b] and ftr(x|η, a, b)=0 otherwise.

7 From πν = 1(I - Pν + 1T 1)-1 as shown in Proposition (2.14.1) [Resnick, 2002].

8λi has mean and variance of 2.0. The 95th percentile of Poisson distribution with λi = 4.0 lies between 7 and 8.

Real data. In the real data experiment, we used car probe data provided by NAVITIME JAPAN Co, Ltd. The dataset is a collection of GPS trajectories of individuals who used a car navigation application in the greater Tokyo area, Japan. We divided the region using an approximately 5km × 5km grid mesh; the total number of mesh cells was approximately 150. We used the data recorded during the period between 2015.4.13 to 2015.4.17 (5 working days in total) in the morning (6:00 am to 10:59 am). The number of unique users per day was, on average, approximately 8000. The data were made by converting the GPS trajectories into sequences of visited mesh cells (states). We excluded the states that appeared fewer than 20 times per day on average and used only episodes containing more than 2 steps. Here we did not use feature φ. Since the data are complete transition data, training data were made by randomly extracting the visited state following the generative process using ZTP, whose parameters are set via a gamma distribution analogous to the synthetic data. Training and test data were made from the logs of one day and that of the next day, respectively. We ran 4 trials using 4 sets of training and test data.

Evaluation Measure. As the evaluation metric, we used the negative test log likelihood computed using the test data given in the form of complete transition data. The negative test log likelihood is defined as (1/Ttest) ∑i,j∈X -˜nij log ˆ pij, where Ttest is the number of total transitions and ˜nij is the number of transitions from state i to state j in the test. ˆ pij is the estimated 1-step transition probability. The results of proposed method using log-linear model for transition probability and Poisson and ZTP for label probability with K=10 are reported 9 . A lower value indicates 1-step transition probability is precisely recovered.

Baseline Methods. Since our method is the first method that can estimate MC parameters from multi-step transition data, we decided to compare the proposed method with the simple baseline methods that use only labeled multi-step transition data. Note that existing methods including HMM and Bayesian nets are not designed to recover 1-step transition probability. Using the log-linear model similar to the proposed method, the parameter of the baseline, ˆ νbase = arg minν {- ∑i,j,k Mijk log{(P k ν )ij}-log P(ν)}, was obtained by numerical optimization. A comparison with this baseline also confirms the usefulness of the proposed method that can handle unlabeled data.

5.2 Results

Quantitative Result. Figure 3 shows the results of the synthetic data experiment. We can confirm that the performance of the proposed method improves as the amount of unlabeled multi-step transition data increases for both cases wherein the labeled multi-step transition data is and is not used; this verifies that the proposed method well handles unlabeled multistep transition data. From Fig. 3b, it is also confirmed that the proposed method outperforms the baseline when challenged

9 In this experiment, "sensitivity of the truncated level K" described later, shows that stable performance is demonstrated when K is set to 10 or more.

(a) w/o Labeled Data (b) w/ Labeled Data (5%) (c) w/ Labeled Data (10%)

Figure 4: Results of real data experiment in which the ratio of episodes used for unlabeled multi-step transition data was varied. The performances (a) without and (b)(c) with labeled multi-step transition data are shown. The labeled multi-step transition data were made from (b) 5% and (c) 10% of the episodes. Lower values are better.

0 1 2 3 4 5 6 7 8 9

10

11

12

13

14

0 1 2 3 4 5 6 7 8 9 1011121314

state-id (xnext)

(a) True

with labeled multi-step transition data. This verifies the usefulness of our method in its ability to handle both labeled and unlabeled data.

Figure 4 shows that similar results were obtained in the real data experiment. Moreover, by comparing the performance of the proposed method with that of the baseline in Fig. 4bc, the degree of improvement is large when the amount of labeled multi-step data is small (5%); this implies that the proposed method is more effective when the amount of labeled data is limited. We can also confirm that the performance of the proposed method improves as the amount of labeled data increases. These results confirm the effectiveness of our method for recovering 1 step transition probability from multi-step transition data.

Sensitivity of the truncation level K. Figure 5 shows the results of the synthetic data experiment when the truncation level (or the number of components) of the proposed method, K, was varied. It shows that the performance improves as the value of K increases, and converges when K is approximately 10. This demonstrate that the proposed method provides stable performance when K is set to 10 or more.

Qualitative Result. Figure 6 illustrates the true and estimated 1-step transition probability from the synthetic data experiment. It confirms that the estimated probability by the proposed method is closer to the true probability. This also implies that the proposed method can accurately estimate MC parameters.

6 Conclusion

This paper proposed a new model and algorithm for estimating MC parameters from labeled and unlabeled multistep transition data for recovering 1-step transition probability from incomplete transition data. We proved the convergence of the algorithm and introduced an approximation scheme for handling infinite mixture models that has strong theoretical support. We also evaluated the effectiveness of the proposed method using synthetic and real data. Remaining future work is to theoretically analyze the performance of the proposed method. Developing Bayesian algorithms is also a promising direction.