ORIGINAL RESEARCH
published: 20 November 2019
doi: 10.3389/frobt.2019.00115

HVGH: Unsupervised Segmentationfor High-Dimensional Time Series
Using Deep Neural Compression andStatistical Generative Model
Masatoshi Nagano 1*, Tomoaki Nakamura 1 , Takayuki Nagai 2,3 , Daichi Mochihashi 4 ,Ichiro Kobayashi 5 and Wataru Takano 6
1

Department of Mechanical Engineering and Intelligent Systems, The University of Electro-Communications, Tokyo, Japan,
Graduate School of Engineering Science, Osaka University, Osaka, Japan, 3 Artificial Intelligence Exploration ResearchCenter, The University of Electro-Communications, Tokyo, Japan, 4 Department of Statistical Inference and Mathematics, The
Institute of Statistical Mathematics, Tokyo, Japan, 5 Advanced Sciences, Graduate School of Humanities and Sciences,
Ochanomizu University, Tokyo, Japan, 6 Center for Mathematical Modeling and Data Science, Osaka University, Osaka, Japan

2

Edited by:
Georg Martius,
Max Planck Institute for Intelligent
Systems, Germany
Reviewed by:
Arash Mehrjou,
Max Planck Institute for Intelligent
Systems, Germany
Dominik M. Endres,
University of Marburg, Germany*Correspondence:
Masatoshi Naganon1832072@edu.cc.uec.ac.jp
Specialty section:
This article was submitted toComputational Intelligence in
Robotics,
a section of the journal
Frontiers in Robotics and AI
Received: 31 March 2019
Accepted: 22 October 2019
Published: 20 November 2019
Citation:
Nagano M, Nakamura T, Nagai T,
Mochihashi D, Kobayashi I and
Takano W (2019) HVGH:
Unsupervised Segmentation for
High-Dimensional Time Series Using
Deep Neural Compression andStatistical Generative Model.
Front. Robot. AI 6:115.
doi: 10.3389/frobt.2019.00115

Humans perceive continuous high-dimensional information by dividing it into meaningful
segments, such as words and units of motion. We believe that such unsupervised
segmentation is also important for robots to learn topics such as language andmotion. To this end, we previously proposed a hierarchical Dirichlet process–Gaussian
process–hidden semi-Markov model (HDP-GP-HSMM). However, an important
drawback of this model is that it cannot divide high-dimensional time-series data.
Furthermore, low-dimensional features must be extracted in advance. Segmentation
largely depends on the design of features, and it is difficult to design effective
features, especially in the case of high-dimensional data. To overcome this problem,
this study proposes a hierarchical Dirichlet process–variational autoencoder–Gaussian
process–hidden semi-Markov model (HVGH). The parameters of the proposed HVGH
are estimated through a mutual learning loop of the variational autoencoder andour previously proposed HDP-GP-HSMM. Hence, HVGH can extract features fromhigh-dimensional time-series data while simultaneously dividing it into segments inan unsupervised manner. In an experiment, we used various motion-capture data todemonstrate that our proposed model estimates the correct number of classes andmore accurate segments than baseline methods. Moreover, we show that the proposed
method can learn latent space suitable for segmentation.
Keywords: motion segmentation, Gaussian process, variational autoencoder, hidden semi-Markov model, motion
capture data, high-dimensional time-series data

1. INTRODUCTION
Humans perceive continuous high-dimensional information by dividing it into meaningful
segments, such as words and units of motion. For example, we recognize words by segmenting
speech waves, and we perceive particular motions by segmenting continuous motion. Humans
learn words and motions by appropriately segmenting continuous information without explicit
segmentation points. We believe that such unsupervised segmentation is also important for robots,in order for them to learn language and motion.

Frontiers in Robotics and AI | www.frontiersin.org

1

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

Furthermore, studies have proposed methods of detectingchange points in time-series data in an unsupervised manner(Yamanishi et al., 2002; Lund et al., 2007; Liu et al., 2013; Haber
et al., 2014). These are the methods of finding points withdifferent fluctuations based on previous observations. Therefore,
these methods assume that similar temporal patterns are repeatedbetween the change points. On the other hand, in this study, we
consider that the segments comprise not only repeated patterns
but also an arbitrary pattern. Thus, the change points do notnecessarily indicate the segment boundaries.
In some studies, segmentation is formulated stochasticallyusing hidden Markov models (HMMs) (Beal et al., 2002; Fox
et al., 2011; Taniguchi et al., 2011; Matsubara et al., 2014).
However, it is difficult for HMMs to represent complicatedmotion patterns. Instead, we use Gaussian processes inour model, which is a type of non-parametric model thatcan better represent complicated time-series data compared
to HMMs. Figure 2B shows how an HMM represents the
trajectory of data points shown in Figure 2A. The HMM
represents the trajectory using five Gaussian distributions.
However, one can see that the details of the trajectory are
lost. On the other hand, in HDP-GP-HSMM (Figure 2C), thetrajectory can be represented continuously using two Gaussian
processes (GPs). We confirmed that our GP-based model canestimate segments more accurately than HMM-based methods(Nagano et al., 2018).
In some studies, the number of classes is estimated byintroducing a hierarchical Dirichlet process (HDP) into an HMM(Beal et al., 2002; Fox et al., 2007). An HDP is a method of
estimating the number of classes by assuming an infinite numberof classes. Fox et al. extended an HDP–HMM to develop a socalled sticky HDP-HMM, which prevents over-segmentation byincreasing the self-transition probability (Fox et al., 2007).
Among methods of combining statistical models with neuralnetworks, a method of classifying complicated data using a GMM
and VAE was proposed (Johnson et al., 2016). In contrast, ourproposed HVGH is a model that combines a statistical modelwith a VAE for segmenting high-dimensional time-series data.

In this paper, we define the segments as arbitrary temporal
patterns that appear multiple times in the time-series data.
Additionally, we have proposed a method to extract the segments
by capturing such a nature stochastically. One of our previous
methods is the hierarchical Dirichlet process–Gaussian process–hidden semi-Markov model (HDP-GP-HSMM) (Nagano et al.,2018). HDP-GP-HSMM is a non-parametric Bayesian model thatis a hidden semi-Markov model, the emission distributions ofwhich are Gaussian processes (MacKay, 1998), and it facilitates
the segmentation of time-series data in an unsupervised manner.
In this model, segments are continuously represented using a
Gaussian process. Moreover, the number of segmented classes
can be estimated using hierarchical Dirichlet processes (Teh et al.,2006). The Dirichlet processes assume an infinite number ofclasses. However, only a finite number of classes are actually used.
This is accomplished by stochastically truncating the number ofclasses using a slice sampler (Van Gael et al., 2008).
However, our HDP-GP-HSMM cannot handle highdimensional data, and feature extraction is required to reducethe dimensionality in advance. Indeed, segmentation largely
depends on this feature extraction, and it is difficult to extracteffective features, especially in the case of high-dimensional data.
To overcome this problem, this study introduces a variational
autoencoder (VAE) (Kingma et al., 2013) to HDP-GP-HSMM.
Thus, the model we propose in this paper is a hierarchical
Dirichlet process–variational autoencoder–Gaussian process–hidden semi-Markov model (HVGH1 ). Figure 1 shows an
overview of HVGH. The observation sequence is compressedand converted into a latent variable sequence by the VAE,
and the latent variable sequence is divided into segments by
HDP-GP-HSMM. Furthermore, parameters learned by HDPGP-HSMM are used as the hyperparameters for the VAE. In thisway, the parameters are optimized in a mutual learning loop, andappropriate latent space for segmentation can be learned by the
VAE. In experiments, we evaluated the efficiency of the proposed
HVGH using real motion-capture datasets. The experimental
results show that HVGH achieves a higher accuracy compared tobaseline methods.
Many studies on unsupervised motion segmentation havebeen conducted. However, heuristic assumptions are used inmany of them (Lioutikov et al., 2015; Wächter et al., 2015; Takano
et al., 2016). Wächter et al. proposed a method for segmenting
object-manipulation motion in robots and used contact between
the end-effector and the object as a segmentation clue (Wächteret al., 2015). Lioutikov et al. proposed a method that requires
candidates for the segmentation points in advance (Lioutikov
et al., 2015). In addition, the method proposed by Takano
et al. used errors between the predicted and actual values ascriteria for segmentation (Takano et al., 2016). Moreover, somemethods use motion features such as the zero velocity of joint
angles (Fod et al., 2002; Shiratori et al., 2004; Lin et al., 2012).
However, this assumption typically induces over-segmentation(Lioutikov et al., 2015).

2. HIERARCHICAL DIRICHLETPROCESS–VARIATIONAL
AUTOENCODER–GAUSSIAN
PROCESS–HIDDEN SEMI-MARKOV
MODEL (HVGH)
Figure 3 shows a graphical model of our proposed HVGH,
which is a generative model of time-series data. In this model,
cj (j = 1, 2, · · · , ∞) denotes the classes of the segments, where
the number of classes is assumed to be countably infinite.
β denotes an infinite-dimensional multinomial distribution,
which is generated from the GEM distribution (Pitman, 2002),
parameterized by γ . GEM denotes the co-authors Griffiths,
Engen, and McCloskey—with the so-called stick-breaking
process (SBP) (Sethuraman, 1994). The probability π c denotes
the transition probability, which is generated by the Dirichlet

1 HVGHstands for hierarchical Dirichlet process, variational autoencoder,
Gaussian process, and hidden semi-Markov model.

Frontiers in Robotics and AI | www.frontiersin.org

2

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

FIGURE 1 | Overview of the generative process of the HVGH.

FIGURE 2 | Representing a trajectory: (A) Observed data points, (B) representation by HMM, and (C) representation by HDP-GP-HSMM.

1998) with the parameter φc as follows:

process (Teh et al., 2006), parameterized by β:β ∼ GEM(γ ),π c ∼ DP(η, β).

(1)cj ∼ P(c|cj−1 , π c ),

(3)

(2)

Zj ∼ G P(Z|φcj ),(4)where φc represents the parameter of the Gaussian process
corresponding to the class c and is a set of segments classified
into the class c in the learning phase. The segment X j is generated
from the segment of the latent variables Zj by using the decoder
pdec of the VAE:

γ and η are the concentration parameters of the Dirichlet
processes; the smaller the value of the concentration parameter,
the sparser the generated distribution is. The process bywhich the probability distribution is constructed through a
two-phase Dirichlet process is called a hierarchical Dirichlet
process (HDP) (Teh et al., 2006). HDP is described in detail in
Nagano et al. (2018).
The class cj of the j-th segment is determined by the class of the(j − 1)-th segment and transition probability π c . The segment oflatent variables Zj is generated by a Gaussian process (MacKay,

Frontiers in Robotics and AI | www.frontiersin.org

X j ∼ pdec (X|Zj ).

(5)

The observed sequence s = X 1 , X 2 , · · · , X J is assumed to begenerated by connecting segments X j sampled by the above
generative process. Similarly, the sequence of the latent variables

3

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

distribution with the parameters µ and σ 2 :s̄ = Z1 , Z2 , · · · , ZJ is obtained by connecting the segments of thelatent variables Zj . In this paper, the i-th data point included in
X j is described as xji , and the i-th data point included in Zj isdescribed as z ji . If the characters represent what they obviously
are, we omit their subscripts.
The generative process of the observed sequence s describedabove is summarized in Algorithm 1. This pseudo code
represents the generative process, and it is difficult to directly
implement this code because the number of classes is infinite. The
details on how to address this problem are described in section 3.

p(ẑ|î, φc ) ∝ N (z|µ, σ 2 ),(6)

µ = kT C−1 t c ,(7)

σ 2 = ρ − kT C−1 k.

(8)Here, C is a matrix having the following elements:
C(ip , iq ) = k(ip , iq ) + ω−1 δpq ,(9)where k(·, ·) denotes the kernel function and ω denotes a
hyperparameter that represents noise in the observations. k is avector with the elements k(ip , î), and ρ is k(î, î). φc represents a setof the segments of latent variables that are classified into the classc, and t c and ic are the vectors where a respective latent variable
zi and time step i in φc are arranged. For example, assuming thatthe latent variables have one dimension, and segments Z1 , Z2 , · · ·are classified into class c, we can compute t c and ic as follows:

2.1. Gaussian Process (GP)In this paper, each class represents a continuous trajectory bylearning the emission zi of time step i using a Gaussian process(MacKay, 1998). In the Gaussian process, given t c and ic , which
are the vectors of the latent variable zi that is classified intoclass c and its time step i, respectively (details are explainedlater), the predictive distribution of ẑ of time step î is a Gaussian

φc = {Z1 , Z2 , · · · } = {{z11 , z12 , z13 , · · · }, {z21 , z22 , z23 , · · · }, · · · },(10)t c = [z11 , z12 , z13 , · · · , z21 , z22 , z23 , · · · ]T ,T

ic = [1, 2, 3, · · · 1, 2, 3, · · · ] .

(11)
(12)

Once t c and ic are computed, they are shared to compute theprobability of (ẑ, î) being classified into class c. A Gaussian
process can represent complicated time-series data owing to thekernel function. In this paper, we used the following kernel
function:1k(ip , iq ) = θ0 exp(− θ1 ||ip − iq ||2 ) + θ2 + θ3 ip iq ,(13)2where θ∗ denotes the parameters of the kernel, which are fixed
for all classes in our experiments. The reason why we selectthis kernel is that the motions are generally smooth and, hence,
we consider the latent variable to also be temporally smooth.
Figure 4 illustrates the samples from various kernels: linear,
exponential, periodic, and radial basic function (RBF) kernels(Bishop, 2006). As can be observed in this figure, it is difficult
for the linear and periodic kernels to represent a non-linear andnon-periodic pattern, and the sample of the exponential kernel is
not smooth. On the other hand, the RBF kernel can represent asmooth temporal pattern. Therefore, we use the kernel based on
the RBF kernel, which is generally used for Gaussian processes(Bishop, 2006). However, it is not evident as to which kernel isthe most appropriate. Moreover, the appropriate kernel depends
on the task. This issue will be considered in future work because
it exceeds the scope of this paper.
Additionally, if the observations are composed ofmultidimensional vectors, we assume that each dimension is
independently generated. Therefore, the predictive distribution
G P(z|φc ) by which the emission z = (z(1) , z(2) , · · · ) of time stepi is generated using a Gaussian process of class c is computed asfollows:

FIGURE 3 | Graphical model of HVGH. The white nodes represent
unobserved variables, and the gray node represents the high-dimensional
observed sequence that is obtained by concatenating segments.

Algorithm 1: Generative process of s by HVGH
1: Draw β ∼ GEM(γ )2: for c = 1, · · · , ∞ do
3:
Draw π c ∼ DP(η, β)
4: end for
5:
6: s = ǫ (an empty sequence)
7: for j = 1, · · · , J do

Draw cj ∼ P(cj |cj−1 , π c ),Draw Zj ∼ G P(Zj |φcj ),10:Draw X j ∼ pdec (X j |Zj ),11:
Append X j to s.
12: end for
8:9:

Frontiers in Robotics and AI | www.frontiersin.org

GP (z|φc ) = p(z(1) |i, φc,1 )p(z(2) |i, φc,2 )p(z(3) |i, φc,3 ) · · ·= N (z(1)(14)|µ1 , σ12 )N (z(2) |µ2 , σ22 )N (z(3) |µ3 , σ32 ) · · ·(15)

4November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

FIGURE 4 | Four temporal patterns sampled from the Gaussian processes, where (A) linear : ip iq , (B) exponential : exp(−|ip − iq |), (C) periodic : exp(cos |ip − iq |), and(D) RBF : exp(−|ip − iq |2 ) are used.

By using this probability, the latent variables can be classified into
the classes. Moreover, because each dimension is independently
generated, the mean vector µc (i) and the variance–covariancematrix 6c (i) of GP (z ji |φc ) are represented as follows:µc (i) = (µ1 , µ2 , µ3 , · · · ), 2
σ1 0 02
6c (i) =  0 σ2 0  ,
..
0 0 .

(16)
(17)
FIGURE 5 | Overview of the variational autoencoder (VAE). The
low-dimensional latent variable z j is obtained by compressing the observed
data point xj through the encoder network. x′j is an observation reconstructedby the decoder from the latent variable z j .

where (µ1 , µ2 , µ3 , · · · ) and (σ12 , σ22 , σ32 , · · · ), respectively,
represent the mean and variance of each dimension. HVGH is amodel in which the VAE and GP influence each other mutuallywith the use of µc (i) and 6c (i) as the prior distribution of
the VAE.

6enc (x) that are estimated by using encoder networks from input
x is used as qenc (z):2.2. Overview of the Variational
Autoencoder
In this paper, we compress a high-dimensional time-series
observation into low-dimensional latent variables using the VAE(Kingma et al., 2013). The VAE is a neural network that can learn
the correspondence between a high-dimensional observation x
and the latent variable z. Generally, in a probabilistic model, theposterior distribution of z can be expressed as follows:
pdec (x|z)p(z)p(z|x) =.
p(x)qenc (z) = N (z|µenc (x), 6enc (x)).

The latent variable z is stochastically determined by this
distribution, and x′ is generated through decoder networks pdec :z ∼ qenc (z),(20)
x′ ∼ pdec (x|z).

(18)

(21)

The parameters of the encoder and decoder are determined
to maximize the likelihood by using the variational Bayesian
method. Generally, the prior distribution of the VAE is a
Gaussian distribution, the mean of which is the zero vector
0, and the variance–covariance matrix is the identity matrix e.
On the other hand, HVGH uses a Gaussian distribution whoseparameters are µc (i) and 6c (i) of class c into which z ji is

However, if a neural network that has expressive power is used
for the generative model pdec (x|z), Equation (18) cannot beanalytically derived. To solve this problem, in the VAE, p(z|x)
is approximated by qenc (z). Figure 5 shows an overview of the
VAE. A Gaussian distribution with a mean µenc (x) and variance

Frontiers in Robotics and AI | www.frontiersin.org

(19)

5

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

3.1. Parameter Inference of
HDP-GP-HSMM

classified. As a result, latent space suitable for segmentation canbe constructed. By using this VAE, a sequence of the observations = X 1 , X 2 , · · · , X J is converted into a sequence of the latentvariables s̄ = Z1 , Z2 , · · · , ZJ through the encoder.

3.1.1. Blocked Gibbs Sampler
In HDP-GP-HSMM, segments and classes of latent variables are
determined by sampling. For efficient estimation, we utilize a
blocked Gibbs sampler (Jensen et al., 1995), in which all segments
and their classes in one observed sequence are sampled. First,
all sequences of the latent variables are randomly divided into
segments and randomly classified into classes. Next, the segmentsof latent variables Znj (j = 1, 2, · · · , Jn ) obtained by segmenting
the n-th sequence s̄n are excluded from the training data,
and the parameters of the Gaussian process φc and transition
probabilities P(c|c′ ) are updated. The segments of latent variables
and their classes are sampled as follows:

3. PARAMETER INFERENCE
The log likelihood of HVGH is expressed as follows:
log p(X 1 , · · · , X J , c1 , · · · , cJ )YZ
= logp(Zj , cj )p(X j |Zj )dZj
Zj
j
YZ
= log
GP (Z j |φc )P(cj |cj−1 , π c ) p(X j |Zj ) dZ j .
Zj |
j
{z
} | {z }HDP−GP−HSMM

(22)
(23)
(24)

(Zn,1 , · · · , Zn,Jn ), (cn,1 , · · · , cn,Jn ) ∼ p((Z1 , Z2 , · · · , ZJ ),(c1 , c2 , · · · , cJ )|s̄n ).

VAE

(25)

The parameters of the Gaussian process of each class φc andtransition probability P(c|c′ ) are updated by using the sampled
segments and their classes. The parameters are optimized by
iterating this procedure. Algorithm 2 is the pseudo code of theblocked Gibbs sampler. Ncnj and Ncnj ,cn,j+1 are parameters tocompute the transition probability in Equation (29). However, it
is difficult to compute Equation (25) because an infinite numberof classes is assumed. To overcome this problem, we use aslice sampler to compute these probabilities by stochasticallytruncating the number of classes.
Moreover, the probabilities of all possible patterns of segments
and classifications are required in Equation (25), and thesecannot be computed naively owing to the large computationalcost. To compute Equation (25) efficiently, we utilize forward
filtering–backward sampling (Uchiumi et al., 2015).

In Equation (24), the first and second factors are computedin HDP-GP-HSMM, and the third factor is computed in VAE.
It is difficult to directly maximize Equation (22); therefore,
HDP-GP-HSMM and the VAE are alternately optimized, andthe parameters that approximately maximize Equation (22)are computed. Figure 6 depicts an outline of the parameter
estimation for HVGH. A sequence of observations s =X 1 , X 2 , · · · , X J is converted into a sequence of latent variabless̄ = Z1 , Z2 , · · · , ZJ by the VAE. Then, through HDP-GP-HSMM,
the sequence of latent variables s̄ is divided into segments oflatent variables Z1 , Z2 , · · · , ZJ , and the parameters µc (i) and6c (i) of the predictive distribution of z are computed. Thispredictive distribution is used as a prior distribution of the
VAE. Thus, the parameters of the VAE and HDP-GP-HSMM are
mutually optimized.

3.1.2. Slice Sampler
In the HDP, we assumed that the number of classes is countably
infinite. However, it is difficult to compute Equation (25) because
c can have infinite possibilities. To overcome this problem, we use

Algorithm 2: Blocked Gibbs Sampler
1: Repeat until convergence
2: for n = 1 to N do
3:for j = 1 to Jn do
4:
Ncnj − = 1
5:
Ncnj ,cn,j+1 − = 1
6:Delete segments Znj from φcnj
7:
end for8:// Sampling segments and their classes
9:Zn∗ , cn∗ ∼ p(Z∗ , c∗ |sn )10:
for j = 1 to Jn do
11:
Ncnj + +12:
Ncnj ,cn,j+1 + +13:Append segments Znj to φcnj
14:
end for
15: end for

FIGURE 6 | Overview of parameter estimation for HVGH. The parameters arelearned by a mutual learning loop of the VAE and HDP-GP-HSMM.

Frontiers in Robotics and AI | www.frontiersin.org

6

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGHa slice sampler (Van Gael et al., 2008) to stochastically truncate
the number of classes. In the slice sampler, an auxiliary variable
uj that follows the distribution for each time step j is introduced:p(uj |cj−1 , cj ) =ξ (0 < uj < πcj−1 ,cj )πcj−1 ,cj

segments. Plen (k|λ) represents a Poisson distribution with a mean
parameter λ:Plen (k|λ) =

,

(26)

λk e−λ.
k!

(28)This corresponds to the distribution of the segment lengths. In
addition, P(c|c′ , β, η) is the transition probability, which can becomputed as follows:

where ξ (A) = 1 if condition A is true; otherwise, it is 0. By
truncating the classes with a transition probability πcj−1 ,cj thatis less than uj , the number of classes becomes finite, as shownin Figure 7.

P(c|c′ , β, η) =Nc′ c + ηβc′,Nc′ + η(29)

3.1.3. Forward Filtering–Backward Samplingwhere Nc′ represents the number of segments of class c andNc′ c denotes the number of transitions from c′ to c. In addition,k′ and c′ are the length and class of possible segments befores̄t−k : k , respectively, and these probabilities are marginalizedout in Equation (27). Moreover, α[t][k][∗] = 0 if t − k <0, and α[0][0][∗] = 1.0. Equation (27) can be recursively
computed from α[1][1][∗] using dynamic programming, asshown in Figure 8A. This figure depicts an example of computing
α[t][2][2], which is the probability that the two samples beforet become a segment having the class c. In this case, the lengthis two and, therefore, all the segments with end points t −2 can potentially transit to this segment. In α[t][2][2], these
possibilities are marginalized out. Finally, the length of segments
and their classes can be sampled through backward samplingfrom t = T:

The number of classes can be truncated by slice sampling.
Consequently, forward filtering–backward sampling (Uchiumi
et al., 2015) can be applied to compute Equation (25). In forward
filtering, the probability that k samples s̄t−k : k before time step tform a segment of class c is as follows:
α[t][k][c]= G P(s̄t−k : k |φc )Plen (k|λ)

C̄
K X
X

P(c|c′ , β, η)α[t − k][k′ ][c′ ] ,k′ =1 c′ =1(27)
where C̄ denotes the maximum number of classes estimatedby slice sampling and K denotes the maximum length ofk, c ∼ p(z j |st−k : t )α[t][k][c]P(cj−1 |c).

(30)

Figure 8B depicts an example of backward sampling. The lengthof segment k1 and its class c1 are sampled from the probabilities
α[t][∗][∗]. If k1 = 2, k2 and c2 are sampled from α[t − 2][∗][∗].
By iterating this procedure until t = 0, the segments and theirclasses can be determined. Algorithm 3 shows the pseudo-codeof forward filtering–backward sampling with slice sampling.

3.2. Parameter Inference of the VAE

FIGURE 7 | Slice sampling truncates the number of classes by thresholding
πc′ ,c .

The parameters of the encoder and decoder of VAE are estimatedto maximize the likelihood p(x). However, it is difficult to

FIGURE 8 | (A,B) Forward filtering–backward sampling.

Frontiers in Robotics and AI | www.frontiersin.org

7

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGHmaximize the likelihood directly. Instead, the normal VAE
maximizes the following variational lower limit:
Z
L(xji , z ji ) = qenc (z ji |xji ) log pdec (xji |z ji )dz ji−DKL (qenc (z ji |xji )||p(z ji |0, e)),¯

the variance–covariance matrix is e. DKL (qenc (z ji |xji )||p(z ji |0, e))¯is the Kullback–Leibler divergence, and this functions as aregularization term. On the other hand, in HVGH, the meanµc (i) and the variance–covariance matrix 6c (i) are used as theparameters of the prior distribution. These are the parameters of
the predictive distribution of class c into which z ji is classified,
and they are estimated by HDP-GP-HSMM:
Z
L(xji , z ji ) = qenc (z ji |xji ) log pdec (xji |z ji )dz ji

(31)

R
where
qenc (z ji |xji ) log pdec (xji |z ji )dz ji
represents
thereconstruction error. Moreover, p(z ji |0, e) is a prior distribution¯
of z ji , which is a Gaussian distribution whose mean is 0, and−DKL (qenc (z ji |xji )||p(z ji |µc (i), 6c (i))).

(32)

Figure 9 illustrates the difference in prior distributions between
Equations (31, 32). In the normal VAE using Equation (31), theprior distribution is N (0, e) against all data points, as shown¯
in Figure 9A. On the other hand, the parameters of the priordistribution of HVGH are computed by the Gaussian process,as shown in Figure 9B. Because the GP restricts data points thathave closer time steps to being more similar values, zji becomesa similar value to zj,i−1 and zj,i+1 . Therefore, the latent spacelearned by the VAE can reflect the characteristics of time-series
data. Moreover, these parameters have different values dependingon the class of the data point. Therefore, the latent space can also
reflect the characteristics of each class.

Algorithm 3: Forward Filtering–Backward Sampling
1: // Slice sampling
2: for j = 1 to Jn do
3:uj ∼ p(uj |cj−1 , cj )4: end for
5: C̄ = maxj (count(πcj−1 ,cj > uj ))6: // Forward filtering
7: for t = 1 to T do
8:for k = 1 to K do
9:for c = 1 to C̄ do
10:Compute α[t][k][c]11:
end for12:
end for
13: end for
14: // Backward sampling
15: t = T, j = 0, c0 = 0
16: while t > 0 do
17:
j=j+1
18:
k, c ∼ p(z j |st−k : t )α[t][k][c]P(cj−1 |c)19:
z j = st−k : t
20:
cj = c21:t =t−k22: end while
23: Jn = j
24: return (z Jn , z Jn −1 , · · · , z 1 ), (cJn , cJn −1 , · · · , c1 )4. EXPERIMENTS
To validate the proposed HVGH, we applied it to several typesof time-series data. For comparison, we used HDP-GP-HSMM(Nagano et al., 2018), HDP-HMM (Beal et al., 2002), HDPHMM+NPYLM (Taniguchi et al., 2011), BP-HMM (Fox et al.,2011), and Autoplait (Matsubara et al., 2014) as baseline methods.

4.1. Experimental Setup
To evaluate the validity of the proposed method, we used thefollowing four motion-capture datasets.
• Chicken dance: We used a sequence of motion-capture dataof a human performing a chicken dance from the CMUFIGURE 9 | Influence of prior distribution. The orange region represents the standard deviation. (A) The same prior distribution is used for any data points in thenormal VAE. (B) The distribution is varied depending on the time step and the class of the data point in the VAE in HVGH.

Frontiers in Robotics and AI | www.frontiersin.org

8

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

FIGURE 10 | Four unit motions included in the chicken dance: (A) beaks, (B)wings, (C) tail feathers, and (D) claps.
FIGURE 12 | Seven unit motions included in exercise motion 1: (A) jumping
jack, (B) twist, (C) arm circle, (D) bend over, (E) knee raise, (F) squatting, and(G) jogging.

λ = 14.0, θ0 = 1.0, θ1 = 1.0, θ2 = 0.0, and θ3 = 16.0, whichwere empirically determined for the segmentation of the 4-fps
sequences. For the chicken dance exclusively, we set λ to half thatof the others because its motion-capture data was shorter thanthe others. To train the VAE, we used 1/4 of all the data points
as a mini batch, Adam (Kingma et al., 2017) was used for the
optimization, and the optimization was iterated 150 times. To
train HDP-GP-HSMM, the blocked Gibbs sampler was iterated
10 times to converge the parameters. Furthermore, the mutual
learning loop of the VAE and HDP-GP-HSMM was iterated untilthe variational lower limit converged.

FIGURE 11 | Seven unit motions included in the “I’m a little teapot” dance: (A)short and stout, (B) bending knee, (C) spread arms, (D) handle, (E) spout, (F)steam up, and (G) pour.

4.2. Evaluation Metrics
To evaluate the segmentation accuracy, we used fourmeasures: the normalized Hamming distance, precision,
recall, and F-measure.
The normalized Hamming distance represents the error rateof the classification of the data points, and it is computed
as follows:

Graphics Lab Motion Capture Database2 . The dance includes
four motions, as shown in Figure 10.
• “I’m a little teapot” dance (teapot dance): We also used
two sequences from the teapot dance motion-capture datafrom subject 29 in the CMU Graphics Lab Motion Capture
Database3 . These sequences include seven motions, as shownin Figure 11.
• Exercise motion 1: To determine the validity against more
complicated motions, we used three sequences of exercise
motion-capture data from subject 13 in the CMU Graphics
Lab Motion Capture Database. These sequences include sevenmotions, as shown in Figure 12.
• Exercise motion 2: Furthermore, we used three sequences ofdifferent exercises from the motion-capture data from subject
14 in the CMU Graphics Lab Motion Capture Database. These
sequences include 11 motions, as shown in Figure 13.

ND(c, c) =(33)

where c and c, respectively, represent sequences of estimatedclasses and correct classes in the data points in the observed
sequence. Moreover, D(c, c) represents the Hamming distance
between two sequences, and |c| is the length of the sequence.
Therefore, the normalized Hamming distance ranges from zero
to one, and smaller values indicate that the estimated classes are
more similar to the ground truth.
To compute the precision, recall, and F-measure, we evaluatedboundary points (boundaries between segments) as true positives(TPs), true negatives (TNs), false positives (FPs), and false
negatives (FNs), as shown in Figure 14. A TP is assigned tothe points that are correctly estimated as boundary points. An
estimated boundary point is treated as correct if the estimatedboundary is within the error range, as shown in Figure 14, Frame(2). The error range is defined as ±ψ% of the sequence length,and ψ represents the percentage of the error range. A TN is
assigned to the points that are correctly estimated not to beboundary points, as shown in Figure 14, Frame (3). Conversely,
FPs and FNs are assigned to points that are falsely estimated asboundary points, as shown in Figure 14, Frame (10), and falsely

To reduce computational cost, all the sequences werepreprocessed by down sampling to 4 fps. These motioncapture datasets included the directions of 31 body parts, eachof which was represented by a three-dimensional Euler angle.
Therefore, each frame was constructed in 93-dimensional
vectors. We used sequences of 93-dimensional vectors as input.
Moreover, HVGH requires hyperparameters, and we set them to
2 http://mocap.cs.cmu.edu/: subject 18, trial 15.
3 http://mocap.cs.cmu.edu/: subject 29, trials 3 and 8.

Frontiers in Robotics and AI | www.frontiersin.org

D(c, c),|c|

9

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

FIGURE 13 | Eleven unit motions included in exercise motion 2: (A) jumping jack, (B) jogging, (C) squatting, (D) knee raise, (E) arm circle, (F) twist, (G) side reach, (H)boxing, (I) arm wave, (J) side bend, and (K) toe touch.

FIGURE 14 | Example of segmentation evaluation: TP is assigned to the boundary (2) because the estimated boundary is within the error range from the trueboundary.

estimated not to be boundary points, as shown in Figure 14,
Frame (6), respectively. From these boundary evaluations, theprecision, recall, and F-measure are computed as follows:
NTP,
NTP + NFPNTP,R =NTP + NFN
2P · R
F =
,
P+R
P =

(34)
(35)
(36)

where NTP , NFP , and NFN represent the number of boundary
points estimated as TP, FP, and, FN, respectively.
Figure 15 depicts the results of a preliminary experiment todetermine the error range. The horizontal axis represents thepercentages of error range ψ, and the vertical axis represents
the average F-measures of the segmentation of four datasets usedin the experiments. The details of the datasets are described in
section 4.1. Figure 16 shows the result of the average Hamming
distance of the four datasets used in the experiments.
To support the evaluation, we used three random baseline
methods (A-C). The random baselines were computed given thenumber of segments as follows: a sequence is divided into thespecified number of segments by using a uniform distribution,and classes of the segments are randomly sampled from theuniform distribution. The random baselines (A-C) represent theresults of baseline segmentation using the correct number ofFrontiers in Robotics and AI | www.frontiersin.org

FIGURE 15 | Variation of F-measure with the error range.

segments, double the number, and half the number, respectively.
The sequences are divided by iterating this procedure 100 times,
and the values in the figures represent the averages of the 100
segmentation trials. As shown in Figure 15, in a smaller error
range, the F-measure of the random baseline (B) is greater than

10

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

TABLE 1 | Segmentation results for the chicken dance.
Hamming Precision Recall
F# of estimateddistancemeasureclasses
HVGH

0.23

0.86

0.86

0.86

4

VAE +
HDP-GP-HSMM

0.31

1.0

0.71

0.83

4

VAE + HDP-HMM

0.73

0.15

1.0

0.26

11

VAE +
HDP-HMM+NPYLM

0.74

0.85

0.79

0.81

15

VAE + BP-HMM

0.33

1.0

0.86

0.92

3

VAE + Autoplait

0.66

0.0

0.0

0.0

1

Random baseline (A)

0.55

0.55

0.41

0.47

4

Random baseline (B)

0.51

0.45

0.69

0.54

4

Random baseline (C)

0.60

0.64

0.20

0.30

3

FIGURE 16 | Hamming distance of each method.
TABLE 2 | Segmentation results for the teapot dance.
Hamming Precision Recall
F# of estimateddistancemeasureclassesthat of the random baseline (A). This is because the number ofboundary points of the random baseline (B) is greater than thatof the random baseline (A), the more boundary points of (B)are likely to be within the error range, and TP increases. On theother hand, in a larger error range, the F-measure of the randombaseline (B) is less than that of the random baseline (A). This isbecause the more boundary points of the random baseline (A) arealso within the error range in the case of a larger error range, andTP increases. Moreover, the precision of the random baseline (B)decreases and recall increases with an increase in FP because thenumber of boundary points is greater than that for the randombaseline (A). In contrast, the F-measure of the random baseline(C) is less than that of the random baseline (A). This is because
precision increases and recall decreases with an increase in FN.
In the case of the percentages of error range where the F-measure
is saturated, the F-measure is lower in both cases in which thenumber of segments is larger or smaller because of the trade-off
relationship between recall and precision. These results indicatethat F-measure reflects the accuracy of boundary points as wellas the correctness of the number of segments. From Figure 15,
we can see that the F-measure begins to saturate at the 5% error
range in all methods except for the random baselines; therefore,
we use a 5% error range in the subsequent experiments.

0.31

0.74

0.83

0.79

7

VAE +
HDP-GP-HSMM

0.36

0.80

0.71

0.75

8

VAE + HDP-HMM

0.75

0.10

1.0

0.17

16

VAE +
HDP-HMM+NPYLM

0.61

0.58

1.0

0.74

21

VAE + BP-HMM

0.34

0.50

0.86

0.63

10

VAE + Autoplait

0.75

0.0

0.0

0.0

1

Random baseline (A)

0.59

0.54

0.46

0.49

6

Random baseline (B)

0.59

0.42

0.72

0.53

7

Random baseline (C)

0.65

0.59

0.21

0.31

5

TABLE 3 | Segmentation results for the exercise motion 1: subject 13.
Hamming Precision Recall
F# of estimateddistancemeasureclasses

4.3. Segmentation of Motion-capture DataFirst, we applied baseline methods to the 93-dimensional
time-series data. However, the baseline methods were
not able to segment the 93-dimensional time-series dataappropriately, because of high dimensionality. Therefore, weapplied the VAE with the same parameters as HVGH, andsequences of three-dimensional latent variables were used forsegmentation of the baseline methods. Tables 1–4 show theresults of segmentation on each of the four motion-capture
datasets. The random baselines in the tables indicate theresults of the random segmentation, which are describedin section 4.2.
VAE+HDP-GP-HSMM and VAE+BP-HMM were able tosegment the motion-capture data from the chicken dance andteapot dance. However, in the results with exercise motion
Frontiers in Robotics and AI | www.frontiersin.org

HVGH

HVGH

0.27

0.66

0.93

0.75

14

VAE +
HDP-GP-HSMM

0.41

0.53

0.93

0.67

11

VAE + HDP-HMM

0.79

0.05

1.0

0.09

10

VAE +
HDP-HMM+NPYLM

0.76

0.32

1.0

0.48

34

VAE + BP-HMM

0.57

0.29

1.0

0.45

7

VAE + Autoplait

0.76

0.0

0.0

0.0

2

Random baseline (A)

0.60

0.45

0.38

0.41

7

Random baseline (B)

0.59

0.36

0.62

0.46

7

Random baseline (C)

0.64

0.51

0.21

0.29

5

obtained using VAE+BP-HMM, the value of the normalizedHamming distance was larger and the F-measure was smaller
than those for the dance motions. This is because simple anddiscriminative motions were repeated in the chicken dance andteapot dance. Therefore, BP-HMM, which is an HMM-basedmodel, was able to segment them. In contrast, the Gaussian
process used in HVGH and HDP-GP-HSMM is non-parametric,
11

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGHmotion. However, it is reasonable to estimate the stationary
state as a unit of motion. Further, dividing a symmetricalmotion into two classes was not erroneous, because the observed
values for the left- and right-sided motion were different.
Therefore, we conclude that HVGH yielded better results inthis case.
Figure 17 illustrates the segmentation results for exercise
motion 2. In this graph, the horizontal axis represents time
steps, the color reflects motion classes, and the top bar is theground truth of the segmentation. It is clear that the segments
and their classes estimated by HVGH are the most similar to theground truth.
Moreover, we compared the VAE with other dimensional
compression methods in HDP-GP-HSMM. Table 5 presents theresults of the segmentation of exercise motion 2 obtained using
the methods in which dimensional compression is performed
through principal component analysis (PCA) (Pearson, 1901)
and independent component analysis (ICA) (Hyvärinen et al.,2000) instead of VAE. PCA and ICA are generally used for
dimensional compression. We used the general FastICA4 as the
ICA implementation, and their three-dimensional output was

making it possible to represent complicated motion patterns
in the exercise data. Moreover, HVGH achieved more accuratesegmentation than HDP-GP-HSMM. We believe that this isbecause the appropriate latent space for the segmentation wasconstructed by using the predictive distribution of the GP as theprior distribution of the VAE in HVGH.
Furthermore, the number of motion classes in the chicken
dance and teapot dance was correctly estimated by HVGH. In
the exercise motion, larger numbers were estimated because
their sequences included complicated motions. In the case of
exercise motion 1, 14 classes were estimated by HVGH—more
than the correct number seven. This is because the stationary
state was estimated as a unit of motion, and symmetricalmotion was separately classified as left-sided and right-sided
motion in different classes. Moreover, 13 classes—more thanthe correct number 11—were estimated by HVGH in exercise
motion 2. Again, this is because stationary motion was estimatedas one motion and because the symmetrical motion shown in
Figure 13J was divided into two classes: left- and right-sided

TABLE 4 | Segmentation results for the exercise motion 2: subject 14.
Hamming Precision Recall
F# of estimateddistancemeasureclasses
HVGH

0.23

0.50

1.0

0.66

13

VAE +
HDP-GP-HSMM

0.39

0.46

1.0

0.63

14

VAE + HDP-HMM

0.82

0.03

1.0

0.07

25

VAE +
HDP-HMM+NPYLM

0.75

0.23

0.86

0.36

42

VAE + BP-HMM

0.79

0.48

0.81

0.55

4

VAE + Autoplait

0.79

0.0

0.0

0.0

3

Random baseline (A)

0.62

0.46

0.41

0.43

9

Random baseline (B)

0.63

0.37

0.66

0.47

11

Random baseline (C)

0.66

0.51

0.19

0.28

7

TABLE 5 | Segmentation results of comparison with other compression methods.
Hamming Precision Recall
F# of estimateddistancemeasureclasses
HVGH

0.23

0.50

1.0

0.66

13

VAE +
HDP-GP-HSMM

0.39

0.46

1.0

0.63

14

PCA +
HDP-GP-HSMM

0.49

0.44

0.77

0.56

7

ICA +
HDP-GP-HSMM

0.56

0.48

0.66

0.54

7

4 https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.
FastICA.html

FIGURE 17 | Segmentation results for exercise motion 2.

Frontiers in Robotics and AI | www.frontiersin.org

12

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

A

B

C

FIGURE 18 | Latent space of the VAE: (A–C) respectively represent the first and second, first and third, and second and third dimension of the latent variables. The
color of each point, which is latent variable reflects the correct motion class.

A

B

C

FIGURE 19 | Latent space of the HVGH: (A–C) respectively represent the first and second, first and third, and second and third dimension of the latent variables. The
color of each point, which is latent variable reflects the correct motion class.

high-dimensional output and low-dimensional latent variables.
Thus, we qualitatively analyze them by comparing the latentvariables with motions, where each dimension roughly represents
the following motions:used identical to the latent variables of VAE. In the case of PCA
and ICA, the min-max normalization, in which the values arenormalized to a range from –1 to 1, was applied for the sequenceof latent variables, as with Nagano et al. (2018).
Additionally, we described the results of HVGH
and VAE+HDP-GP-HSMM for comparison. PCA andICA compress the high-dimensional data into lowdimensional data by linear transform, and complicatedhigh-dimensional time-series data cannot be converted
into low-dimensional data appropriate for segmentationby using these methods. On the other hand, VAE nonlinearly compresses high-dimensional data by using a neuralnetwork and enables the conversion of the high-dimensional
time-series data into low-dimensional data appropriate
for segmentation.
With regard to exercise motion 1, Figure 18 shows thelatent variables estimated by the VAE, and Figure 19 shows
the latent variables learned by mutual learning with HVGH. In
these figures (a-c), respectively, represent the first and second,
first and third, and second and third dimensions of the latentvariables. The color of each point reflects the correct motion
class. In Figure 18, latent variables do not necessarily reflect
the motion class, because they were estimated with the VAE
exclusively. In contrast, in Figure 19, the latent variables in thesame class have more similar values. This means that latent spaceis estimated to represent the features of unit motions. However,
it is difficult to fully understand the meaning of latent spacebecause the VAE represents a non-linear relationship between
Frontiers in Robotics and AI | www.frontiersin.org• Positive region in the first dimension: right and left motionsof the arms
• Negative region in the first dimension: up and down motionsof the arms
• Positive region in the second dimension: twisting of theupper body to the left
• Negative region in the second dimension: twisting of theupper body to the right• Positive region in the third dimension: motions of the legs
• Negative region in the third dimension: bending of theupper body
From these results, we conclude that HVGH can estimate
the correct number of classes and accurate segmentsfrom high-dimensional data by using the proposed mutual
learning loop.

5. CONCLUSION
In this paper, we proposed HVGH, which segments highdimensional time-series data by mutual learning of a VAE and
HDP-GP-HSMM. In the proposed method, high-dimensional
vectors are converted into low-dimensional latent variables
13

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGHrepresenting features of unit motions with the VAE. Using theselatent variables, HVGH achieves accurate segmentation. The
experimental results showed that the segments, their classes, andthe number of classes could be estimated correctly using the
proposed method. Moreover, the results showed that HVGH iseffective with various types of high-dimensional time-series datacompared to a model where the VAE and HDP-GP-HSMM areused independently.
However, the computational cost of HVGH is very highbecause it takes O(N 3 ) to learn N data points using a Gaussian
process, and this is repeated in the mutual learning loop.
Because of this problem, HVGH cannot be applied to largescale time-series data. We plan to reduce the computational
cost by introducing the approximation method for theGaussian process proposed in Nguyen-Tuong et al. (2009)
and Okadome et al. (2014).
Moreover, to simplify the computation, we assumed thatthe dimensions of the observation were independent, and we
consider this assumption reasonable because the experimental
results demonstrated that the proposed method works well.
However, the dimensions of the observation are not actuallyindependent, and the dependency between the dimensions willneed to be considered to model more complicated whole-body
motion. We believe that multi-output Gaussian processes canbe used to represent dependencies between dimensions (Álvarez
et al., 2010; Dürichen et al., 2014).
In HVGH, we do not consider the time warp of the timeseries data. Although GP can deal with motions whose speeds

are slightly different by estimating the variance, motions whose
speeds are considerably different are classified into different
classes. Therefore, we will investigate the robustness of HVGH
against time warp in a future study and extend it to a method
considering time warp.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data
can be found here: http://mocap.cs.cmu.edu/.

AUTHOR CONTRIBUTIONS
MN, TNak, TNag, and DM conceived, designed, anddeveloped the research. MN and TNak performed
the experiment and analyzed the data. MN wrote
the manuscript with support from TNak, TNag, DM,
and IK. DM, IK, and WT supervised the project. All
authors discussed the results and contributed to thefinal manuscript.

ACKNOWLEDGMENTS
This work was supported by JST CREST Grant Number
JPMJCR15E3 and JSPS KAKENHI Grant Number JP18H03295.
We thank IEEE for permission to reuse figures in(Nagano et al., 2018).

REFERENCES

Johnson, M. J., Duvenaud, D. K., Wiltschko, A., Adams, R. P., and Datta, S. R.
(2016). “Composing graphical models with neural networks for structured
representations and fast inference,” in Advances in Neural Information
Processing Systems (Barcelona), 2946–2954.
Kingma, D. P., and Ba, J. (2017). ADAM: a Method for stochastic optimization.
arXiv preprint arXiv:1412.6980v9.
Kingma, D. P., and Welling, M. (2013). Auto-encoding variational bayes. CoRR,
abs/1312.6114.
Lin, J. F. S., and Kulić, D. (2012). “Segmenting human motion for automatedrehabilitation exercise analysis,” in Conference of IEEE Engineering in Medicine
and Biology Society (San Diego, CA), 2881–2884.
Lioutikov, R., Neumann, G., Maeda, G., and Peters, J. (2015). “Probabilistic
segmentation applied to an assembly task,” in IEEE-RAS International
Conference on Humanoid Robots (Seoul), 533–540.
Liu, S., Yamada, M., Collier, N., and Sugiyama, M. (2013). Changepoint detection in time-series data by relative density-ratio estimation.
arXiv[Preprint].arXiv:1203.0453v2
Lund, R., Wang, X. L., Lu, Q. Q., Reeves, J., Gallagher, C., and Feng, Y. (2007).
Changepoint detection in periodic and autocorrelated time series. J. Climate
20, 5178–5190. doi: 10.1175/JCLI4291.1
MacKay, D. J. (1998). “Introduction to Gaussian processes,” in NATO ASI Series F
Computer and Systems Sciences, Vol. 168 (Springer-Verlag), 133–166.
Matsubara, Y., Sakurai, Y., and Faloutsos, C. (2014). “Autoplait: automatic miningof co-evolving time sequences,” in ACM SIGMOD International Conference on
Management of Data (Snowbird, UT), 193–204.
Nagano, M., Nakamura, T., Nagai, T., Mochihashi, D., Kobayashi, I., and Kaneko,
M. (2018). “Sequence pattern extraction by segmenting time series data using
GP-HSMM with hierarchical dirichlet process,” in IEEE/RSJ International
Conference on Intelligent Robots and Systems (Madrid), 4067–4074.
Nguyen-Tuong, D., Seeger, M., and Peters, J. (2009). “Local Gaussian
process regression for real time online model learning and control,” inAdvances in Neural Information Processing Systems 21 (Vancouver, BC),
1193–1200.

Álvarez, M., Luengo, D., Titsias, M., and Lawrence, N. D. (2010). “Efficient
multioutput Gaussian processes through variational inducing kernels,” inProceedings of the Thirteenth International Conference on Artificial Intelligenceand Statistics (Sardinia), 25–32.
Beal, M. J., Ghahramani, Z., and Rasmussen, C. E. (2002). “The infinite
hidden Markov model,” in Advances in Neural Information Processing Systems(Vancouver, BC), 577–584.
Bishop, C. M. (2006). Pattern Recognition and Machine Learning. New York, NY:
Springer, 291–319.
Dürichen, R., Pimentel, M. A., Clifton, L., Schweikard, A., andClifton, D. A. (2014). Multitask Gaussian processes for multivariate
physiological time-series analysis. IEEE Trans. Biomed. Eng. 62, 314–322.
doi: 10.1109/TBME.2014.2351376
Fod, A., Matarić, M. J., and Jenkins, O. C. (2002). Automated derivationof primitives for movement classification. Auton. Robots 12, 39–54.
doi: 10.1023/A:1013254724861
Fox, E. B., Sudderth, E. B., Jordan, M. I., and Willsky, A. S. (2007). The sticky HDPHMM: Bayesian nonparametric hidden Markov models with persistent states.
in Arxiv preprint.
Fox, E. B., Sudderth, E. B., Jordan, M. I., and Willsky, A. S. (2011). Joint modeling ofmultiple related time series via the beta process. arXiv preprint arXiv:1111.4226.
Haber, D., Thomik, A. A., and Faisal, A. A. (2014). “Unsupervised Time Series
Segmentation for High-Dimensional Body Sensor Network Data Streams,”
in IEEE International Conference on Wearable and Implantable Body Sensor
Networks (Zurich), 121–126.
Hyvärinen,
A.,
and
Oja,
E.
(2000).
Independent
component
analysis: algorithms and applications. Neural Netw. 13, 411–430.
doi: 10.1016/S0893-6080(00)00026-5
Jensen, C. S., Kjærulff, U., and Kong, A. (1995). Blocking Gibbs sampling in verylarge probabilistic expert systems. Int. J. Hum. Comput. Stud.. 42, 647–666.
doi: 10.1006/ijhc.1995.1029

Frontiers in Robotics and AI | www.frontiersin.org

14

November 2019 | Volume 6 | Article 115

Nagano et al.

HVGH

Uchiumi, K., Tsukahara, H., and Mochihashi, D. (2015). “Inducing Word and Partof-Speech with Pitman-Yor Hidden Semi-Markov Models,” in Proceedings of
the 53rd Annual Meeting of the Association for Computational Linguistics andthe 7th International Joint Conference on Natural Language Processing (Beijing),
1774–1782.
Van Gael, J., Saatci, Y., Teh, Y. W., and Ghahramani, Z. (2008). “Beam Sampling
for the Infinite Hidden Markov Model,” in International Conference on Machine
Learning (Helsinki), 1088–1095.
Wächter, M., and Asfour, T. (2015). “Hierarchical segmentation of manipulation
actions based on object relations and motion characteristics,” in International
Conference on Advanced Robotics (Istanbul), 549–556.
Yamanishi, K., and Takeuchi, J. I. (2002). “A unifying framework for detecting
outliers and change points from non-stationary time series data,” in ACM
SIGKDD International Conference on Knowledge Discovery and Data Mining(Edmonton, AB), 676–681.

Okadome, Y., Urai, K., Nakamura, Y., Yomo, T., and Ishiguro, H.
(2014). Adaptive LSH based on the particle swarm method with theattractor selection model for fast approximation of Gaussian process
regression. Artif. Life Robot. 19, 220–226. doi: 10.1007/s10015-0140161-1
Pearson, K. (1901). LIII. On lines and planes of closest fit to systems ofpoints in space. Philos. Magaz. J. Sci. 2, 559–572. doi: 10.1080/14786440109
462720
Pitman, J. (2002). Poisson-Dirichlet and GEM invariant distributions for splitand-merge transformations of an interval partition. Combin. Probabil. Comput.
11,501–514. doi: 10.1017/S0963548302005163
Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Stat. Sin. 4,
639–650.
Shiratori, T., Nakazawa, A., and Ikeuchi, K. (2004). “Detecting dance motion
structure through music analysis,” in IEEE International Conference on
Automatic Face and Gesture Recognition (Seoul), 857–862.
Takano, W., and Nakamura, Y. (2016). Real-time unsupervised segmentationof human whole-body motion and its application to humanoid robot
acquisition of motion symbols. Robot. Auton. Syst 75, 260–272.
doi: 10.1016/j.robot.2015.09.021
Taniguchi, T., and Nagasaka, S. (2011). “Double articulation analyzer for
unsegmented human motion using Pitman–Yor language model and infinite
hidden Markov model,” in IEEE/SICE International Symposium on System
Integration (Kyoto), 250–255.
Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M. (2006).
Hierarchical Dirichlet processes. J. Am. Stat. Assoc. 101, 1566–1581.
doi: 10.1198/016214506000000302

Frontiers in Robotics and AI | www.frontiersin.org

Conflict of Interest: The authors declare that the research was conducted in theabsence of any commercial or financial relationships that could be construed as apotential conflict of interest.
Copyright © 2019 Nagano, Nakamura, Nagai, Mochihashi, Kobayashi and Takano.
This is an open-access article distributed under the terms of the Creative Commons
Attribution License (CC BY). The use, distribution or reproduction in other forums
is permitted, provided the original author(s) and the copyright owner(s) are credited
and that the original publication in this journal is cited, in accordance with accepted
academic practice. No use, distribution or reproduction is permitted which does notcomply with these terms.

15

November 2019 | Volume 6 | Article 115