A Bayesian Hybrid Approach to Unsupervised Time Series Discretization
Yoshitaka Kameya∗ , Gabriel Synnaeve† , Andrei Doncescu‡ , Katsumi Inoue§ , Taisuke Sato∗∗ Tokyo Institute of Technology, Japan, {kameya,sato}@mi.cs.titech.ac.jp† Grenoble University, France, gabriel.synnaeve@gmail.com‡ LAAS-CNRS, France, andrei.doncescu@laas.fr§ National Institute of Informatics, Japan, ki@nii.ac.jp

Abstract—Discretization is a key preprocessing step inknowledge discovery to make raw time series data applicableto symbolic data mining algorithms. To improve the comprehensibility of the mined results, or to help the induction stepof the mining algorithms, in discretization, it is natural toprefer having discrete levels which can be mapped into intuitive
symbols. In this paper, we aim to make smoothing of the data
points along with the time axis, and make binning or clustering
at the measurement axis. In particular, we propose a hybrid
discretization method based on variational Bayes, in which theoutput of one discretization method is smoothly exploited as
hyperparameters of another probabilistic discretization modelsuch as a continuous hidden Markov model. The experimentswith artificial and real datasets exhibit the usefulness of thishybrid approach.
Keywords-discretization; continuous hidden Markov models;variational Bayes

I. INTRODUCTION
In a knowledge discovery process, the data are often
obtained as time series of real-valued measurements by
experiments or from sensors. Discretization [1], [2] is then akey preprocessing step to make these time series applicableto symbolic data mining algorithms, e.g. frequent patternmining or inductive logic programming. In discretization,
to improve the comprehensibility of the mined results, or tohelp the induction step of the mining algorithms, it is naturalto prefer having discrete levels which can be mapped into
intuitive symbols like “high,” “medium,” or “low.” Thesesymbols can also be basic elements in qualitative/logicalreasoning [3], [4].
One frequent problem in time series discretization is thatthe data just include raw measurements and have no extra
clue to know the meaningfulness of a discretization. In thiscase, discretization needs to be performed in an unsupervised
fashion, which has been less explored [2], and such unsupervised methods would inherently rely on some assumption
or model behind the data. We may see that discretization ofa time series is a segmentation or clustering process in a
two-dimensional (measurement-time) space, preserving the
temporal behavior of the time series. In this paper, we aim
to make smoothing of the data points along with the timeaxis, and make binning or clustering at the measurement
axis. This strategy basically works, but as we will see later,

it is not always easy to find a desired discretization onlywith a single underlying segmentation/clustering model.
From this background, hybrid approaches are practically attractive where we effectively combine a couple ofpromising heterogeneous methods for smoothing, binning orclustering. In this paper, we propose a hybrid discretization
method based on variational Bayes [5], [6], in which theoutput of one discretization method is smoothly exploitedas hyperparameters (as “prior knowledge”) of another probabilistic discretization model. Also variational Bayes issuitable by itself for unsupervised discretization in that it
is known as robust against noises and provides a principled
way of determining the plausible number of discrete levels(i.e. model selection). Our main tool is continuous (density)
hidden Markov models (HMMs) [7], and the experiments
show that the output of Persist [8] or SAX [9] can be agood guide for building a robust HMM against noisy timeseries.
The remainder of this paper is structured as follows. In
Section II, we describe the previous approaches to unsupervised time series discretization and clarify our motivation.
Section III and IV respectively describe the proposed method
and the experimental results. In Section V, we conclude thepaper and mention the future work.
II. PREVIOUS APPROACHESIn this section, we make a review on unsupervised discretization with a little detailed description on the methods
our hybrid method uses. Also we present a reproduced resultof a comparative experiment by Mörchen and Ultsch [8]and discuss the advantages and disadvantages of the existing
methods to make our motivation clear.
Before starting, let us make some preliminaries. In thispaper, we consider to discretize a univariate time seriesx = (x1 , x2 , . . . , xT ) into x̃ = (k1 , k2 , . . . , kT ) where T
is the length of x, each xt is the measurement at discrete
time t, and kt is a discrete level into which xt is converted(t ∈ [1, T ]). We use positive integers to indicate the discretelevels, and will map a larger measurement into a higher
discrete level. K denotes the number of possible levels(hence kt ∈ [1, K]). Also let µ̂(x), σ̂ 2 (x), and ζq (x) denotethe sample mean, the sample variance and the q-percentile

value (e.g. ζ50 (x) is the median) of the measurements
appearing in a time series x, respectively.
A. Binning, clustering and smoothing
In the literature of discretization, many algorithms havebeen proposed for supervised situations, but unsupervised
discretization has been less explored [2], and also for thecase of time series discretization. In a simplest way, binningor clustering methods, such as equal width binning, equal
frequency binning, K-means, Gaussian mixture models andso on, can be applied as stand-alone discretizers (e.g. [1],[10]), though the temporal information is lost. On the other
hand, while the temporal characteristics of interest (e.g. periodicity, sharp peaks/valleys, long trends and so on) would
differ according to the purpose, following [8], we focus on
the time series whose flat portions (referred to as “enduring
states” in [8]) are crucial to the goal of the application.
Hence, we aim to reveal the potentially flat portions bysmoothing (by removing noises), and then to map theminto discrete levels. For instance, Geurts’s discretization [11]makes a smoothing based on a regression tree where thetime axis is recursively segmented so that the measurementsin each segmented time interval are as close as possible.
Piecewise aggregate approximation, which will be describedin the next section, is another smoothing method. Also as
done in [10], we may apply smoothing filters like SavitzkyGolay filters in advance to the noisy time series.
B. SAX
One may consider from the above that it is reasonable
to combine the binning/clustering/smoothing methods, sequentially or simultaneously. Symbolic Aggregate approXimation (SAX) [9] is a well-known algorithm in which
the measurements in a time series are first smoothed ineach of equal-width segments (frames) at the time axis bypiecewise aggregate approximation (PAA), and the smoothedmeasurements are then grouped into K equal frequency
bins under the assumption that the measurements follow a
Gaussian distribution.
More formally speaking, in PAA, a raw time series x oflength T is compressed into a time series x̄ = (x̄1 , . . . , x̄T 0 )where T 0 < T . For simplicity, we suppose T = wT 0 where
w is some positive integer which indicatesthe frame width1 .
∑wt0
1
Each x̄t0 is computed as x̄t0 = w t=w(t0 −1)+1 xt , the average of the measurements in the t0 -th segment (t0 ∈ [1, T 0 ]).
Then, SAX classifies x̄t0 into the k-th bin [βk−1 , βk ) where
−1
−1
k0
0
βk0 = Φµ̂(x),σ̂2 (x) ( K ) (k ∈ [1, K], k ∈ [0, K]). Here Φµ,σ 2
is the inverse cumulative distribution function of a Gaussian
distribution N (µ, σ 2 ). Finally we construct a discretized
time series x̃ = (k1 , k2 , . . . , kT 0 ) where x̄t0 ∈ [βkt0 −1 , βkt0 ).
We can see βk ’s (k ∈ [1, K − 1]) as the breakpoints at themeasurement axis, and these breakpoints can also be usedlater as a byproduct.
1 See Section 3.5 of [9] for the case in which T is not divisible by T 0 .

In the original description, SAX takes as input a standardized time series, and then it only depends on K and T 0
above, which are specified by the user. In this sense, SAX
is not data-adaptive. On the other hand, the compression by
PAA brings a good approximation of the original series anda significant improvement in efficiency at the later tasks,
called dimension reduction [9].
C. Persist
The algorithm named Persist [8] makes variable-widthbinning at the measurement axis adaptively to the input
time series. In this algorithm, we first regard the eventthat a measurement xt at time t falls into the k-th bin[βk−1 , βk ) as the event that the time series stays at anenduring state sk at time t (k ∈ [1, K], t ∈ [1, T ]).
Accordingly a time series is considered to move around in
the state space S = {s1 , s2 , . . . , sK }. Then, for a state space
S, we introduce a heuristic score, named the persistence∑K
def
1
score, as Persistence(S) = K
k=1 Persistence(sk ).
Each Persistence(sk ) is defined as sgn(P̂ (sk |sk )− P̂ (sk ))·SKL(P̂ (sk |sk ), P̂ (sk )), where p̂(sk |sk ) is the estimated self
transition probability at the state sk , p̂(sk ) is the estimatedprobability of staying at sk , and SKL(q, q 0 ) is the symmetric Kullback-Leibler divergence between two Bernoulli
distributions {q, 1 − q} and {q 0 , 1 − q 0 }. The persistencescore is based on a Markov process, and measures thetotal persistence (or enduringness) of the state space S withrespect to the input time series. Under this setting, Persist
tries to find a state space, or a set of variable-width bins,
S that makes Persistence(S) as high as possible. After thebins found, similarly to SAX, the time series is discretizedand the bin boundaries or the breakpoints βk can be used
as a byproduct.
The higher persistence leads to the higher smoothnessof the discretized time series, so we can say that Persist
performs binning at the measurement axis and smoothingalong the time axis simultaneously. Also by definition, we
can use the persistence score to measure the plausibility of
the number K of discrete levels. Due to the space limit, weomit the algorithmic details.
D. Continuous hidden Markov models
HMMs are one of the standard tools for analyzing sequence data in various application fields such as speech
recognition [7], [12]. As in [8], continuous HMMs can also
be used for discretization with K levels, in which each hidden state sk corresponds to the k-th discrete level. Similarly
to Persist, we have a state space S = {s1 , s2 , . . . , sK } anda time series is considered to move around in this space. A
continuous HMM has the initial state distribution p(Q1 =sk ), the transition distributions p(Qt+1 = sk0 |Qt = sk ) andthe emission distributions (densities) p(Xt = xt |Qt = sk ),where Qt is a random variable which takes the state at timet and Xt is a random variable which takes the measurement
at time t.
In HMMs considered in this paper2 , the emission distribution at the k-th state p(Xt = xt |Qt = sk ) follows aGaussian distribution N (µk , σk2 ). In our case, this indicatesthat the k-th level has its own mean µk and variance σk2of the measurements. Without loss of generality, we assume
that µk ≤ µk0 if k ≤ k 0 . Discretization is then performedby the Viterbi algorithm [7] which takes as input a raw timeseries x = (x1 , x2 , . . . , xT ) and returns the most probablestate sequence z ∗ = (sk1 , sk2 , . . . , skT ), from which we canconstruct the discretized time series x̃ = (k1 , k2 , . . . , kT ).
Under the maximum likelihood framework, the parameters of the above distributions are estimated from the input
time series x by the forward-backward algorithm [7]. The
forward-backward algorithm can be seen as a temporalextension of probabilistic clustering of the measurements
based on a univariate Gaussian mixture model, and theestimated positions (µ̂k ) and the shapes (σ̂k2 ) of the emission
Gaussian distributions are crucial in discretization. In theViterbi algorithm, on the other hand, the probabilities ofself transitions work as weights to keep the HMM staying
at the same state as long as possible, and consequently playa role of smoothing.
E. Experiment with the enduring-state dataset
In this section, to discuss the characteristics of the existing
discretization methods above, we present a reproduced resultof the comparative study in [8] with an artificial dataset,
which we hereafter call the enduring-state dataset. In addition to accuracy adopted in [8], we introduce normalized
mutual information (NMI) [13] as an evaluation criterion on
predictive performance. NMI is frequently used in evaluationof a clustering result.
In the enduring-state dataset, raw time series of length
1,000 are generated by a state machine which randomlychanges its state after a random duration. See [8] for moredetails. At each state, the data points (the measurements) are
generated with Gaussian noises around the mean proper tothe state. The generation process is thus close to a hidden
Markov process, but additionally, some of the data points
are replaced with outliers. The ratio of these outliers varies
from 0% to 10%. Fig. 1 above shows a time series with fivestates and 5% outliers. As in Persist and HMMs, each statecorresponds to a discrete level. Hereafter the state sequenceobtained in the sampling process of a raw time series iscalled the answer sequence, and the output of a discretizer
is called the predicted sequence.
The goal here is to see how well the discretizers recover
the answer sequence from the noisy time series. We pickup six methods to compare from [8]: equal width binning (EQW), equal frequency binning (EQF), SAX, Persist,
2 HMMs in [7] uses a Gaussian mixture as the emission distribution. Ourversion is a special case with only one mixture component.

Figure 1. (above) An enduring-state time series. The plots with a thickline indicate the answer sequence. (below) Predictive performance of theexisting methods in the enduring-state dataset with five levels.

clustering by Gaussian mixture models (GMM), continuous
HMMs (HMM). We tested these methods on 100 time series
for each number K of discrete levels and ratio R of outliers.
In the forward-backward algorithm for HMMs, the meansand variances of the emission distributions are initialized
as the sample mean µ̂(x) (with small Gaussian noises) andthe sample variance σ̂ 2 (x) of the raw time series x, andwe chose the one achieving the highest likelihood from 100
reinitializations. In SAX, we picked up the frame width from{1, 2, 3, 5, 10, 20, 50} that works best for each pair of K
and R. Fig. 1 below shows the median accuracy (left) andthe median NMI (right) for the time series with five discretelevels and various ratios of outliers. The error bars indicate
the 95% median absolute deviation (MAD) t confidence
interval [14].
Under the accuracy, as reported in [8], a direct use of binning or clustering (EQW, EQF or GMM) works poorly, andSAX’s normality assumption at the measurement axis seems
not to fit to this dataset. Also it is observed that continuous
HMMs work nearly perfectly for the time series without
outliers (this is not surprising from the way of generatingtime series), but their performance quickly degenerates as
the outliers increase. Persist, on the other hand, is foundto be stable (i.e. the variances of accuracies are small) androbust against outliers. Interestingly, the results under NMI
show a different tendency. That is, the degeneration of the
HMMs’ performance is similar to the others, and HMMs
constantly outperform them.
The difference in tendencies comes from the nature of

evaluation criteria and the nature of the discretization methods. The accuracy is defined on exact matchings of the discrete levels between the answer sequence and the predictedsequence (i.e. we check if kt and kt0 are equal or not for eachtime t). So it is crucial under the accuracy to identify (thebreakpoints or the means of) the discrete levels in the answer
sequence. Contrastingly, NMI does not depend on suchexact matchings and just indicates the degree of overlapping
between two groupings of data points (that is why NMI isused for the evaluation of a clustering result). The answer
sequences are of course smooth, and thus NMI emphasizes
the smoothness of the predicted sequence. Indeed, after acloser look into the predicted sequences, we found that,in Persist, many non-outliers (only perturbed by Gaussian
noises) go across the boundaries of the bins, and HMMs
frequently fail to identify the answer discrete levels for thetime series with many outliers. From these observations andthe descriptions in Sections II-C and II-D, we summarize
that, in this dataset, Persist is robust in identifying theunderlying discrete levels thanks to the persistence score’s
global nature, while HMMs are good at local smoothingamong the neighboring measurements. Next, we propose aprincipled way to exploit these promising characteristics atthe same time.
III. HYBRIDIZATION VIA VARIATIONAL BAYESTo exploit the advantages of the heterogeneous discretizers, i.e. Persist (or another method) and continuous HMMs,
we propose a hybrid method based on variational Bayes(VB) [5], [6]. Bayesian learning in general is known to berobust against noises, and VB is a fast approximate Bayesianlearning that makes available a dynamic programming based
procedure such as the forward-backward algorithm. Although there have been a couple of previous works which
apply VB to continuous HMMs [15], [12], to the best of
our knowledge, our hybrid method is the first attempt to
introduce VB to discretization tasks. The rest of this section
first explains how VB is applied to continuous HMMs, andthen, in Section III-D, describes how we hybridize Persist
and HMMs in the VB framework.
A. Model description
In Section II-D, we have explained that a continuous
HMM has the initial state distribution p(Q1 = sk ), thetransition distributions p(Qt+1 = sk0 |Qt = sk ) and theemission distributions p(Xt = xt |Qt = sk ). Hereafter wedenote the state at time t by zt , and abbreviate p(Xt = xt )and p(Qt = zt ) as p(xt ) and p(zt ), respectively. Then,p(x, z) is the joint distribution where x is a time series and z
is the sequence of states the HMM stays at while generating
x. In HMMs, this joint distribution is factored as follows:
(∏
) (∏
)
T −1
T
p(x, z) = p(z1 )t=1 p(zt+1 |zt )t=1 p(xt |zt ) .
(1)

Here p(z1 ) follows a categorical distribution {π1 , π2 , . . . ,πK }, where πk is the probability of the state sk being chosen
as the initial state. Similarly, when zt = sk , p(zt+1 |zt ) follows a categorical distribution {πk1 , πk2 , . . . , πkK }, where
πkk0 is the probability of the state sk0 being chosen as
the next state. Besides, when zt = sk , p(xt |zt ) follows aGaussian distribution N (µk , σk2 ). Then, πk , πkk0 , µk andσk (k, k 0 ∈ [1, K]) are considered as parameters of the joint
distribution, and jointly denoted by θ.
In Bayesian approaches, we consider the extended joint
distribution p(x, z, θ), which is the product of the prior distribution p(θ) and the likelihood p(x, z|θ). For mathematicalconvenience, we use conjugate priors to define the priordistribution. That is,∏Kp(θ) = p(π1 , . . . , πK ) k=1 p(πk1 , . . . , πkK )p(µk )p(λk )(2)where we introduce the precision λk = 1/σk2 and assumep(π1 , . . . , πK )p(πk1 , . . . , πkK )p(µk )p(λk )∼ D(π1 , . . . , πK |α1 , . . . , αK )(3)
∼ D(πk1 , . . . , πkK |αk1 , . . . , αkK )(4)
∼ N (µk |mk , (τ λk )−1 )(5)∼ G(λk |a, b).

(6)Here D(ρ1 , . . . , ρK |α1 , . . . , αK ) and G(λ|a, b) denote theDirichlet distribution and the Gamma distribution, respectively (see the appendix for their definitions). Here αk , αkk0 ,mk , τ , a and b (k, k 0 ∈ [1, K]) are the parameters of the priordistribution, called hyperparameters and jointly denoted byφ.
B. VB learning for continuous HMMs
In variational Bayes (VB), we first consider the logdef
arithm∑ ∫of the marginal likelihood L(x) = log p(x) =log z Θ p(x, z, θ)dθ. Here we fix the hyperparameters
φ and omit them for the moment. Now x is said to beincomplete since we cannot know the state sequence z
from x, and then, some approximation method is required
to compute L(x). We consider an approximation of L(x)via VB. To be more specific, we introduce the variationaldef ∑ ∫
p(x,z,θ)free energy F [q] =z Θ q(z, θ) log q(z,θ) dθ, which isa functional of the test distribution q(z, θ). It can be shownthat F [q] is a lower bound of L(x), and hence the maximized
free energy F [q ∗ ] can be seen as a good approximation ofL(x). Then, as we will see in the next section, the maximizerq ∗ is used to find a desirable discretization.
In VB, we assume q(z, θ) ≈ q(z)q(θ), and obtain ageneric form of variational Bayesian EM (VB-EM) algorithm as an iterative procedure consisting of the following
two updating rules:(∫)
VB-E step : q(z) ∝ exp Θ q(θ) log p(x, z|θ)dθ(7)∑
VB-M step : q(θ) ∝ p(θ) exp ( z q(z) log p(x, z|θ)) (.8)

Now we can derive the algorithm specific to continuous
HMMs by substituting the distribution form of continuous
HMMs to the generic VB-EM procedure above. As a result,
we obtain the adjusted hyperparameters φ∗ that specify q ∗ :αk , αkk0 , mk , τ k , ak and bk . The derived VB-EM algorithm
is presented in the appendix.
To determine the number K ∗ of discrete levels, we
further extend the joint distribution to p(x, z, θ, K) where
K indicate the number of discrete levels. Then we find
K = K ∗ that maximizes the probability p(K|x). Assumingp(K) is uniform, and having p(K|x) = p(x|K)p(K)/p(x)from the Bayes’ theorem, this maximization is equivalent to
the maximization of the logarithm∑ ∫of the marginal likelihood
L(x) = log p(x|K) = log z Θ p(x, z, θ|K)dθ. Again,
L(x) is approximated by the variational free energy, whichis thus used as the score on the plausibility of the numberof discrete levels.
C. Finding the discretized time series
From the property L(x) − F [q] = KL(q(z, θ), p(z, θ|x)),where KL denotes the Kullback-Leibler divergence, finding
the maximizer q ∗ of F [q] leads to a good approximationof p(z, θ|x), the posterior distribution of hidden state sequence and the parameters. Averaging by this approximatedposterior distribution, we compute the predictive
distribution∫of the (hidden) state sequence p(z|x) ∝ Θ q ∗ (θ)p(x, z|θ)dθ
for the input series x. Finally we get the most probable state
sequence (the discretized time series) z ∗ = argmaxz p(z|x).
Note here that p(z|x) is analytically obtained but cannot
be computed in a dynamic programming fashion [6]. To
remedy this problem, we take a heuristic approach known
as reranking [16]. That is, we first find top-n ranked state
sequences {z1 , . . . , zn } by the n-best Viterbi algorithm (weused n = 10) working on an HMM whose initial state distribution, transition distributions, and emission distributions
are independently averaged by q ∗ (θ)3 . Then, we computep(zi |x) exactly for each i ∈ [1, n], and finally obtain z ∗ asargmaxzi :i∈[1,n] p(zi |x).
D. Hyperparameter settings for hybridization
As in Sections II-D and II-E, in discretization, it is crucial
to find the positions (µk ) and the shapes (σk2 ) of the emission
Gaussian distribution which capture well the input series
x. The posterior distribution of µk (i.e. q ∗ (µk )) follows aGaussian distribution whose mean is the adjusted mk , andhence mk eventually affects the quality of discretization. In
the VB-M step, we obtain mk by
mk := (τ mk + T k xk )/(τ + T k ),(9)

3 Consequently, the initial state distribution, the transition distributions,and the emission distributions
respectively follow a categorical distribu∑∗∗
α }, categorical distributions {πkktion {π∑0 =k = αk /` ` k∈[1,K]αkk0 / ` αk` }k0 ∈[1,K] and student’s t-distributions (details omitted).

where mk and τ are the hyperparameters, xk and T k can
be interpreted as the mean of real values emitted whilestaying at the state sk , and the expected counts of staying
at sk , respectively. We can see from Eq. 9 that mk is anweighted sum of mk and xk , and that mk can be controlled
by hyperparameters mk and τ , which are the mean from ourprior knowledge and its weight, respectively.
In combining HMMs with Persist, we determine mk using
the breakpoints βk obtained as a byproduct by Persist,and the 5% and 95% percentile values of the input x.
0
Specifically, we compute mk := (βk−1+ βk0 )/2 where0= ζ95 (x).
βk0 = βk for k ∈ [1, K − 1], β00 = ζ5 (x) and βK
The weight τ is set to balance the prior knowledge andthe input series. The hyperparameters αk and αkk0 play thesame role, and we fixed a ≈ 21 and b = 12 σ̂ 2 (x) in our
experiments. Note that we can also combine HMMs with
SAX using the breakpoints βk , which are SAX’s byproduct.
The hybridization above is surely simple, but is flexible since
the hyperparameters are the only connection point between
the discretizers to be combined.
IV. EXPERIMENTS
A. The enduring-state dataset revisited
To test our hybrid method, we conducted three experiments with artificial and real datasets. First, we test our
method with the enduring-state dataset, described in Section II-E. The experimental settings are the same except that
we set αk = 1 and αkk0 = 1, and chose the hyperparameterτ from {0.5, 1, 5, 10, 20, 50, 70, 100} that works best foreach pair of the number K of discrete levels and the ratio R
of outliers. The results with five levels are shown in Fig. 2 inwhich “HMM+P” and “HMM+S” respectively indicate theHMM combined with Persist and the one combined with
SAX. As expected, these hybrids outperform the original
single methods under both accuracy and NMI. For the othercases under accuracy, by Wilcoxon’s rank sum test with thesignificance level 0.01, “HMM+P” is shown to be better
than the original Persist, except several cases with a large
number of levels and many outliers4 . This superiority was
also confirmed under NMI for all cases.
B. Time series classification
The task in the second experiment is supervised timeseries classification, where we classify a whole time series
into one of the predefined classes. The 1-nearest neighbor(1-NN) classifiers are√∑often used with the Euclidean disdef
0 2tance ∆(x, x0 ) =t∈[1,T ] (xt − xt ) as a dissimilarity
measure between two raw time series x and x0 . In thisexperiment, to see how the discretization methods affect theclassification performance, we replace ∆ with its discretized
4 More precisely, in the cases (K, R) = (6, 7%), (6, 8%), (6, 10%),(7, 4%), (7, 7%), (7, 8%), (7, 9%) and (7, 10%), there is no significant
difference.

Table I
P REDICTIVE PERFORMANCE IN TIME SERIES CLASSIFICATION .
Error rates (%) for the CC dataset
K
2
3
4
5
6
7
8

SAX
*24.61±0.94*4.73±0.55*5.08±0.58*3.03±0.42*2.52±0.38*1.57±0.30*1.40±0.30

K
2
3
4
5
6
7
8

SAX
*23.64±0.82
5.51±0.48
4.41±0.48
2.84±0.38
1.85±0.33
1.40±0.24
1.43±0.27

Persist24.82±1.01
18.58±0.97
17.17±1.07
14.03±0.84
10.67±0.79
11.73±0.84
11.40±0.71

HMM38.15±1.71
8.70±0.90
16.48±1.67
19.58±1.63
15.68±2.1514.65±2.65
12.92±2.48

HMM+SAX36.75±0.72
8.00±0.58
10.38±0.80
25.08±0.70
25.47±0.73
6.23±0.62
4.40±0.62

HMM+Persist
36.70±0.70
8.55±0.60
11.65±1.1614.95±1.52
23.33±0.60
11.05±1.316.47±0.92

Error rates (%) for the CBF dataset

Figure 2. Predictive performance of Persist, HMMs and two hybrids inthe enduring-state dataset with five discrete levels. The plots for Persist andHMMs are the same as those in Fig. 1.

√∑
def
0 2
˜
version ∆(x,x0 ) =
t∈[1,T ] (kt − kt ) , where (k1 , . . . ,kT ) and (k10 , . . . , kT0 ) are the discretized series of x and x0 ,respectively5 .
We worked on two widely-used artificial datasets called
the control chart (CC) dataset and the cylinder-bell-funnel(CBF) dataset6 . The CC dataset contains 600 time series oflength 60 and has six classes, and the CBF dataset contains
798 time series of length 128 and has three classes. Alltime series in each dataset are standardized in advance.
We conducted 10 times 10-fold stratified cross-validation.
In each fold, we first train a discretizer with the training
time series, and then discretize all time series by the trained
discretizer7 . Finally we apply the 1-NN classifier in a usual
manner. In SAX, we chose the frame width w from {1,2, 3, 5, 10, 15, 20, 25, 30} that works best. Also wetrained HMMs under αk = 1, αkk0 = 1 and τ chosenfrom {1, 10, 100, 1000} with 10 reinitializations. Table I
shows the error rates (%) with 95% confidence interval under
Student’s t-distribution with various numbers K of discretelevels. In each row, the smallest error rate is marked with the5∆
˜ indicates that the penalty cost of mismatching between two discretelevels is the square of their difference. In SAX, a dissimilarity measure,called MINDIST, is proposed [9]. However, we do not present the resultswith MINDIST here, since the classifier based on MINDIST worked poorly
with fewer discrete levels (e.g. with three levels).
6 The CC dataset is available from the UCI Machine Learning Repository(http://archive.ics.uci.edu/ml/), and the CBF dataset was generated by theauthors based on the description in [11].
7 This procedure is necessary for data-adaptive discretizers such as Persist
and HMMs. We extended Persist and HMMs in a straightforward way tohandle more than one time series.

Persist25.99±0.77
12.29±0.80
4.01±1.28
1.05±0.25
2.83±0.45
1.90±0.30
1.99±0.44

HMM
27.55±1.03
2.43±0.33*0.76±0.19
1.23±0.27
0.80±0.21
1.20±0.25*0.83±0.21

HMM+SAX27.50±1.02
2.29±0.31*0.76±0.19*0.98±0.24*0.74±0.19
0.90±0.26
1.01±0.25

HMM+Persist
25.81±1.01*1.74±0.27
1.24±0.24
0.99±0.22
0.83±0.20*0.82±0.20
0.88±0.19‘*’ symbol. Without discretization (i.e. using ∆), the error
rates (%) are 7.90±0.62 for CC and 1.85±0.27 for CBF.
In the results for the CC dataset, SAX outperforms theother methods including the method without discretization.
This seems to be due to the effect of smoothing by PAA [1],[9], and more importantly to the assumption made by themethods other than SAX — time series have flat portions
which are crucial to the goal of the application. Indeed,
the time series with long increasing/decreasing trend incur more errors to these methods. On the other hand, itshould be remarked that the two hybrids “HMM+SAX” and“HMM+Persist” surely improve the performance of 1-NN
over the original HMMs in most cases. Also with three
discrete levels, the performance with HMM-based methods
is comparable to the one without discretization.
For the CBF dataset, the performance with two hybrids isbetter than the one without discretization, and than the one
with SAX or Persist especially when we have fewer discretelevels. This result is important since from the viewpoint ofcomprehensibility, the number of discrete levels should be
preferred to be small. Besides, we observed that, in bothdatasets, a large τ is sometimes needed for hybrids, since
the training set contains numerous measurements (e.g. 600×
0.9 × 60 = 32, 400 in the case of CC).
C. The muscle dataset
In the last experiment, following [8] again, we discretize
the time series on the muscle activation of a professional
inline speed skater [10]8 , focusing on determining the number of discrete levels. In this dataset, the measurement is
the muscle activation calculated from the original EMG
8 The dataset is included in the package of Persist’s MATLAB implementation (http://www.mybytes.de/persist.php).

-2200τ = 500 without PAA-30000PAA
weight = 500
weight = 1-40000-50000

-2400

-2600

τ = 1 with PAA under w = 20-60000-2800-70000
2

3

4

5

6

7

-3000
8

Number of discrete levels

Figure 3. (left-above) The original and discretized time series. (left-below) The time series compressed by PAA and discretized. (right) The variational
free energy with different number of discrete levels.

(Electromyography) as the logarithm of the energy, andnearly 30,000 measurements are recorded. The plots with athin line in Fig. 3 left-above shows the original time series.
Since the time series contains numerous measurements,
to determine the number of discrete levels, we should set
relatively large numbers to τ , αk and αkk0 , the weights forprior knowledge (Section III-D). Otherwise, the estimatednumber K ∗ of discrete levels turns to be counterintuitively
large. This overfitting problem is illustrated in Fig. 3 right,
where the x-axis indicates the number K of discrete levels
and the y-axis indicates the variational free energy (VFE)with K levels. The plots “weight = 1” and “weight = 500”
show the VFEs under τ = 1 and τ = 500, respectively9 .
Also, we can have a similar effect by compressing the timeseries with PAA. The plots “PAA” indicate the VFEs under
the frame width w = 20 and τ = 1.
The plots with thick lines in Fig. 3 left indicate thediscretized sequence in the case of τ = 500 without PAA(above), and in the case of τ = 1 with PAA under w = 20
(below). In both cases, K ∗ is estimated as four, but one ofthe discrete levels (whose mean is at −17.7) appears only
rarely. This result coincides with the fact that the expert
prefers K ∗ to be three [8], and the discretized sequencewith PAA clearly shows the high activation at the end ofeach cycle, which means the last kick to the ground to move
forward. Note here that this result has been obtained using
Persist [8], but we reached a similar result by a more generalapproach via variational Bayes.

of data points in a two dimensional (measurement-time)
space. As examined by experiments, VB provides a principled and flexible way to build a robust discretizer and to
determine the number of discrete levels, which would finallylead to a discretization meaningful for symbolic data mining
or other AI systems. Also the hybrid discretizer we built
is data-adaptive, and hence is expected to produce a morecomprehensible discretization. To the best of our knowledge,
the proposed method is the first attempt to introduce VB for
discretization tasks.
There are some remaining works. We need to apply the
proposed method to other datasets in time series classification. Computational issues are important as well. In
particular, the computation time of HMM-related algorithms
is quadratic to the number of discrete levels. This problem isremedied by a compression method like PAA, and as mentioned in Section IV-B, fewer discrete levels are preferredfrom the viewpoint of comprehensibility. Besides, there arecases where we wish to discretize multivariate time series
with common discrete levels. Synnaeve et al. [4] realized itby extending HMMs with the notion of parameter tying [7].
ACKNOWLEDGMENT
This work is supported in part by Grant-in-Aid for Scientific Research (No. 20240016) from Ministry of Education,
Culture, Sports, Science and Technology of Japan.
R EFERENCES

V. CONCLUSION[1] C. Daw, C. Finney, and E. Tracy, “A review of symbolic analysis of experimental data,” Review of Scientific Instruments,
vol. 74, no. 2, pp. 915–930, 2003.

We proposed a hybrid method for unsupervised timeseries discretization based on variational Bayes (VB). By hybridization, we can benefit from heterogeneous discretizers
which have their own assumption/model on the distribution[2] J. Dougherty, R. Kohavi, and M. Sahami, “Supervised andunsupervised discretization of continuous features,” in Proc.
of the 12th Int’l Conf. on Machine Learning (ICML-95), 1995,
pp. 194–202.

9 In this experiment, we always set (τ + 1) to the other weight hyperparameters αk and αkk0 , and tried 500 reinitializations in the VB-EM
algorithm.

[3] R. King, S. Garrett, and G. Coghill, “On the use of qualitative reasoning to simulate and identify metabolic pathways,”
Bioinformatics, vol. 21, no. 9, pp. 2017–2026, 2005.

[4] G. Synnaeve, A. Doncescu, and K. Inoue, “Kinetic models
for logic-based hypothesis finding in metabolic pathways,” in19th Int’l Conf. on Inductive Logic Programming (ILP-09),
2009.
[5] H. Attias, “A variational Baysian framework for graphicalmodels,” in Advances in Neural Information Processing Systems 12 (NIPS-99), 1999, pp. 209–215.
[6] M. J. Beal, “Variational algorithms for approximate Bayesianinference,” Ph.D. dissertation, University College London,
2003.
[7] L. Rabiner and B.-H. Juang, Fandamentals of Speech Recognition. Prentice Hall, 1993.
[8] F. Mörchen and A. Ultsch, “Optimizing time series discretization for knowledge discovery,” in Proc. of the 11th
ACM SIGKDD Int’l Conf. on Knowledge Discovery and Data
Mining (KDD-05), 2005, pp. 660–665.
[9] J. Lin, E. Keogh, L. Wei, and S. Lonardi, “Experiencing SAX:
a novel symbolic representation of time series,” Data Mining
and Knowledge Discovery, vol. 15, no. 2, pp. 107–144, 2007.
[10] F. Mörchen, A. Ultsch, and O. Hoos, “Extracting interpretablemuscle activation patterns with time series knowledge mining,” Int’l J. of Knowledge-based and Intelligence Engineering Systems, vol. 9, no. 3, pp. 197–208, 2005.
[11] P. Geurts, “Pattern extraction for time series classification,” inProc. of the 5th European Conf. on Principles and Practiceof Knowledge Discovery in Databases (PKDD-01), 2001, pp.
115–127.
[12] S. Watanabe, Y. Minami, A. Nakamura, and N. Ueda, “Application of variational Bayesian approach to speech recognition
recognition,” in Advances in Neural Information Processing
Systems 15 (NIPS-2002), 2002, pp. 1237–1244.
[13] C. D. Manning, P. Raghavan, and H. Schütze, Introduction to
Information Retrieval. Cambridge University Press, 2008.
[14] M. O. Abu-Shawiesh, F. M. Al-Athari, and H. F. Kittani,“Confidence interval for the mean of a contaminated normal
distribution,” J. of Applied Science, vol. 9, no. 15, pp. 2835–
2840, 2009.
[15] S. Ji, B. Krishnapuram, and L. Carin, “Variational Bayesfor continuous hidden Markov models and its applicationto active learning,” IEEE Trans. on Pattern Analysis andMachine Intelligence, vol. 28, no. 4, pp. 522–532, 2006.
[16] M. Collins, “Discriminative reranking for natural language
parsing,” in Proc. of the 17th Int’l Conf. on Machine Learning(ICML-2000), 2000, pp. 175–182.

A PPENDIXThe Dirichlet and the gamma distributions are defined as:∑KΓ(αk ) ∏K
def
αk −1
k=1
∏
D(ρ1 , . . . , ρK |α1 , . . . , αK ) =K
k=1 ρkk=1

G(λ|a, b)def=Γ(αK )1a a−1
exp(−bλ),Γ(a) b λwhere Γ is the gamma function. In the VB-EM algorithm, we
initialize αk = αk +ε, αkk0 = αkk0 +ε, mk := mk +ε, ak :=a, bk := b and τ k := τ for each k ∈ [1, K]. Here each of
ε’s is a distinct small random noise, and the hyperparameters
αk , αkk0 , mk , τ , a and b are specified in advance. Then, inthe VB-E step, we compute the expectations U k , T kk0 , xk
and S k by the following procedure:
(∏
) (∏
)
T −1
T
q(z) := q(z1 )t=1 q(zt+1 |zt )t=1 q(xt |zt )∑q(sk ) := exp (Ψ(αk ) − Ψ ( k αk ))∑q(sk0 |sk ) := exp (Ψ(αkk0 ) − Ψ ( k0 αkk0 ))( 1())
q(xt |sk ) := exp − 2 log 2π + τ1k + log bk − Ψ(ak ) ·)( ak(xt − mk )2
exp − 2bk∑
U k :=q(z)z:z =s
∑ 1 k
T kk0 :=z q(z) {t ∈ [2, T ] | zt−1 = sk , zt = sk0 }∑∑1
xk := Tz q(z)
t:t∈[1,T ],zt =sk xt
k
∑
∑
1
2
S k := Tz q(z)
t:t∈[1,T ],zt =sk (xt − xk ) ,kdef ∑where T k =k0 T kk0 and Ψ is the digamma function:
def dΨ(x) = dx log Γ(x). The VB-E step can be computedby the forward-backward algorithm on an HMM whose
initial state distribution, transition distributions and emission
distributions are q(z1 ), q(zt+1 |zt ) and q(xt |zt ), respectively.
In the VB-M step, on the other hand, we update thehyperparameters as below.

αk
mk
αkk0
ak
τk
bk

:= αk + U k:= (τ mk + T k xk )/(τ + T k ):= αkk0 + T kk0
1:= a + T k
2:= τ + T k
)
(
1
τTk
2
:= b +(mk − xk ) + S k
2 τ + Tk

We iterate the VB-E step and the VB-M step alternately until
the convergence of the variational free energy derived as:
F [q] :=log+∑

z q(z) +

∑

∑

∑

Γ(
Γ(αk )∑k αk )
k log Γ(αk ) + log Γ(
αk )

∑

k

∑
Γ(
Γ(αkk0 )∑k0 αkk0 )k,k0 log Γ(αkk0 ) +k log Γ(αkk0 )

∑

0

k
∑
U
(Ψ(α
)
−
Ψ
(
α
))
k
k
k
k
∑k
∑
+ k,k0 T kk0 (Ψ(αkk0 ) − Ψ ( ` αk` ))
)
∑ (
ba
k)
+ k 21 log ττ + log Γ(aΓ(ak ) + log bakk∑+ k ab k ((bk − b) − 12 τ (mk − mk )2 )∑ k
+ 21 k T k ( τ1k + log bk − Ψ(ak )).

+