Estimating People Flow from Spatio-temporal Population
Data via Collective Graphical Mixture Models
TOMOHARU IWATA, HITOSHI SHIMIZU, FUTOSHI NAYA, and NAONORI UEDA, NTT
Communication Science Laboratories
Thanks to the prevalence of mobile phones and GPS devices, spatio-temporal population data can be obtained
easily. In this paper, we propose a mixture of collective graphical models for estimating people flow fromspatio-temporal population data. The spatio-temporal population data we use as input are the number ofpeople in each grid cell area over time, which are aggregated information about many individuals; to preserve
privacy, they do not contain trajectories of each individual. Therefore, it is impossible to directly estimate
people flow. To overcome this problem, the proposed model assumes that transition populations are hidden
variables, and estimates the hidden transition populations and transition probabilities simultaneously. The
proposed model can handle changes of people flow over time by segmenting time-of-day points into multiple
clusters, where different clusters have different flow patterns. We develop an efficient variational Bayesianinference procedure for the collective graphical mixture model. In our experiments, the effectiveness of the
proposed method is demonstrated by using four real-world spatio-temporal population data sets in Tokyo,
Osaka, Nagoya and Beijing.
CCS Concepts: • Information systems → Data mining;
Additional Key Words and Phrases: collective graphical models, mixture models, variational Bayes, spatiotemporal data, population data
ACM Reference format:
Tomoharu Iwata, Hitoshi Shimizu, Futoshi Naya, and Naonori Ueda. 2017. Estimating People Flow from
Spatio-temporal Population Data via Collective Graphical Mixture Models. ACM Trans. Spatial Algorithms
Syst. 9, 4, Article 39 (March 2017), 18 pages.
DOI: 0000001.0000001

1

INTRODUCTION

Analyzing spatio-temporal population data is important in various fields including disaster management [20], marketing [5], public health [3], and urban planning [24]. For example, population
information is useful for planning the locations of new stores. Census data have long been usedas spatial population data. However, since conducting a census is time-consuming and costly,
census data are not updated frequently; e.g. the US census is conducted every ten years. In recent
years, thanks to the prevalence of mobile phones and global positioning system (GPS) devices,spatio-temporal population data can be obtained easily. For example, mobile spatial statistics [22]contain the hourly population in 500 meter grid squares in Japan’s urban areas, which are calculated
A part of this work has been achieved by “Research and Development on Fundamental and Utilization Technologies forSocial Big Data,” the Commissioned Research of National Institute of Information and Communications Technology (NICT),JAPAN.
Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee
provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and thefull citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored.
Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires
prior specific permission and/or a fee. Request permissions from permissions@acm.org.
© 2017 Copyright held by the owner/author(s). Publication rights licensed to ACM. 2374-0353/2017/3-ART39 $15.00
DOI: 0000001.0000001
ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39

39:2

T. Iwata et al.

from mobile network operational data. Similar data can be obtained from the GPS data of navigation
apps [23], Wi-Fi [13], geotagged tweets 1 or photos 2 and check-in data of location sharing services
such as foursquare 3 .
In this paper, we consider the task of estimating people flow from spatio-temporal population data.
The estimated people flow can be used for a wide variety of applications, which include simulating
people movement in the case of a disaster, detecting anomalous people movement, predicting
the future spatial population given the current spatial population, and designing transportation
systems. The spatio-temporal population data we use as input are the populations in each gridcell over time as shown in Figure 4. The output is the people flow between grid cells over time
as shown in Figure 7. The spatio-temporal population data consist of aggregated information
about many individuals. They are aggregated to preserve privacy or because of the difficulty of
tracking individuals over time. For instance, the mobile spatial statistics [22] are preprocessedfor privacy protection and no one can follow a particular user, which allows the mobile phone
operating company to publish population data calculated based on information about 60 million
mobile phone users. If the trajectories for each individual are given, people flow can be estimatedstraightforwardly by counting the number of people who moved between grid cells, i.e. transition
population. However, with aggregated data, it is impossible to directly determine the size of thetransition population.
To overcome this problem of modeling individual behavior given aggregated data, we propose a
mixture of collective graphical models. The proposed model assumes that individuals move according
to transition probabilities that depend on their locations and time points. Since the transition
populations are not given, we treat them as hidden variables. The hidden transition populations
relate to observed populations in each cell; the population at a cell is equal to the sum of transition
populations from the cell, and the population at a cell in the next time point is equal to the sumof transition populations to the cell. By using these relations that represent flow conservation as
constraints, the hidden transition populations and transition probabilities are inferred simultaneously. The proposed model can handle changes of people flow over time by segmenting time-of-day
points into multiple clusters, where different clusters have different transition probabilities. For ourtask of modeling people flow, incorporating time information is crucial. For example, people movefrom suburbs to the city center to work in the morning, they return to the suburbs from the city
center in the evening after work, and they do not travel so often in the middle of the night. The
proposed model is an extension of the collective graphical models [17, 18] for handling change overtime. Since the proposed model does not use trajectory information, when we have enough amountof population data that approximate the true population distribution, we can estimate people flow.
We develop an efficient variational Bayesian inference procedure for the proposed model. In
existing literature on collective graphical models, the EM algorithm is used to obtain point estimatesof the parameters [6, 10, 17, 21]. Instead of obtaining point estimates, by estimating distributions ofparameters based on the variational Bayesian framework, we can alleviate overfitting especiallywhen data are sparse and models are flexible as in the case of mixture models [1].
The major contributions of this paper include the following:• It is the first attempt to estimate people flow from spatio-temporal population data without
tracking.
• We propose a mixture of collective graphical models for handling behavior change overtime.
1 http://twitter.com
2 http://www.flickr.com
3 http://foursquare.com

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

Collective Graphical Mixture Models

39:3
Table 1. Notation

Symbol
N ti
M ti jT
IJi
Ss (t )дsr (i, j)Kz s (t )θ ki j
ϕk
τk
ηk
θ ki
N
M

Description#people at grid cell i at time t, N ti ≥ 0#people who moved from grid cell i to grid cell j at time t, M ti j ≥ 0#time points, t ∈ {1, · · · ,T }#grid cells, i ∈ {1, · · · ,I }neighbors of grid cell i#time-of-day indicestime-of-day index of time t, s (t ) ∈ {1, · · · ,S }sth time-of-day дs ∈ [0, 24]
relative position index of grid cell j from grid cell i#clusters
cluster index of time t, k ∈ {1, · · · ,K }probability that a person in grid cell i moves to grid cell j during a time-of-day point
P
that belongs to cluster k, θ ki j ≥ 0, j ∈J i θ ki j = 1
PKprobability that a time-of-day index belongs to cluster k, ϕ k ≥ 0, k=1
ϕk = 1time-of-day mean of cluster k, τ k ∈ [0, 24]
time-of-day precision of cluster k, η k > 0
set of transition probabilities at grid cell i during a time-of-day point that belongs tocluster k, θ ki = (θ ki j ) j ∈J iI )T
set of population data over time and grid cells, N = ((N ti )i=1
t=1
I )
T −1
set of transition population data over time and grid cells, M = (((M ti j )i=1
j ∈J i ) t=1• We develop a variational Bayesian inference procedure for collective graphical models.
The paper is organized as follows: In Section 2, we propose collective graphical mixture models for
estimating people flow from spatio-temporal population data, and present its variational Bayesianinference procedure in Section 3. In Section 4, we demonstrate the effectiveness of the proposed
method by using real spatio-temporal population data obtained in Tokyo, Osaka, Nagoya andBeijing. In Section 5, we outline related work. Finally, we present concluding remarks and future
work in Section 6.
2

PROPOSED MODEL

I ) T , where
Suppose that we have population data over time for each of I grid cells, N = ((N ti )i=1
t=1
N ti is the number of people in grid cell i at time t. We also have neighbor information for eachgrid cell. Let Ji ⊆ {1, · · · ,I } be the set of neighbor cells of grid cell i, and Ji always contains
grid cell i itself, i ∈ Ji . We assume that people do not move to its non-neighbor cells, or we canignore those people because their size is small. Our task is, given aggregated population data
I ) T , to estimate transition population M = (((M ) I )T −1
N = ((N ti )i=1
ti j i=1 j ∈J i ) t=1 in the time period
t=1
of the given aggregated data, where M ti j is the number of people who move from grid cell i togrid cell j at time t. Our notation is summarized in Table 1. We use non-bold typefaces for scalarvariables, and bold typefaces for vector and set variables. The time index t and time-of-day indexs (t ) can be different since the given data might be taken at multiple days. Suppose that the tth time
point is Oct 1st 10:00, and t 0th time point is Oct 3rd 10:00. Since these time points have the sametime-of-day, their time-of-day indices are the same s (t ) = s (t 0 ), but t , t 0.
We model the probability of people moving to the neighbor grid cells using a mixture of collective graphical models. With the proposed model, the transition patterns are assumed to change

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:4

T. Iwata et al.

depending on the time-of-day, and the time-of-day points are divided into K clusters according totheir transition patterns. Let θ ki j be the probability that a person in grid cell i moves to grid cell jduring a time-of-day point that belongs to cluster k. The relationship between observed data N andtransition probability θ ki j can be modeled by introducing latent variable M ti j , which represents
the number of people who move from grid cell i to grid cell j at time t. In particular, when thetime-of-day of t is assigned to cluster k, the probability of transition population Mti = (M ti j ) j ∈J i
given the current population N ti and transition probability θ ki = (θ ki j ) j ∈J i is given by the following
multinomial distribution,
Y M
N ti !
p(Mti |N ti ,θ ki ) = Q
θ ki jti j .
(1)j ∈J i M ti j ! j ∈Ji

The population in a cell N ti and the transition population between cells Mti have the following
two relations,X
N ti =
M ti j ,(2)
j ∈J iwhich represents the fact that the population at grid cell i is the same as the sum of the transition
populations from grid cell i, andX
N t+1,i =M t ji ,(3)
j ∈J iwhich represents the population in grid cell i at the next time point is the same as the sum of thetransition populations to grid cell i at the current time point. These relations are used as constraints
in the inference as described in Section 3.
For the prior distribution of the transition probability, we use the non-symmetric Dirichlet
distribution. With the proposed model, the hyperparameters of the Dirichlet distributions areshared among different grid cells when the relative positions are the same. Let r (i, j) be the index
of the relative position of grid cell j from grid cell i. Figure 1 shows an example of the relative
position index when the neighbors are its surrounding eight cells with the addition of the cell itself.
Then, the Dirichlet prior for the transition probaiblity is given as follows,Y αΓ(ᾱ i )−1p(θ ki |α i ) = Q
θ kirj(i, j) ,
(4)
Γ(α
)
r
(i,
j)
j ∈J ij ∈J iP
where α i = (α r (i, j) ) j ∈J i is the Dirichlet hyperparameter, and ᾱ i = j ∈J i α r (i,j) . The parameter α r
represents the prospect of moving in the r th relative position, which is shared by all cells andall time points, e.g. ᾱα 1i is the prior probability of staying at the current point when r (i,i) = 1. By
estimating the Dirichlet hyperparameters α = (α 1 ,α 2 , · · · ) from the given data, we can robustly
calculate the population flow especially in cells with sparse data by sharing hyperparameters among
different cells.
People flows at temporally close time points are considered to be similar. For example, transition
probabilities at 14:00 are similar to those at 15:00. To model this property, we assume that the
time-of-day points in cluster k are distributed as Gaussian with mean τ k and precision η k , where
precision is the inverse of variance, as follows, η  12 ηkkp(дs |τ k ,η k ) =exp − |дs − τ k | 2 ,(5)2π2where дs is the sth time-of-day. For example, when observations are given for every 30 minutes,
д1 = 0, д2 = 0.5, д48 = 23.5, and the total number of time-of-day indices is S = 48. We use continuoustime-of-day дs for clustering time-of-day points based on the continuous time-of-day values, and
ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

Collective Graphical Mixture Models

a

b

c

d

e

f

g

h

i

39:5
r(e,a)=9 r(e,b)=2 r(e,c)=3r(e,d)=8r(e,e)=1r(e,f)=4r(e,g)=7 r(e,h)=6(a) name for each grid cellr(e,i)=5(b) relative position indices from the center grid cellFig. 1. (a) name for each grid cell and (b) their relative position indices from the center grid cell “e”, r (e, ·).
cluster3
cluster4

cluster2

0:00

4:00

cluster1

8:00

12:00

16:00

20:00

24:00

Fig. 2. Example of Gaussian distributions over time-of-day. The mean of cluster k is τ k , and the precision isdη k .

use time-of-day index s for representing the point to be clustered. Figure 2 shows an example ofGaussian distributions over time-of-day. By assuming the Gaussian distribution for each cluster, we
can cluster time-of-day points depending on their interval of time-of-day as well as their transition
probabilities.
The proposed model assumes the following generative process for the transition population overI ) T −1 ,time M = ((Mti )i=1
t=1
(1) Draw cluster proportions, ϕ ∼ Dirichlet(β )(2) For cluster k = 1 to K(a) For grid cell i = 1 to I(i) Draw transition probability, θ ki ∼ Dirichlet(α i )(b) Draw precision of time-of-day, η k ∼ Gamma(a,b)(c) Draw mean of time-of-day, τ k ∼ Normal( f , (dη k ) −1 )(3) For time-of-day index s = 1 to S(a) Draw cluster, z s ∼ Discrete(ϕ)(b) Draw time-of-day, дs ∼ Normal(τ z s ,η z−1s )(4) For time points t = 1 to T − 1(a) For grid cell i = 1 to I(i) Draw transition population, Mti ∼ Multinomial(θ z s (t ) i , N ti )K is cluster proportions, and s (t ) represents the index of time-of-day at t. We
Here, ϕ = (ϕ k )k=1
assume that the transition probability is the same when time-of-day is the same even if theirdays are different. We use conjugate priors, i.e. Gaussian-Gamma priors for the Gaussian meanand precision, and Dirichlet priors for Discrete or multinomial parameters, which enables us toanalytically derive equations for updating model parameters in the variational Bayesian inference.
In particular, the prior of Gaussian mean τ k and precision η k is Gaussian-Gamma distribution

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:6

T. Iwata et al.

K

α
a

M

θ

N

I

s

η

T

b
τ

f

g

z

φ

β

S

d

Fig. 3. Graphical model representation of the proposed mixture of collective graphical models.

Normal( f , (dη k ) −1 )Gamma(a,b), and the prior of multinomial parameters θ ki and ϕ is Dirichlet
distribution Dirichlet(α i ) and Dirichlet(β ), respectively, where f , d, a, b, α i and β are hyperparameters. Given spatio-temporal population N, time-of-day indices s and hyperparameters α , β,a,b,d, f ,S , model parameters
the joint probability of transition population M, time-of-day data g = (дs )s=1
I
K
K
K
S is given asΘ = ((θ ki )i=1 )k=1 , ϕ, τ = (τ k )k=1 , η = (η k )k=1 , and cluster assignments z = (z s )s=1
follows,p(M, g, Θ,ϕ, z,τ ,η|N, s,α , β,a,b,d, f )= p(M|N, Θ, z, s)p(Θ|α )p(z|ϕ)p(ϕ|β )p(g|z,τ ,η)p(τ | f ,d,η)p(η|a,b).

(6)

Figure 3 shows a graphical model representation of the proposed model, where the shaded andunshaded nodes indicate observed and latent variables, respectively. We consider that the time-ofday index s is observed since it is determined from the time of the population data.
3

INFERENCE

We present a variational Bayesian inference procedure for the proposed collective graphical mixture
model. The observed variables are spatio-temporal population N and time-of-day data g. The
unknown variables are transition probabilities Θ, cluster assignments z, cluster proportions ϕ,time mean τ , time precision η, transition populations M, and transition hyperparameters α . We
approximate the posterior distributions of all unknown variables except for M and α . We obtain apoint estimate of M and α by maximizing the lower bound of the log marginal likelihood sinceapproximate posterior distributions of M and α cannot be analytically obtained.
We approximate the following true posterior distributionsp(Θ,ϕ, z,τ ,η|M, N, g, s,α , β,a,b,d, f ),(7)

by using the mean-field familyq(Θ,ϕ, z,τ ,η) =K Y
IY
k=1 i=1q(θ ki )q(ϕ)

S
Y
s=1

q(z s )KYq(τ k ,η k ),

(8)

k=1

which are called variational distributions. The following lower bound of the log marginal likelihood
is obtained by applying Jensen’s inequality,
log p(M|N, s,α , β,a,b,d, f ) ≥ E[log p(M|N, Θ, z, s)] + E[log p(Θ|α )] + E[log p(z|ϕ)]+ E[log p(ϕ|β )] + E[log p(g|z,τ ,η)] + E[log p(τ | f ,d,η)] + E[log p(η|a,b)]+ H[q(Θ)] + H[q(ϕ)] + H[q(z)] + H[q(τ )] + H[q(η)] ≡ F ,ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

(9)Collective Graphical Mixture Models

39:7

where E[f ] represents the expectation of f with variational distributions q, and H[q] represents
the entropy of q. The complete form of the lower bound is described in Appendix A.
We iteratively optimize each of the variational distributions and unknown variables while keeping
the others fixed by maximizing the lower bound (9). Since we use conjugate priors, the variational
distributions to be estimated have the same form as their priors as follows,0q(θ ki ) = Dirichlet(α ki
),

(10)

q(ϕ) = Dirichlet(β 0 ),(11)q(τ k ) = Normal( f k0, (d k0 η k ) −1 ),(12)q(η k ) = Gamma(a k0 ,b k0 ),(13)

0 = (α 0 )K . Letwhere α ki
, and β 0 = (β k0 )k=1
ki j j ∈J i

q sk ≡ q(z s = k ),(14)

be the variational probability that the cluster assignment of time-of-day s is k. Its update equation
is given bylog q sk ∝ Ψ(β k0 ) − Ψ( β¯0 ) +I X
X
i=1 j ∈J i

00
M ti j [Ψ(α kij ) − Ψ(ᾱ ki )]

a0
1
1
1
(15)
+ Ψ(a k0 ) − log b k0 − 0 − k0 (дs − f k0 ) 2 ,2
2
2d k 2b k∂ log Γ(x )0 =P0¯0 P K β 0 .
where Ψ(·) is the digamma function defined by Ψ(x ) = ∂x , ᾱ kij ∈J i α ki j , and β =k=1 k
The update equations of unknown parameters in the variational distributions (10)–(13) are given by
T
−1
X

0
α ki
j = α r (i, j) +q s (t )k M ti j ,

(16)

t=1

β k0 = β +d k0 = d +SX

q sk ,(17)

q sk ,

(18)

s=1
S
X
s=1

f k0 = d k0−1 (d f +SX

q sk дs ),

(19)

s=1

PS
a k0 = a +s=1 q sk + 1

2

,

(20)

S

b k0 = b +1Xdq sk (дs − f k0 ) 2 + ( f k0 − f ) 2 .
2 s=1
2

(21)

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:8

T. Iwata et al.

We obtain a point estimate of transition populations M by maximizing the lower bound (9). In
the lower bound, only the first term E[log p(M|N, Θ, z, s)] (30) depends on M, whose approximation
is given by
Z X
E[log p(M|N, Θ, z, s)] =q(Θ)q(z) log p(M|N, Θ, z, s)dΘ
z

∝

T
−1 X
I X
X

K
X

t=1 i=1 j ∈J i

k=1≈

−log M ti j !+M ti j

T
−1 X
I X
X

0
0
q s (t )k [Ψ(α kij ) −Ψ(ᾱ ki )]M ti j − M ti j log M ti j + M ti j

t=1 i=1 j ∈J i

K
X
k=1



0
0
q s (t )k [Ψ(α kij ) − Ψ(ᾱ ki )] ≡ L,

(22)

where log N ti ! is omitted from the second line since it does not depend on M, and Stirling’s
approximation, log M! ≈ M log M − M, is used at the third line as in [17]. The constraints (2) and (3)might not hold in real-world data because errors in sensing are inevitable and populations would
be modified to preserve privacy. We incorporate them as soft constraints, and try to minimize thesquared difference between the left- and right-hand sides in each of (2) and (3). Then, the objectivefunction to be maximized becomes as follows,
L0 = L −T −1 IT −1 I
X
X
λ XX
λ XX
|N ti −
|N t+1,i −M ti j | 2 −M t ji | 2 ,2 t=1 i=1
2 t=1 i=1
j ∈J
j ∈J
i

(23)

i

where the second and third terms in the left-hand side correspond to the soft constraints of (2)and (3), respectively, and λ ≥ 0 is a hyperparameter for controlling the penalty for violating
the constraints. Since integer programming problems require a prohibitively large amount ofcomputational time, we allow integer-valued variable M to take any real value, which enables us tomaximize the objective function efficiently. We maximize the objective function L 0 with respect to
M by using the quasi-Newton method [9] with non-negative bound constraint M ≥ 0. The gradientof the objective function L 0 is given by
K
X
X
X
∂L 0
0
0
)] + λ(N ti −M ti j ) + λ(N t+1,i −M t ji ).
q s (t )k [Ψ(α ki) − Ψ(ᾱ ki= − log M ti j +j∂M ti j
j ∈J
j ∈J
k=1

i

i

(24)
The objective function L 0 is convex with respect to M since the first and third terms in (22) arelinear, the second term in (22), M ti j log M ti j , is convex, and the second and third terms in (23) areconvex quadratic functions. Therefore, the global optimum solution for M is obtained when theother parameters are fixed. Note that Mt can be optimized in parallel since the objective function
L 0 is decomposed into the sum of T − 1 terms each of which contains only the transition population
at t, Mt .
We estimate the transition hyperparameters α by maximizing the lower bound (9) using the
fixed-point iteration method described in [12]. The update rule is given by
PK PI P0
k=1 i=1 j ∈J i I(r (i, j) = r )[Ψ(α ki j ) − Ψ(α r )]α r ← α r PK PI P,
(25)
0
k=1 i=1 j ∈J i I(r (i, j) = r )[Ψ(ᾱ ki ) − Ψ(ᾱ )]where I(A) is the indicator function, i.e. I(A) = 1 is A is true and I(A) = 0 otherwise.
We can obtain a local optimum for the model parameters. Algorithm 1 summarizes the variational Bayesian inference procedure for the proposed collective graphical mixture model. The end
ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

Collective Graphical Mixture Models

39:9

ALGORITHM 1: Variational Bayesian inference procedure for the proposed collective graphical mixture
model.
I , the number of clusters K,
Input: spatio-temporal population data N, neighbor information (Ji )i=1
hyperparameters a,b, f ,d, β,λ0 ) I ) K , β 0 , ( f 0 ) K , (d 0 ) K , (a 0 ) K , (b 0 ) K ,Output: parameters of variational distributions ((α ki
i=1 k=1
k k=1
k k=1
k k=1
k k=1
K ) S , transition population M, transition hyperparameters α((q sk )k=1
s=1
K ) S , M, α ;initialize ((q sk )k=1
s=1
repeat0 by (16) for k = 1 to K, i = 1 to I , j ∈ J ;update α kiijupdate β k0 by (17) for k = 1 to K;
update d k0 by (18) for k = 1 to K;
update f k0 by (19) for k = 1 to K;
update a k0 by (20) for k = 1 to K;
update b k0 by (21) for k = 1 to K;
update q sk by (15) for s = 1 to S, k = 1 to K;
update M by maximizing (24);update α r by (25) for all relative position indices;
until;

condition can be based on the number of iterations, or the convergence of the following valueF−T −1 IT −1 I
X
X
λ XX
λ XX
|N ti −
M ti j | 2 −|N t+1,i −M t ji | 2 ,2 t=1 i=12
t=1 i=1
j ∈J
j ∈J
i

(26)

i

which is the sum of the lower bound of the log marginal likelihood (9) and soft constraints thatare used in (24). In the experiments, we fixed the hyperparameters as follows: β = 10 −2 ,a = 1,b =1, f = 12,d = 1. It is impossible to tune optimal soft constraint parameter λ for estimating flow,
since our task is unsupervised and any people flow data are not given. Therefore, λ is set by using
the predictive performance of the population of the next time point in the training data for eachdata set and for each number of clusters.
The computational complexity for an iteration with the proposed method is O (T I JK ); it increases
linearly with respect to the number of time points T , the number of grid cells I , the numberof neighbors J , and the number of clusters K. Note that it does not depend on the population
P PI
size Tt=1 i=1
N ti since by using the framework of collective graphical models and continuous
approximation of population M, we only deal with the sufficient statistics instead of individual
behavior data. Therefore, the proposed method is applicable to large quantities of population data.
0 and qThe transition probability at grid cell i at time t is estimated using estimated α ki
sk by
θˆti j =α ẑ0 s (t ) i jᾱ ẑ0,(27)s (t ) i

where ẑ s (t ) = arg maxk q s (t )k is the estimated cluster assignment of time t. The expected number ofpeople who move from grid cell i to grid cell j at time t is calculated by θˆti j N ti . Then, the populationof the next time point given the current population Nti is predicted byXN̂ t+1,i =θˆti j N ti ,(28)
j ∈J i

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:10T. Iwata et al.

since the sum of the transition population to i is equal to the population at i at the next timestep.
4
4.1

EXPERIMENTS
Data

We evaluated the proposed method using real-world spatio-temporal population data sets obtainedin Tokyo, Osaka, Nagoya and Beijing. Figure 4 shows examples of gridded population at 0:00, 6:00,
12:00 and 18:00 on a day in Tokyo, Osaka, Nagoya and Beijing data. With all of the data sets, gridcells are square, and the neighbors are its surrounding eight cells with the addition of the cell itself.
For the grid cells at corners or edges, the outside cells are removed from their neighbor sets.
The original data of the Tokyo, Osaka and Nagoya data sets 4 contained latitude and longitude
information for each person obtained every five minutes that were interpolated using railway androad information [15] from geotagged tweets. We created gridded spatio-temporal population data,
where the time interval was 30 minutes, the size of each grid cell was 10km × 10km, and there were
16 × 14 grid cells. The Tokyo data consisted of data on July 1st, July 7th, October 7th, October 13th,
December 16th, and December 22nd of 2013, where the numbers of people were 6,432, 9,166, 6,822,
10,134, 6,646 and 10,338, respectively. The Osaka data consisted of data on August 8th, August 11th,
September 16th, September 22th, December 24th, and December 29nd of 2013, where the numberof people were 2,256, 3,034, 2,999, 3,569, 2,487 and 3,480 respectively. The Nagoya data consisted ofdata on July 22th, July 28th, September 16th, September 22th, December 24th, and December 29ndof 2013, where the numbers of people were 929, 1,332, 1,148, 1,460, 975 and 1,570, respectively.
The Beijing data consist of population in Beijing from 3rd February to 7th of 2008, which wereobtained from T-Drive trajectory data [25, 26], which contained trajectories of 10,357 taxies. We
created gridded spatio-temporal population data and with 15-minute time intervals, 2km × 2km
size grid, and 20 × 16 grid cells in total. The length of the area was 40km in the north-south direction,and 32km in the east-west direction.
We used the data sets whose original data contain trajectories of twitter users or taxies so that
we can evaluate the people flow estimation performance. We did not use trajectory information for
the estimation.
4.2

Results

We evaluated the proposed variational Bayesian collective graphical mixture model (VCGMM)
in terms of the estimation performance of the people flow. For comparing methods, we used thefollowing four methods: CGM, ICGM, VCGM and STAY. The CGM is a collective graphical modelobtained by a maximum likelihood estimation with the same soft constraints with the proposedmodel, where the transition probability for each grid cell is fixed over time. The ICGM is theinhomogeneous transition probability CGM, which assumes that the transition probability for eachgrid cell is different across different time points. The VCGM is a collective graphical model obtainedwith a variational Bayesian inference. When there is one mixture, the proposed method corresponds
to the VCGM. CGM and VCGM assume that people flows do not change over time, while VCGMM
and ICGM assume that they change over time. The STAY method predicts the population at thenext time point from the current population; it assumes that all people do not move, and stay intheir current cells. With the proposed VCGMM, we used ten clusters.

4 SNS-based People Flow Data, Nightley, Inc., Shibasaki & Sekimoto Laboratory, the University of Tokyo, Micro Geo Data

Forum, People Flow project, and Center for Spatial Information Science at the University of Tokyo, http://nightley.jp/archives/1954
ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

Collective Graphical Mixture Models
0:00

39:11

6:00

12:00

18:00

(a) Tokyo(b) Osaka(c) Nagoya(d) Beijing
Fig. 4. Spatio-temporal population data (a) in Tokyo on July 1st 2013, (b) in Osaka August 8th 2013, (c) in
Nagoya on July 22th 2013, and (d) in Beijing on February 3rd 2008. Darker colors represent higher population
densities in each grid cell.

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:12

T. Iwata et al.
Table 2. Normalized absolute errors when predicting people flow averaged over all time points.

Tokyo
Osaka
NagoyaBeijing

VCGMM
0.167
0.250
0.250
0.470

ICGM
0.176
0.265
0.281
0.500

VCGM
0.208
0.280
0.292
0.479

CGM
0.208
0.280
0.291
0.479

STAY
0.192
0.272
0.269
0.532

For a evaluation measurement, we used the following normalized absolute error on people flow,
P I PT −1 P∗
j ∈J i |M ti j − M ti j |i=1 t=1,(29)
P I PT −1
i=1 t=1 N ti
∗ is the true number of transition people who moved from grid cell i to the jth neighbor atwhere M ti
jtime t, and M ti j is its estimation. The true number of transition people, which is the test data, wasobtained from the original data that contain individual trajectories. Note that for estimation we used
population data for each grid cell and each time points, but did not use transition population data
between grid cells. Table 2 shows the error of VCGMM, ICGM, VCGM, CGM and STAY methods.
The proposed method, VCGMM, achieved the lowest error, which implies its effectiveness on people
flow estimation. Since the ICGM estimates the transition probability using only the data at thattime point, its performance was worse than the proposed method. The errors by the VCGM andCGM were almost the same. The VCGM and CGM were better than the STAY with the Beijing data,
but worse with the Tokyo, Osaka and Nagoya data. It is because that many people stay in the same
cell with the Tokyo, Osaka and Nagoya data, and many people move to different cells with the
Beijing data, which is reasonable since the Beijing data is taxi trajectories. The proposed method
achieved the best performance in both types of data by flexibly changing flows depending on time
points.
Figure 5 shows the normalized absolute errors with different numbers of clusters with theproposed VCGMM. As the number of clusters increased, the error decreased. This result indicatesthe importance of using multiple people flow patterns. The error became relatively steady after
seven, three and three with the Tokyo, Osaka and Nagoya data, respectively, and the VCGMM did
not overfit the training data even with many clusters.
Figure 6 shows the estimated cluster proportions q sk provided by the VCGMM with ten clusters.
For each time-of-day index, only one cluster proportion was estimated as nearly one, and theother cluster proportions were estimated as nearly zero. Consecutive time points were assigned
to the same cluster. Nine, seven, six and nine clusters were used with the Tokyo, Osaka, Nagoya
and Beijing data sets, respectively. With the variational Bayesian inference, the proportions forunnecessary clusters become zero, and the proportions for necessary clusters become non-zero,
which helps to avoid overfitting.
Figure 7 shows people flows estimated by the proposed method. With the Tokyo data, From 0:00
to 6:00, there were few flows; this is reasonable since most people are sleeping at home at this time.
From 8:30 to 11:00, people move to the center of Tokyo from the suburbs for work. There were stillflows to the city center from 14:30 to 15:00. From 18:30 to 23:00, flows from the city center to thesuburbs were estimated since at that time people return home from their offices. With the Osakadata, similar temporal flow patterns were extracted. With the Nagoya data, flows to the city center
were discovered in the day time. Since its population size was smaller than that of the Tokyo andOsaka data, estimated flows were not clear. With the Beijing data, the proposed method estimatedthat people move to the center in the morning.

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

Collective Graphical Mixture Models
0.21

39:13
0.280

VCGMM
CGM
ICGM
STAY

0.20

0.270

0.19

0.265
0.260

0.18

0.255
0.250

0.17
1

VCGMM
CGM
ICGM
STAY

0.275

2

3

4

5

6

#clusters

7

8

9

10

0.24512

3(a) Tokyo

5

6#clusters

7

8

9

10

(b) Osaka

0.30

0.54

VCGMM
CGM
ICGM
STAY

0.29
0.28

VCGMM
CGM
ICGM
STAY

0.53
0.52
0.51

0.27

0.50
0.49

0.26

0.48

0.25
0.241

4

0.47
2

3

4

5

6

#clusters(c) Nagoya

7

8

9

10

0.461

2

3

4

5

6

#clusters

7

8

9

10

(d) Beijing

Fig. 5. Normalized absolute errors with different numbers of the clusters with the proposed VCGMM.

5 RELATED WORKGreat interest is being shown in developing methods to analyze people trajectories for a wide
variety of applications, such as disaster management [19, 20], marketing [5], public health [2, 3],urban planning [14, 24, 27], traffic forecasting [8], traffic anomaly detection [11], transportation
system management [4], and travel route recommendation [7]. These methods require the trajectories of individuals. However, because such trajectories are private information, or because
tracking individuals is difficult, only aggregated population data might be available as mobile spatial
statistics [22]. The proposed method can estimate people flows using aggregated data without
trajectory data, and it would make existing trajectory data mining methods applicable even when
trajectory data are unavailable.
Collective graphical models have been proposed and used for modeling contingency tables [18],bird migration [16], and infection [6], but have not been applied to population data. The proposedmodel is based on collective flow diffusion models [6]. We extend the collective flow diffusion
models to mixture models for handling changes in flow patterns over time. To our knowledge,
the proposed model is the first mixture model consisting of collective graphical models. To inferhidden variables in collective graphical models, several algorithms have been proposed, such as
Markov chain Monte Carlo (MCMC) [18], maximum a posteriori (MAP) [17], the expectation
propagation with Gaussian approximation [10], and the belief propagation [21]. The variational
Bayesian inference procedure presented in this paper is applicable not only to a mixture of collective
graphical models, but also to other collective graphical models.
The CGMs in [10, 17, 18, 21] consider noisy observations, where a noise model, such as the Poission distribution, is assumed. The proposed model does not consider noise models for observations.
However, by using the soft constraints for flow conservation, the proposed model can handle noisy
observations, where flow is not strictly preserved.
ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:14

T. Iwata et al.

(a) Tokyo(b) Osaka(c) Nagoya(d) Beijing
Fig. 6. Estimated cluster proportions q sk provided by the proposed method. White represents probability
one, and black represents zero. The cluster index is sorted by the start time point.

6

CONCLUSION

We have proposed a mixture of collective graphical models for estimating people flows fromspatio-temporal population data, and developed a variational Bayesian inference procedure for the
proposed model. We confirmed experimentally that our proposed method can estimate flow thatchanges over time, and predict the population at the next time point more precisely than existing
methods.
Although our results are encouraging, our framework can be further improved in a number ofways. Firstly, we plan to incorporate the spatial correlation of flows by using location dependent
ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

Collective Graphical Mixture Models

0:00–6:00

39:15

8:30–11:00

14:30–15:00

18:30–23:00

11:00–16:00

16:30–20:00

17:00–18:00

18:30–21:30

13:00–16:45

14:30–21:15(a) Tokyo

0:00–6:30

8:30–10:30(b) Osaka

0:00–7:00

9:00–16:30(c) Nagoya

0:00–7:45

8:30–12:45(d) Beijing

Fig. 7. Estimated people flows for each time-of-day cluster obtained with the proposed method using the (a)Tokyo, (b) Osaka, (c) Nagoya and (d) Beijing data. When an estimated transition probability is higher than athreshold, an arrow is drawn in that direction.
ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:16

T. Iwata et al.

mixture models as we did to incorporate time correlation in this paper. Secondly, we would liketo extend the proposed model to use information about day of week. Thirdly, we will study theeffectiveness of the proposed method with different grid sizes, time intervals and neighborhood
settings. Finally, it is important to develop a method for selecting neighbors that are appropriate
for the given data.
REFERENCES[1] Matthew J. Beal and Zoubin Ghahramani. 2003. The variational Bayesian EM algorithm for incomplete data: withapplication to scoring graphical model structures. Bayesian Statistics 7 (2003), 453–464.
[2] Vittoria Colizza, Alain Barrat, Marc Barthélemy, and Alessandro Vespignani. 2007. Predictability and epidemic
pathways in global outbreaks of infectious diseases: the SARS case study. BMC Medicine 5, 1 (2007), 34.
[3] Stephen Eubank, Hasan Guclu, VS Anil Kumar, Madhav V Marathe, Aravind Srinivasan, Zoltan Toroczkai, and Nan
Wang. 2004. Modelling disease outbreaks in realistic urban social networks. Nature 429, 6988 (2004), 180–184.
[4] Yong Ge, Hui Xiong, Alexander Tuzhilin, Keli Xiao, Marco Gruteser, and Michael Pazzani. 2010. An energy-efficientmobile recommender system. In Proceedings of the 16th ACM SIGKDD International Conference on Knowledge Discovery
and Data Mining. ACM, 899–908.
[5] Ronald L Hess, Ronald S Rubin, and Lawrence A West. 2004. Geographic information systems as a marketing
information system technology. Decision Support Systems 38, 2 (2004), 197–212.
[6] Akshat Kumar, Daniel Sheldon, and Biplav Srivastava. 2013. Collective diffusion over networks: Models and inference.
In Proceedings of International Conference on Uncertainty in Artificial Intelligence.
[7] Takeshi Kurashima, Tomoharu Iwata, Go Irie, and Ko Fujimura. 2010. Travel route recommendation using geotags inphoto sharing sites. In Proceedings of the 19th ACM International Conference on Information and Knowledge Management.
ACM, 579–588.
[8] Marco Lippi, Matteo Bertini, and Paolo Frasconi. 2010. Collective traffic forecasting. Machine Learning and Knowledge
Discovery in Databases (2010), 259–273.
[9] Dong C Liu and Jorge Nocedal. 1989. On the limited memory BFGS method for large scale optimization. Mathematical
Programming 45, 1-3 (1989), 503–528.
[10] Li-Ping Liu, Daniel Sheldon, and Thomas G Dietterich. 2014. Gaussian Approximation of Collective Graphical Models.
In Proceedings of International Conference on Machine Learning.
[11] Wei Liu, Yu Zheng, Sanjay Chawla, Jing Yuan, and Xie Xing. 2011. Discovering spatio-temporal causal interactions intraffic data streams. In Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data
Mining. ACM, 1010–1018.
[12] Thomas Minka. 2000. Estimating a Dirichlet distribution. (2000).
[13] A.B.M. Musa and Jakob Eriksson. 2012. Tracking unmodified smartphones using Wi-Fi monitors. In Proceedings of the10th ACM Conference on Embedded Network Sensor Systems. ACM, 281–294.
[14] Guande Qi, Xiaolong Li, Shijian Li, Gang Pan, Zonghui Wang, and Daqing Zhang. 2011. Measuring social functionsof city regions from large-scale taxi behaviors. In 2011 IEEE International Conference on Pervasive Computing andCommunications Workshops. IEEE, 384–388.
[15] Yoshihide Sekimoto, Ryosuke Shibasaki, Hiroshi Kanasugi, Tomotaka Usui, and Yasunobu Shimazaki. 2011. PFlow:
Reconstructing people flow recycling large-scale social survey data. IEEE Pervasive Computing 10, 4 (2011), 0027–35.
[16] D. Sheldon, M.A.S. Elmohamed, and D. Kozen. 2007. Collective inference on Markov models for modeling bird
migration. In Advances in Neural Information Processing Systems. 1321–1328.
[17] Daniel Sheldon, Tao Sun, Akshat Kumar, and Tom Dietterich. 2013. Approximate inference in collective graphicalmodels. In Proceedings of International Conference on Machine Learning.
[18] Daniel R Sheldon and Thomas G Dietterich. 2011. Collective graphical models. In Advances in Neural Information
Processing Systems. 1161–1169.
[19] Xuan Song, Quanshi Zhang, Yoshihide Sekimoto, Teerayut Horanont, Satoshi Ueyama, and Ryosuke Shibasaki. 2013.
Modeling and probabilistic reasoning of population evacuation during large-scale disaster. In Proceedings of the 19th
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 1231–1239.
[20] Xuan Song, Quanshi Zhang, Yoshihide Sekimoto, and Ryosuke Shibasaki. 2014. Prediction of human emergency
behavior and their mobility following large-scale disaster. In Proceedings of the 20th ACM SIGKDD International
Conference on Knowledge Discovery and Data Mining. ACM, 5–14.
[21] Tao Sun, Daniel Sheldon, Akshat Kumar, and EDU SG. 2015. Message Passing for Collective Graphical Models. In
Proceedings of International Conference on Machine Learning. 853–861.

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

Collective Graphical Mixture Models

39:17[22] Masayuki Terada, Tomohiro Nagata, and Motonari Kobayashi. 2012. Population estimation technology for mobile
spatial statistics. NTT DOCOMO Technical Journal 14, 3 (2012), 10–15.
[23] Apichon Witayangkurn, Teerayut Horanont, and Ryosuke Shibasaki. 2013. The Design of Large Scale Data Management
for Spatial Analysis on Mobile Phone Dataset. Asian Journal of Geoinformatics 13, 3 (2013).
[24] Jing Yuan, Yu Zheng, and Xing Xie. 2012. Discovering regions of different functions in a city using human mobility
and POIs. In Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining.
ACM, 186–194.
[25] Jing Yuan, Yu Zheng, Xing Xie, and Guangzhong Sun. 2011. Driving with knowledge from the physical world. In
Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 316–324.
[26] Jing Yuan, Yu Zheng, Chengyang Zhang, Wenlei Xie, Xing Xie, Guangzhong Sun, and Yan Huang. 2010. T-drive:
driving directions based on taxi trajectories. In Proceedings of the 18th SIGSPATIAL International Conference on Advancesin Geographic Information Systems. ACM, 99–108.
[27] Yu Zheng, Yanchi Liu, Jing Yuan, and Xing Xie. 2011. Urban computing with taxicabs. In Proceedings of the 13th
International Conference on Ubiquitous Computing. ACM, 89–98.

A

VARIATIONAL LOWER BOUND

In this appendix, we give the the lower bound of the log marginal likelihood (9). The first term isgiven byE[log p(M|N, Θ, z, s)]
T
−1 X
I 
K

X
X
X
X
0
0
=
log N ti ! −log M ti j ! +M ti jq s (t )k [Ψ(α ki
)
−
Ψ(
ᾱ
)]
,
j
ki
t=1 i=1

j ∈J ij ∈J i(30)

k=1where we used the following expected value of the log of a multinomial parameter under a Dirichlet
distribution,000 ) [log θ ki j ] = Ψ(αEq(α ki
ki j ) − Ψ(ᾱ ki ).

(31)

By using the above equation, the second, third and fourth terms become
E[log p(Θ|α )] =

K X
I 
X

log Γ(

k=1 i=1+

Xj ∈J iX

α r (i, j) ) −j ∈J iX

log Γ(α r (i, j) )j ∈J i00(α r (i, j) − 1)[Ψ(α kij ) − Ψ(ᾱ ki )] ,E[log p(z|ϕ)] =S X
KX
s=1 k=1

q sk [Ψ(β k0 ) − Ψ( β¯0 )],E[log p(ϕ|β )] = log Γ(βK ) − K log Γ(β ) + (β − 1)KX
[Ψ(β k0 ) − Ψ( β¯0 )].

(32)

(33)

(34)

k=1

The fifth, sixth and seventh terms are given by
E[log p(g|z,τ ,η)] =

K
0 S
 T
1 X
Tk a Xq sk (дs − f k0 ) 2 − log 2π ,Tk [Ψ(a k0 ) − log b k0 ] − 0 − k02
d k b k s=1
2

(35)

k=1

E[log p(τ | f ,d,η)] =

K

da 0
1 X
d
d
Ψ(a k0 ) − log b k0 + log− 0 − 0k ( f k0 − f ) 2 ,2
2π d k
bk

(36)

k=1

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

39:18T. Iwata et al.
E[log p(η|a,b)] = −K log Γ(a) + Ka log b + (a − 1)

K
X
k=1

K
X
a0
b0
a k0
(Ψ( k ) − log k ) − b
,
2
2
b k0

(37)

k=1

where we used the following expected value of a positive parameter under a Gamma distribution,
Eq(η k ) [log η k ] = Ψ(a k0 ) − log b k0 .
The eighth term is the entropy of the Dirichlet distribution, and is given by,
K X
I X

X
X
0
0
0
¯0 ki ) −¯0 ki )] .
H[q(Θ)] =log Γ(α ki
)
−
log
Γ(
α
(α
−
1)[Ψ(α
)
−
Ψ(
α
j
ki j
ki j
k=1 i=1 j ∈J i

(38)

(39)

j ∈J iThe ninth term is the entropy of the discrete distribution, and is
H[q(z)] = −SXq sk log q sk .

(40)

s=1

The tenth term is the entropy of the Dirichlet distribution, and is
H[q(ϕ)] =KX
k=1

Γ(β k0 ) − log Γ( β¯0 ) −

K
X
k=1

(β k0 − 1)[Ψ(β k0 ) − Ψ( β¯0 )].

(41)

The eleventh term is the entropy of the Gaussian distribution, and isH[q(τ )] =

K

K

k=1

k=1

1X
K
1X
(log 2π + 1) −log d k0 −(Ψ(a k0 ) − log b k0 ).
2
2
2

(42)

The twelfth term is the entropy of the Gamma distribution, and isH[q(η)] =KX[log Γ(a k0 ) − (a k0 − 1)Ψ(a k0 ) − log b k0 + a k0 ].
k=1

Received February 2007; revised March 2009; accepted June 2009

ACM Transactions on Spatial Algorithms and Systems, Vol. 9, No. 4, Article 39. Publication date: March 2017.

(43)