Sequence Pattern Extraction by Segmenting Time Series Data
Using GP-HSMM with Hierarchical Dirichlet Process
Masatoshi Nagano1,∗ , Tomoaki Nakamura1 , Takayuki Nagai1 , Daichi Mochihashi2 ,Ichiro Kobayashi3 and Masahide Kaneko1

I. I NTRODUCTION
Humans perceive continuous information by dividing it
into significant segments. For example, we can recognize
words by segmenting a speech wave and perceive unit motions by segmenting continuous motion. Humans can learn
words and unit motions by appropriately segmenting thecontinuous information without explicit segmentation points,
even though there are an infinite number of possibilities forthe segmentation points. We consider that such an ability isalso important for robots so that they can learn motions andlanguage flexibly.
To this end, we previously proposed a Gaussian processbased hidden semi-Markov model (GP-HSMM). This model
is a hidden semi-Markov model (HSMM) with a Gaussian
process emission distribution, which makes it possible tosegment the time-series data in an unsupervised manner [1].
Fig. 1 depicts the generative process of the GP-HSMM.
Each state of the HSMM represents a class of a segment,
and the segments are generated from the Gaussian process
corresponding to their class. Similar motion patterns are
1 Department of Mechanical Engineering and Intelligent Systems, The
University of Electro-Communications, Chofu-shi, Japan
2 Department of Mathematical Analysis and Statistical Inference, Instituteof Statistical Mathematics, Tachikawa, Japan
3 Department of Information Sciences, Faculty of Sciences, Ochanomizu
University, Bunkyo-ku, Japan
Correspondence∗ : n1832072@edu.cc.uec.ac.jp

Classes(Gaussian Process)Observation value

Abstract— Humans recognize perceived continuous information by dividing it into significant segments such as words andunit motions. We believe that such unsupervised segmentationis also an important ability that robots need to learn topics suchas language and motions. Hence, in this paper, we propose amethod for dividing continuous time-series data into segments
in an unsupervised manner. To this end, we proposed a method
based on a hidden semi-Markov model with Gaussian process(GP-HSMM). If Gaussian processes, which are nonparametric models, are used, unit motion patterns can be extracted
from complicated continuous motion. However, this approachrequires the number of classes of segments in the time-series
data in advance. To overcome this problem, in this paper,
we extend GP-HSMM to a nonparametric Bayesian model byintroducing a hierarchical Dirichlet process (HDP) and propose
the hierarchical Dirichlet processes-Gaussian process-hidden
semi-Markov model (HDP-GP-HSMM). In the nonparametric
Bayesian model, an infinite number of classes is assumed andit becomes difficult to estimate the parameters naively. Instead,the parameters of the proposed HDP-GP-HSMM are estimatedby applying slice sampling. In the experiments, we use various
synthetic and motion-capture data to show that our proposedmodel can estimate a more correct number of classes andachieve more accurate segmentation than baseline methods.

time

Fig. 1.

Overview of the generative process of the proposed method.

generated from the same Gaussian process, and the observed
continuous sequence is generated by combining them. By theestimating the parameters of the GP-HSMM, the time-series
data can be segmented.
However, a GP-HSMM requires the number of classes
into which the continuous sequence is divided. This is a big
drawback because the number of classes is usually not known
in advance. In this paper, to overcome this problem, we extend GP-HSMM to the nonparametric Bayesian model HDPGP-HSMM by introducing a hierarchical Dirichlet process(HDP) [2], which makes it possible to estimate the numberof classes. The nonparametric Bayesian model assumes aninfinite number of classes; however, a finite number of theseclasses are actually used for representing the training data.
As a result, the number of classes can be estimated. By estimating the parameters of the HDP-GP-HSMM, the segments,
their class, and the number of classes can be estimated at thesame time. However, forward filtering–backward sampling,which is the parameter estimation method of GP-HSMM,
cannot be naively applied to HDP-GP-HSMM because thismodel assumes an infinite number of classes. To apply it to
HDP-GP-HSMM parameter estimation, we utilize the slice
sampler proposed in [3] to limit the potential number ofclasses. We evaluate the efficiency of the proposed HDPGP-HSMM using synthetic data and real motion-capture
data. The experimental results show that HDP-GP-HSMM
achieves a higher accuracy than some baseline methods.
Many studies for unsupervised motion segmentation havebeen conducted. However, heuristic assumptions are usedin many of them [4], [5], [6]. Wächter et al. proposed amethod for the segmentation of object manipulation motions
by robots and used contact between the end-effector and theobject as a segmentation clue. The method proposed by Lioutikov et al. requires candidates for the segmentation points

１𝑣1

1 − 𝑣1"𝑣2

1 − 𝑣2𝑣3

０

１

1 − 𝑣3𝑣4

Beta(1, 𝛾)

1 − 𝑣4

…

!
(a)

𝜷

"

"

・・・

Fig. 3.

!
(b)

Stick Breaking Process (SBP).

!
(c)

𝜷
𝝅𝒄 ~𝑫𝑷(𝜼, 𝜷)Fig. 2. Illustration of representation of trajectory: (a) Observed data points,(b) representation by HMM and (c) representation by GP-HSMM.

in advance. In addition, the method proposed by Takano etal. utilizes the errors between the predicted and actual values
as criteria for segmentation. Moreover, some methods use
motion features such as the zero velocity of joint angles [7],[8], [9]. However, it is not clear that these heuristics are valid
for all motions. In contrast, we do not use such heuristics andformulate the problem of the segmentation using a stochastic
model. In some studies, the segmentation is formulated
stochastically, and hidden Markov models (HMMs) are used[10], [11], [12], [13]. However, it is difficult for the HMMs
to represent complicated motion patterns, and we instead
use Gaussian processes, which are a type of non-parametricmodel that can represent more complicated time-series data
than HMMs.
Regarding to the HMMs, HMM using HDP is also proposed [10], [14]. Beal et al. proposed HDP-HMM that has
infinite number of states [10]. Moreover, Fox et al. extended
it to sticky HDP-HMM [14], which can prevent over-segment
by increasing self-transition probability. Fig. 2(b) depicts
the representation of trajectory of data points (Fig. 2(a))
by HMM. HMM represents trajectory using five Gaussian
distributions and the number of Gaussian distributions can be
estimated by using HDP. However, one can see that detailsof the trajectory are lost. On the other hand, in GP-HSMM,
trajectory can be represented continuously by using two GPs
and, moreover, HDP makes it possible to estimate the number
of GPs. Thus, we consider HDP-GP-HSMM can represents
complex trajectories more accurately.
II. N ONPARAMETRIC BAYESIAN MODEL
In this study, the number of motion classes are estimatedby utilizing a non-parametric Bayesian model. In the parametric Bayesian model case, the number of classes is finite,
the parameters of a multinomial distribution are determinedby a finite-dimensional Dirichlet distribution, and, finally,
the classes are sampled from this multinomial distribution.
In contrast, in the non-parametric Bayesian model case, an
infinite number of classes are assumed, and the classes are

𝝅𝒄𝟐

𝝅𝒄𝟏
Fig. 4.

・・・

𝝅𝒄∞

Hierarchical Dirichlet Process (HDP).

sampled from the infinite-dimensional multinomial distribution. Only a finite number of classes are actually used,
and, therefore, their number can be estimated. In the nonparametric Bayesian model, an infinite-dimensional multinomial distribution must be constructed, and one of the methods
for it is the stick breaking process (SBP).
A. Stick Breaking Process (SBP)The SBP [3] is a Dirichlet process that can construct an
infinite-dimensional multinomial distribution. In the SBP, theparameters of infinite-dimensional multinomial distribution
β are constructed by repeatedly breaking a stick, the lengthof which is one, with ratio vk sampled from a beta distribution, as follows (Fig. 3):
vk

∼

βk

=

Beta(1, γ)
k−1∏vk(1 − vi )(k = 1, · · · , ∞),(1)(k = 1, · · · , ∞).

(2)

i=1

This process is represented as β ∼ GEM(γ), where GEM
stands for Griffiths, Engen and McCloskey [15].
B. Hierarchical Dirichlet Process (HDP)
Using a Dirichlet process, HMMs with countably manystates can be constructed. However, if the SBP is naively
used, the destinations reachable from each state have different states. In the case of models like HMMs, all state have toshare the destination states and have different probabilities oftransitioning to each state. To construct such a distribution,
the distribution β that was generated by the SBP is sharedwith all states as a base measure, and transition probabilities
π c , which are different in each state c, are generated byanother Dirichlet process.
π c ∼ DP(η, β).

(3)

Algorithm 1 Blocked Gibbs Sampler
𝜂Fig. 5.

Graphical model of HDP-GP-HSMM.

The method in which the probability distribution is constructed by a two-phase Dirichlet process is called an HDP(Fig. 4).
III. HDP-GP-HSMM
Fig. 5 shows a graphical model of the proposed method,which is a generative model of the time-series data. Here,cj (j = 1, 2, · · · , ∞) denotes the classes of the segments,
and we assumed that the number of classes is countably
infinite. Probability π c denotes the transition probability,
which is generated from the β generated by the GEM
distribution (SBP) parameterized by γ and the Dirichlet
process parameterized by η.
β∼GEM(γ),

(4)

πc

∼

DP(η, β).

(5)Class cj of the j-th segment is determined by the class ofthe (j − 1)-th segment and transition probability π c . The
segment is generated by a Gaussian process with parameterX c ascj∼

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

(6)

xj

∼

GP(x|X cj ),(7)where X c denotes a set of segments that are classified intoclass c. Observed time series data is assumed to be generated
by connecting segments sampled by the above generative
process.
A. Gaussian process
In this paper, each class represents a continuous trajectory
found by learning the emission xi of time step i using a
Gaussian process. In the Gaussian process, given the pairs(i, X c ) of time step i and its emission, which are classified
into class c, the predictive distribution of xnew of time step
inew becomes the Gaussian distributionp(xnew |inew , X c , i) ∝ N (kT C −1 i, c − kT C −1 k),(8)

where k(·, ·) denotes the kernel function and C is a matrix
whose elements areC(ip , iq ) = k(ip , iq ) + ω −1 δpq .

(9)

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 xnj from X cnj
7:
end for8:
9:// Sampling segments and their classes
10:xn∗ , cn∗ ∼ P (x|sn )11:
12:
for j = 1 to Jn do
13:
Ncnj + +14:
Ncnj ,cn,j+1 + +15:Append segments xnj to X cnj
16:
end for17: end forHere, ω denotes a hyper parameter that represents noise in
the observation. Additionally, k is a vector whose elements
are k(ip , inew ), and c is k(inew , inew ). A Gaussian process
can represent complicated time-series data owing to thekernel function. In this paper, we used the following kernel
function, which is generally used for Gaussian processes:
1k(ip , iq ) = θ0 exp(− θ1 ||ip − iq ||2 + θ2 + θ3 ip iq ), (10)2where θ∗ denotes the parameters of the kernel.
Moreover, if the observations are composed of multidimensional vectors, we assume that each dimension is independently generated and, further, the predictive distribution
is such that the emission x = (x0 , x1 , · · ·) of time stepi is generated by a Gaussian process of class c, which iscomputed as follows:
GP(x|X c ) =p(x0 |i, X c,0 , I)× p(x1 |i, X c,1 , I)× p(x2 |i, X c,2 , I) · · · .

(11)

Using this probability, the time-series data can be classified
into the classes.
B. Parameter Inference
1) Blocked Gibbs Sampler: The parameters of the proposed model can be estimated by sampling segments andtheir classes. For efficient estimation, we utilize a blocked
Gibbs sampler [16] in which all segments and their classes
in the observed sequence are sampled. First, all observed
sequences are randomly divided into segments and classified
into classes. Next, segments xnj (j = 1, 2, · · · , Jn ) obtainedby segmenting sn are excluded from the training data, andthe parameters of Gaussian processes X c and transition
probabilities P (c|c′ ) are updated. The segments and theirclasses are sampled as follows:(xn,1 , · · · , xn,Jn ), (cn,1 , · · · , cn,Jn ) ∼ P (X, c|sn ).

(12)

Using the sampled segments and their classes, the parameters of Gaussian processes X c and transition probabilities

Algorithm 2 Forward Filtering–Backward Sampling

𝜋𝑐 ′ , 𝑐ucFig. 6.
πc′ ,c .

Slice sampling truncates the number of classes by thresholding

𝛼 𝑡 𝑘=2 𝑐=1
𝑡−2

𝑡−1

𝑡

𝑐
𝑡−5

𝑡−4

𝑡−3

𝑡−2

𝑡−1

𝑡

𝑐

𝑘
𝑘=2

𝑘

𝑡
(a) Example of Forward filtering

Fig. 7.

𝑡(a) Example of Backward sampling

Forward filtering–backward sampling.

P (c|c′ ) are updated. Iterating this procedure, the parameters
can be optimized. Algorithm 1 shows the pseudo code of theblocked Gibbs sampler. In this algorithm, Ncnj , Ncnj ,cn,j+1
is a parameter to compute the transition probability in Eq.
(15). However, we must also compute the probability in Eq.
(12) for all possible patterns of segments and classes, which
is impractical because of the infinite number of classes.
To compute Eq. (12), we utilize forward filtering–backward
sampling and a slice sampler, which are explained in the next
section.
2) Slice Sampler: We assumed that the number of classes
is countably infinite. However, it is difficult to compute Eq.
(12) because c can be infinite possibilities. To overcome thisproblem, we utilize slice sampler that limits the number ofclasses stochastically proposed in [3]. In the slice sampler,
auxiliary variable uj for each time step j is introduced,
and by thresholding transition probability π c by using uj ,the number of classes become finite. The parameter uj issampled from following probability:p(uj |cj−1 , cj ) =ξ(0 < uj < πcj−1 ,cj ),πcj−1 ,cj

(13)

where ξ(A) = 1 if condition A is true; otherwise, it is 0.
Fig. 6 illustrates the slice sampler. By truncating the classes
with a transition probability πcj−1 ,cj that is less than uj , thenumber of classes becomes finite.
3) Forward Filtering–Backward Sampling: The numberof classes can be truncated by slice sampling. As a result,forward filtering–backward sampling [17] can be applied tocompute Eq. (12). Forward filtering-backward sampling canbe considered as MCMC version of the forward-backward
algorithm which is a maximum likelihood estimation. Although the forward probabilities can be computed in the samemanner as forward-backward algorithm, the latent variables
are sampled backward based on forward probabilities in the
backward sampling. In the forward filtering, the probability

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:
7: // Forward filtering
8: for t = 1 to T do
9:for k = 1 to K do
10:
for c = 1 to C̄ do
11:Compute α[t][k][c]12:
end for13:
end for
14: end for
15:
16: // Backward sampling
17: t = T, j = 1
18: while t > 0 do
19:k, c ∼ p(xj |st−k:t )α[t][k][c]20:
xj = st−k:t21:
cj = c
22:
t=t−k
23:
j =j+1
24: end while
25: return (xJn , xJn −1 , · · · , x1 ), (cJn , cJn −1 , · · · , c1 )that k samples st−k:k before time step t form a segment ofclass c is as follows:
α[t][k][c]=GP(st−k:k |X c )


K ∑
C̄

∑
×
p(c|c′ )α[t − k][k ′ ][c′ ] ,
′

′
k =1

c =1

(14)
where C̄ denotes the maximum number of classes estimatedby slice sampling and K denotes the maximum length ofsegments. In addition, p(c|c′ , β, α) is the transition probability, which can be computed as follows:Nc′ c + αβc′(15),Nc′ + αwhere Nc′ and Nc′ c represent the number of segmentsof class c and the number of transitions from c′ to c,respectively. In addition, k ′ and c′ are the length and class of
possible segments before st−k:k , and these probabilities are
marginalized out. Moreover, α[t][k][∗] = 0 if t − k < 0, andα[0][0][∗] = 1.0. Equation (14) can be recursively computedfrom α[1][1][∗] using dynamic programming, as shown in
Fig. 7(a). This figure depicts an example of computing
α[t][2][1], which is the probability that the two samples
before t become a segment whose class is c. In this case, thelength is two and, therefore, all the segments with end points
t−2 can potentially transit to this segment. In α[t][2][1], these
possibilities are marginalized out.
Finally the length of segments and their classes can be
sampled by backward sampling from t = T . Fig. 7(b) depicts
an example of backward sampling. The length of segment k1
and its class c1 are sampled from the probabilities α[t][∗][∗].
If k1 = 3, k2 and c2 are sampled from α[t − 3][∗][∗]. Byp(c|c′ , β, α) =

(1)

(3)

(2)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

…

Frames〇×

Estimated result

×

〇

×

×

×

×

×

Evaluation

TN

TP

TN

TN

TN

FN

TN

Ground Truth

×

×

×

〇

×

×

×

×

×

×

〇

TN

TN

FP

t

Fig. 8. Example of segmentation evaluation: TP is assigned to the boundary (2) because estimated boundary is within error range from true boundary,which is ±5% of the sequence length.
1.5

1.5

1.5

1.5

1.5

1.0

1.0

1.0

1.0

1.0

0.5

0.5

0.5

0.5

0.5

0.0

0.0

0.0

0.0

0.0

-0.5

-0.5

-0.5

-0.5

-0.5

-1.0

-1.0

-1.0

-1.0

-1.0

-1.5

-1.5

-1.5

-1.5

-1.5

(a)

Observation value

1.5
1.0
0.5
0.0

-0.5
-1.0
-1.5

1

11

21

31

41

51

61

71

81

91

101

111

121

131

141

Time step(b)

Observation Value

Fig. 9. Synthetic time series 1: (a) five unit sequences and (b) connectedsequence.
1.2

1.2

1.2

1.0

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.0

0.0

1.2
1.0
0.8
0.6
0.4
0.2
0.0

1

11

21

31

41

51

61

0.2
0.0

(a)

71

81

91

101

111

121

(b)

131

141

151

161

Time Step

Fig. 10. Synthetic time series 2: (a) three unit sequences and (b) connectedsequence.

iterating this procedure until t = 0, the segments and theirclasses can be determined.
Algorithm 2 shows the pseudo code of forward filtering–backward sampling with slice sampling.
IV. E XPERIMENTS
To validate the proposed HDP–GP-HSMM, we applied it
to several types of time-series data. For comparison, we used
HDP-HMM [10], HDP-HMM+NPYLM [11], BP-HMM[12],and Autoplait [13] as baseline methods.
A. Evaluation metrics
We used four measures to evaluate the accuracy of segmentation: normalized Hamming distance, precision, recall,and F-measure.
The normalized Hamming distance represents the error
rate of the classification of the data points and is computed
as follows:D(c, c),N D(c, c) =
(16)
|c|

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 thatthe estimated classes are more similar to the ground truth.
To compute precision, recall, and F-measure, we evaluatedboundary points (boundaries between segments) as true positives (TPs), true negatives (TNs), false positives (FPs), andfalse negatives (FNs), as shown in Fig. 8. TP is assigned tothe points that are correctly estimated as boundary points. An
estimated boundary point is treated as correct if the estimatedboundary was within ±5% of sequence length, as shownin Fig. 8, Frame (2). TN is assigned to the points that are
correctly estimated as not boundary points, as shown in Fig.
8, Frame (3). Conversely, FP and FN are assigned to points
that are falsely estimated as boundary points, as shown in
Fig. 8, Frame (10), and falsely estimated as not boundary
points, as shown in Fig 8, Frame (6), respectively. From theseboundary evaluations, the precision, recall, and F-measureare computed as follows:
NT P
,
(17)
NT P + NF P
NT P
R =
(18)
,
NT P + NF N
2P · R
,
F =
(19)
P +R
where NT P , NF P , and NF N represents the number ofboundary points estimated as TP, FP, and FN respectively.
P=

B. Segmentation of synthetic time-series data
First, we generated the following synthetic time-series data
(shown in Figs. 9 and 10), and applied the segmentation
methods to them.
• Synthetic time series 1: Unit sequences are a seriesof constant values randomly selected from a set of fivevalues. The length of each unit sequence is randomly
generated from a uniform distribution.
k ∼ uniform(9, 13).

(20)

Moreover, normal random values are added to eachsequence as follows:x ∼ x + normal(x, 0.05).

(21)

TABLE I
S EGMENTATION RESULTS OF SYNTHETIC TIME SERIES 1.

HDP-GP-HSMM
HDP-HMMHDP-HMM+NPYLM
BP-HMM
AutoplaitHamming
distance
0.011
0.011
0.48
0.36
0.79

Precision
1.0
1.0
0.53
1.0
0.00

Recall
1.0
1.0
1.0
0.94
0.00

F-measure
1.0
1.0
0.69
0.97
0.00

# of estimatedclasses
5
5
35
4
3

(a)

(b)

(c)

(d)

Fig. 11. Four unit motions included in the chicken dance: (a) beaks, (b)wings, (c) tail feathers, and (d) claps.

TABLE IIS EGMENTATION RESULTS OF SYNTHETIC TIME SERIES 2.

HDP-GP-HSMM
HDP-HMMHDP-HMM+NPYLM
BP-HMM
AutoplaitHamming
distance
0.041
0.13
0.45
0.038
0.35

Precision
1.0
0.43
0.57
1.0
0.17

Recall
1.0
1.0
1.0
1.0
0.15

F-measure
1.0
0.60
0.72
1.0
0.16

# of estimatedclasses
3819
32The sequence shown in Fig. 9 is generated by connecting these unit sequences.
• Synthetic time series 2: The sequence shown in Fig.
10(b) is generated by randomly connecting the three
sequences shown in Fig. 10(a).
Tables I and II show the segmentation results of synthetic
time series 1 and 2. The results in these tables show that theHamming distance of the proposed method is less than 0.1
and, moreover, all the precision, recall, and F-measure values
are 1.0. These results indicate that our proposed HDP-GPHSMM is able to estimate the segments almost completely
correctly. Furthermore, the number of classes were correctly
estimated as five in synthetic time series 1 and three insynthetic time series 2.
HDP-HMM was also able to correctly segment synthetic
time series 1. This is because a simple model with a Gaussian
emission distribution is suitable for this time series. However,
HDP-HMM was not able to correctly segment synthetic timeseries 2 which did not follow a Gaussian distribution. In
contrast, HDP-GP-HSMM and BP-HMM, which are more
complicated models, were able to deal with such data.
C. Segmentation of motion-capture data
To evaluate the validity of the proposed method using
more complicated real data, we applied it to the following
three motion capture data series.
• Chicken dance: We used a sequence of motion capture
data of a human performing a chicken dance from the
CMU Graphics Lab Motion Capture Database.1 The
dance includes four motions, as shown in Fig. 11. The
accelerations of both hands and feet, which were usedin the experiment in [13], were also used in theseexperiments.
• Exercise motion: We also used three sequences of
exercise motion capture data from subject 14 in the
CMU Graphics Lab Motion Capture Database. This
motion includes 11 motions, as shown in Fig. 12. The
two-dimensional position (horizontal and vertical) of
both hands and feet were used.
1 http://mocap.cs.cmu.edu/: subject 18, trial 15

(a)

(b)

(c)

(d)

(e)

(g)

(h)

(i)

(j)

(k)

(f)

Fig. 12. Eleven unit motions included in the exercise motion: (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.

(a)

(e)

(c)

(b)

(f)

(d)

(g)

Fig. 13. Examples of motion included in the karate motion: (a) left lower
guard, (b) right punch, (c) right lower guard, (d) left punch, (e) left upper
guard, (f) right upper guard, and (g) preliminary motion.

Karate motion: To include more complicated motions,
we used karate motion capture from mocapdata.com2 .
Four sequences into which it is divided were used. Fig.
13 shows examples of the unit motions of the karate
data. As for the other data series, the two-dimensional
position of both hands were used.
All the data was preprocessed by down sampling andnormalization. In the case of HDP-GP-HSMM, all the data
was downsampled to reduce computational cost, and minmax normalization, in which the values are normalized to arange from -1 to 1, was applied. In the case of other methods,
the segmentation was affected by normalization and down
sampling for each data and, therefore, we choose a better
normalization from min-max and z-score, and a better sample
rate from fifteen fps, eight fps, four fps and two fps. The zscore normalization makes means and standard deviation of
the normalized values zero and one respectively. The usedvalues for down sampling and normalized methods are shownin Table III. HDP-GP-HSMM requires hyper parameters, andwe used λ = 12.0, θ0 = 1.0, θ1 = 1.0, θ2 = 0.0, θ3 = 16.0,
which were determined empirically for segmentation of 4
fps and min-max normalized data. Blocked Gibbs sampler
was repeated 10 times, which was the number of times toconverge the parameters.
•

2 http://www.mocapdata.com/

■ Left lower guard■ Left upper guard■ Right punch■ Right upper guard■ Right lower guard■ Left punch■ Preliminary motion

Ground Truth
HDP-GP-HSMM
HDP-HMMHDP-HMM+NPYLM
BP-HMM
Autoplait

Sequence 1

Sequence 2

Fig. 14.

Sequence 3

Segmentation results for the karate motion.

TABLE IIIVALUES USED FOR PREPROCESSING .
Chicken dance
4 fps / min-max
8 fps / z-score
15 fps / z-score
15 fps / z-score
Original

Exercise motion
4 fps / min-max
2 fps / z-score
2 fps / min-max
15 fps / z-score
15 fps / z-score― Ground Truth

Karate motion
15 fps / min-max
2 fps / min-max
4 fps / min-max
15 fps / z-score
15 fps / z-score

TABLE IVS EGMENTATION RESULTS FOR THE CHICKEN DANCE .

HDP-GP-HSMM
HDP-HMMHDP-HMM+NPYLM
BP-HMM
AutoplaitHamming
distance
0.13
0.46
0.65
0.16
0.026

Precision
1.0
0.39
0.63
0.70
1.0

Recall
0.86
0.71
0.71
1.0
1.0

F-measure
0.92
0.50
0.67
0.82
1.0

# of estimatedclasses
44844TABLE VS EGMENTATION RESULTS FOR THE EXERCISE MOTION .

HDP-GP-HSMM
HDP-HMMHDP-HMM+NPYLM
BP-HMM
AutoplaitHamming
distance
0.31
0.82
0.63
0.23
0.61

Precision
0.38
0.070
0.61
0.25
0.67

Recall
0.95
1.0
1.0
1.0
0.18

F-measure
0.55
0.13
0.76
0.40
0.28

# of estimatedclasses
10
14
26
18
5

TABLE VI
S EGMENTATION RESULTS FOR THE KARATE MOTION .

HDP-GP-HSMM
HDP-HMMHDP-HMM+NPYLM
BP-HMM
AutoplaitHamming
distance
0.32
0.58
0.61
0.39
0.60

Precision
0.48
0.38
0.67
0.34
0.00

Recall
0.85
1.0
0.31
0.75
0.00

F-measure
0.62
0.55
0.42
0.47
0.00

# of estimatedclasses
791073Tables IV, V, and VI show the results of segmentation oneach of the three motion-capture time series. The results inthese tables show that HDP-GP-HSMM is able to accurately
segment all three sets of motion capture data. BP-HMM andAutoplait are able to segment the chicken dance motioncapture data; however, exercise and karate motions are notsegmented well by these methods. This is because simple anddiscriminative motions were repeated in the chicken dance,
which makes this motion the simplest of the three motionsin this experiment. In contrast, the Gaussian process use in
HDP-GP-HSMM is a non-parametric method, which makes

Number of classes

HDP-GP-HSMM
HDP-HMMHDP-HMM+NPYLM
BP-HMM
Autoplait

Sequence 4― HDP-GP-HSMM

16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0
1

11

21

31

41

51

61

71

81

Estimated times

Fig. 15.

Transition of the number of estimated classes.

it possible to represent the complicated motion patterns in
the exercise and karate motion-capture data.
Moreover, the number of chicken-dance and karate-motion
classes could be correctly estimated by HDP-GP-HSMM. In
the case of exercise motion, 10 classes were estimated by
HDP-GP-HSMM, which was less than the correct number
11 because similar unit motions, such as motion F andmotion H, as shown in Fig. 12, were classified as the same
motion class. However, the number of classes estimated by
HDP-GP-HSMM is still closer than those estimated by theother methods. Therefore, we conclude that HDP-GP-HSMM
yielded better results in this case.
Fig. 14 illustrates the segmentation results for the karate
motion. Here, the horizontal axis represents time steps, color
reflects motion classes, and the top bar represents the groundtruth of segmentation. It is clear that the segments and theirclasses estimated by HDP-GP-HSMM are similar to theground truth. Furthermore, Fig. 15 shows the transition of
the number of estimated classes by HDP-GP-HSMM. In thisgraph, the horizontal axis represents the number of iterationsof the blocked Gibbs sampler, which fluctuates greatly in theearlier phase. However, it eventually converges to the correct
value. Using HDP, the number of classes can be estimatedwhile the segments and their classes are being estimated.
V. C ONCLUSION
In this paper, we proposed the HDP-GP-HSMM, whichmakes it possible to segment time-series data and determine
the number of classes by introducing an HDP into an HSMM

whose emission distribution is a Gaussian process. In the
proposed method, a slice sampler truncates the number ofclasses and makes it possible to apply forward filtering–backward sampling. The experimental results showed thatthe segments, their classes, and the number of classes canbe estimated correctly using the proposed method. Moreover,
the results also showed that the HDP-GP-HSMM is effective
on various types of time-series data including synthetic andreal sequences.
However, the computational cost HDP-GP-HSMM is very
high because it takes O(N 3 ) to learn N data points using a
Gaussian process. Because of this problem, HDP-GP-HSMM
cannot be applied to large datasets. We are planning to reducethe computational cost by introducing an approximation
method for the Gaussian process proposed in [18], [19].
Moreover, we used the assumption where dimensions of
the observation were independent to simplify the computation, and we consider this assumption is reasonable because
the experimental results showed the proposed method works
well. However, the dimensions of the observation are notactually independent and we also consider the dependency
between the dimensions should be needed to model more
complicated whole body motions. We consider that multioutput Gaussian processes can be used to represents dependencies between dimensions[20], [21].
ACKNOWLEDGMENTS
This work was supported by JST CREST Grant Number JPMJCR15E3 and JSPS KAKENHI Grant Number
JP17K12758.
R EFERENCES[1] T. Nakamura, T. Nagai, D. Mochihashi, I. Kobayashi, H. Asoh,and M. Kaneko, “Segmenting Continuous Motions with Hidden Semi-Markov Models and Gaussian Processes,” Frontiersin Neurorobotics, vol.11, no.67, 2017.
[2] Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei,“Hierarchical Dirichlet processes,” Journal of the American
Statistical Association, vol.101, no.476, pp.1566–1581, 2006.
[3] J. V. Gael, Y. Saatci, Y. W. Teh, and Z. Ghahremani, “BeamSampling for the Infinite Hidden Markov Model,” International
Conference on Machine Learning, pp.1088–1095, 2008.
[4] M. Wächter and T. Asfour, “Hierarchical segmentation ofmanipulation actions based on object relations and motion characteristics,” International Conference on Advanced
Robotics(ICAR), pp.549–556, 2015.
[5] R. Lioutikov, G. Neumann, G. Maeda, and J. Peters, “Probabilistic segmentation applied to an assembly task,” IEEE-RAS
International Conference on Humanoid Robots, pp.533–540,
2015.
[6] W. Takano, and Y. Nakamura, “Real-time unsupervised segmentation of human whole-body motion and its application tohumanoid robot acquisition of motion symbols,” Robotics andAutonomous Systems, vol.75, pp.260–272, 2016.
[7] A. Fod, M.J. Matarić, and O.C. Jenkins, “Automated derivation of primitives for movement classification,” Autonomous
Robots, vol.12, no.1, pp.39–54, 2002.
[8] T. Shiratori, A. Nakazawa, and K. Ikeuchi, “Detecting dance
motion structure through music analysis,” IEEE International
Conference on Automatic Face and Gesture Recognition,
pp.857–862, 2004.

[9] J. F. S. Lin, and D. Kulić, “Segmenting human motion for automated rehabilitation exercise analysis,” Conference of IEEE
Engineering in Medicine and Biology Society, pp.2881–2884,
2012.
[10] M. J. Beal,, Z. Ghahramani, and C. E. Rasmussen, “Theinfinite hidden Markov model,” In Advances in Neural Information Processing Systems, pp.577–584, 2001.
[11] Taniguchi. T, and Nagasaka. S, “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, pp.250–255,
2011.
[12] E. B. Fox, E. B. Sudderth, M. I. Jordan, and A. S. Willsky,“Joint modeling of multiple related time series via the beta
process,” arXiv preprint arXiv:1111.4226, 2011.
[13] Y. Matsubara, Y. Sakurai, and C. Faloutsos, “Autoplait: Automatic mining of co-evolving time sequences,” ACM SIGMODInternational Conference on Management of Data, pp.193–
204, 2014.
[14] E. B. Fox, E. B. Sudderth, M. I. Jordan and A. S. Willsky, “The sticky HDP-HMM: Bayesian nonparametric hidden
Markov models with persistent states”, Annals of Applied
Statistics, vol. 5, no. 2A, pp. 1020–1056, 2011[15] J. Pitman, “Poisson-Dirichlet and GEM invariant distributions
for splitand-merge transformations of an interval partition,”
Combinatorics, Probability and Computing, vol.11, pp.501–
514, 2002.
[16] Jensen. Claus S, Kjærulff. Uffe and Kong. Augustine, “Blocking Gibbs sampling in very large probabilistic expert systems,”
International Journal of Human-Computer Studies, vol. 42,
no.6, pp.647–666, 1995.
[17] K. Uchiumi, T. Hiroshi, and D. Mochihashi, “Inducing Wordand Part-of-Speech with Pitman-Yor Hidden Semi-Markov
Models,” Joint Conference of the 53rd Annual Meeting of the
Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing,
pp.1774–1782, 2015.
[18] D. Nguyen-Tuong, J. R. Peters, and M. Seeger, “Local Gaussian process regression for real time online model learningand control,” In Advances in Neural Information Processing
Systems, pp.1193–1200, 2009.
[19] Y. Okadome, K. Urai, Y. Nakamura, T. Yomo, and H. Ishiguro,“Adaptive lsh based on the particle swarm method with theattractor selection model for fast approximation of Gaussian process regression,” Artificial Life and Robotics, vol.19,
pp.220–226, 2014.
[20] Álvarez. Mauricio, Luengo. David, Titsias. Michalis andLawrence. Neil, “Efficient multioutput Gaussian processes
through variational inducing kernels,” Proceedings of the Thirteenth International Conference on Artificial Intelligence andStatistics, pp.25–32, 2010.
[21] Dürichen. Robert, Pimentel. Marco AF, Clifton. Lei,
Schweikard. Achim and Clifton. David A, “Multitask Gaussian
processes for multivariate physiological time-series analysis,”IEEE Transactions on Biomedical Engineering, vol.62, no.1,
pp.314–322, 2015.