Visualization of Anomalies
using Mixture Models
Tomoharu Iwata and Kazumi Saito
NTT Communication Science Laboratories, NTT Corporation
Hikaridai 2–4, Seika-cho, Soraku-gun, Kyoto, 619–0237, Japan
Abstract. Anomaly detection is important to learn from major past
events and to prepare for future crises. We propose a new anomaly detection method that visualizes multivariate data in a 2- or 3-dimensional
space based on the probability of belonging to a mixture component
and the probability of not belonging to any components. It helps tovisually understand not only the magnitude of anomalies but also therelationships among anomalous and normal samples. This may provide
new knowledge in the data, since we can see it from a different viewpoint.
We show the validity of the proposed method by using both an artificial
and an economic time series.

1

Introduction

In recent years, a huge amount of data has been electronically accumulated. It
is important to detect anomalies in data to learn from major past events and toprepare for future crises. Anomaly detection techniques are used for various such
purposes as the detection of unauthorized computer use by examining the user
log data [8] or the detection of anthrax outbreaks by tracking over-the-counter
medication sales [7].
Many anomaly detection methods have been proposed, mainly in the field of
outlier detection [2]. However, most resort to plotting anomaly scores in the data[19]. Such methods can determine to what degree a sample in data is anomalous, however, they cannot grasp relationships among anomalous and normal
samples. Since anomalies happen for various reasons, it is important to know
their relationships to understand their causes.
To overcome this problem, we propose a new method of the visualizationof multivariable data for anomaly detection that uses a mixture model. The
method visualizes data in a 2- or 3-dimensional space based on the probabilityof belonging to a mixture component and the probability of not belonging to
any components. The results of visualization provide not only the magnitude ofanomalies but also the relationship among data. This can be used to find new
knowledge in the data, since we can see it from a different viewpoint.
The remainder of this paper is organized as follows. In the next section, we
present our new method, and in Section 3, we explain how to apply our method
to a time series. In Section 4, we show the validity of our method by using
both an artificial and an economic time series. In the last section, we present
concluding remarks and future works.

2

Proposed Method

In this section we present a new method to visualize data for anomaly detection.
The procedure follows:
1. Given particular data, build a mixture model.
2. Estimate the probability of belonging to each mixture component.
3. Calculate the average 2-sigma value α̃, which is used to detect anomalies.
4. Create an augmented probability vector z by combining the probabilities
estimated in Step 2 and α̃.
5. Visualize z by using a dimension reduction method.
2.1

Mixture Model

Let D = {D(1), . . . , D(N )} be a given data set, where N is the number of data
samples. A mixture model assumes that each sample D(i) is generated by a finite
mixture of underlying probability distributions:p(D(i)|Θ) =KX

P (k)p(D(i)|k, Θ),

(1)

k=1

PK
where k=1 P (k) = 1, P (k) ≥ 0, K is the number of components, and Θ is avector constructed by all parameters. We estimate parameters Θ by maximum
likelihood estimation:Θ̂ = arg max p(D|Θ) = arg max
Θ

Θ

N
Y

p(D(i)|k, Θ),

(2)

i=1

where Θ̂ is the estimated parameters by maximum likelihood method.
We can apply any kind of mixture component, such as a normal distribution
model, a multinomial distribution model, or an autoregressive model. It is important to select a component that can appropriately model the given data andwhose parameters can be easily estimated.
2.2

Augmented Probability Vector

Next, we explain how to obtain an augmented probability vector, which is the
combination of the probability of belonging to each component and not belonging
to any components. We assume that a sample that does not belong to any
components is an anomaly.
We consider the joint probability density of k-th component and i-th sample,
P (k, D(i)), as the probability of belonging to component k. We define the probability vector of i-th sample, z∗ (i), by a vector where the element is the joint
probability density as follows:z∗ (i) = (p(k = 1, D(i)), . . . , p(k = K, D(i)))= (P (k = 1)p(D(i)|k = 1), . . . , P (k = K)p(D(i)|k = K)).

(3)z∗ (i) represents the probabilities of belonging to each component of the i-th sample. Anomalies have low probabilities of belonging to any components. Therefore,
by augmenting a small value α̃ to z∗ and normalizing it so that their sum is 1,
the vector represents the probabilities of belonging to each component (by 1st ∼K-th elements) and not belonging to any components (by (K + 1)-th element).
Then, the vector z(i) is:z(i) ∝ (P (k = 1)p(D(i)|k = 1), . . . , P (k = K)p(D(i)|k = K), α̃),(4)PK+1where k=1 z(i)k = 1. We call this an augmented probability vector. We determine α̃ based on the average 2-sigma value, that is, the average of probability
densities on 2-sigma of all components:
α̃ =

K
X

k=1

P (k)p(µk ± 2σk |k),(5)

where µk and σk stand for the mean and standard deviation of conditional
probability P (D(i)|k). The probability that a sample exists in the 2-sigma area
is 0.9545 in the case of a normal distribution. Therefore, most samples havehigher joint probability densities than the average 2-sigma value, and so we
consider them normal. Since a few samples have lower probabilities of belonging
to any components than an average 2-sigma value, we consider them anomalies.
Here we use an average 2-sigma value, however, an average 3-sigma value orothers can be defined in the same way. To reduce false alarm, we could adopt anaverage 3-sigma value. We can choose α̃ in accordance with the application.
2.3

Visualization

Last, we explain how to visualize the augmented probability vectors and how to
study the results of visualization. An augmented probability vector is a vectorof K + 1 dimension. We need to reduce the dimension of this vector to 2 or 3 tovisualize it. There are many dimensionality reduction methods, such as multidimensional scaling [15], spring model[9], and cross-entropy directed embedding(CE) [18]. We can apply any method to reduce the dimension of augmented
probability vectors, in this paper we used the CE method.
Cross-entropy embedding (CE) is a nonlinear embedding method that embeds samples to minimize energy function based on cross-entropy between similarities among original data and proximity among embedded data. Let si,j be
the similarity between z(i) and z(i), calculated by cosine similarity:si,j =

z(i)0 z(j).
|z(i)||z(j)|(6)

Let ri be the coordinate of i-th sample in the embedding space, and the proximity
between ri and ri can be defined as follows:ρ(ri , rj ) ∝ exp(−1
k ri − rj k2 ).
2(7)

Then, the cross-entropy between similarity si,j and proximity ρ(ri , rj ) is:
Ei,j = si,j log(p(ri , rj )) + (1 − si,j )log(1 − ρ(ri , rj )).

(8)Ei,j attains its minimum value when ρ(ri , rj ) = si,j , that is, when similar samples are closely located and dissimilar samples are located far away. An objectivefunction with a regularization term to be minimized is:
J=

N
−1
X

N
X

i=1 j=i+1

N

Ei,j −µX
k ri k2 ,
2 i=1

(9)

where µ is a predetermined constant. By virtue of CE, samples that have similaraugmented probability vectors are closely located.
As a result of visualization, samples that are classified into the same component are closely located and form clusters. Usually, there are a few anomalies.
Therefore, anomalies are located far from the normal samples that form clusters.
The results tell us to which component the anomalies belong, and elucidate therelationship among data. In this way, the visualization of augmented probability
vectors makes it possible to understand not only the magnitude of anomalies
but also their characteristics and their relationships among data.
It is natural to characterize each sample D(i) by a vector z(i) by focusing on
fitting into each component model in which probabilistic structures of variables
are specified. Another commonly used visualization method is based on valuesof the sample itself, using a dimensionality reduction method such as multidimensional scaling [12]. However, we believe that it is difficult to find someessential characteristics by directly applying such a visualization method because
it does not use the probabilistic structure of variables.

3

Application to a Time Series

Here, we explain in detail how to apply our method to a time series. We used avector autoregressive (AR) mixture model [17]. AR models are widely used fortime series analysis [10].
3.1

AR Mixture Model

Let Y = (y(1), . . . , y(T )) be a given time series, where T is the number oftime points and y(t) is a d-dimensional vector of variables at time t. Also, letx(t) = (y0 (t − 1), . . . , y0 (t − τ ))0 be an input vector at time t to the AR model,
where τ is an order of the AR model and 0 means the transpose of the vector.
An AR mixture model assumes that y(t) is generated from the K mixtureof the AR models:p(y(t)|x(t)) =

K
X

k=1

P (k)p(y(t)|x(t), k),(10)p(y(t)|x(t), k) =(y(t)−Ak x(t))0 (y(t)−Ak x(t))
1
exp(−
),
2σk2
(2πσk2 )d/2

(11)

where Ak is the AR model parameters of component k, and σk2 is the varianceof component k. Here, we assume that the variance of all variables are the sameand independent in each component.
3.2

Estimation of Parameters

The estimation procedure for the parameters is as follows:
Step1 Let K ∗ = 1, and estimate A1 by minimizing the square error:
E1 =

T
X

(y(t) − A1 x(t))0 (y(t) − A1 x(t)).

t=τ +1

Step2 Select s (s ≤ K ∗ ), and then split the parameters:
Ak+1 = As − ∆A.

As = As + ∆A,

Step3 Let K ∗ = K ∗ + 1.
Step4 Estimate the parameters by maximizing the log likelihood:
LK ∗ =

T
X

∗

log

t=τ +1KX

P (k)p(y(t)|x(t), k).

k=1

Step5 If K ∗ < K, then return to Step2.

In Step2, we have three choices: the component s that has the highest P (s), thecomponent s with the highest σs , or randomly. In our experiments, we randomly
selected s for simplicity. In Step4, it is impossible to estimate parameters by
maximizing LK ∗ analytically, so we approximately estimate parameters using
an EM algorithm [6]. The E- and M-Steps in the EM algorithm are:
– E-Step:
P (k)p(y(t)|x(t), k)
P (k|x(t), y(t)) = PK ∗,∗∗k∗ =1 P (k )p(y(t)|x(t), k )– M-Step:

P (k) =

σk2 =

Ak =

Ã

T
X

t=τ +1PT

t=τ +1(12)TX
1
P (k|x(t), y(t)),
T − τ t=τ +1

P (k|x(t), y(t))(y(t) − Ak x(t))0 (y(t) − Ak x(t))
dPT

t=τ +1P (k|x(t), y(t))y(t)x(t)

0

(13)

,

(14)

P (k|x(t), y(t))

!Ã

T
X

t=τ +10

P (k|x(t), y(t))x(t)x(t)!−1

By iterating the above E- and M-steps, we obtain a local optimal solution.

.

(15)

3.3

Average 2-sigma value and Augmented Probability Vector

The average 2-sigma value in an AR mixture model α̃M AR is the average ofprobability densities on 2 sigma of K AR components as follows:K
X

1α̃M AR =
P (k)
exp(−
(2πσk2 )d/2
k=1
=

K
X

P (k)

k=1

Pd

2
j=1 (2σk )
)
2σk2

1
exp(−2d).
(2πσk2 )d/2

(16)

Then the augmented probability vector z(t) in the AR mixture model is:z(t) = (P (k = 1)p(y(t)|x(t), k = 1), . . . , P (k = K)p(y(t)|x(t), k = K), α̃M AR ).
(17)
Clearly, α̃M AR is independent of time t.

4

Evaluation of Our Proposed Method

4.1

Artificial Time Series

We evaluated our method by applying it to an artificial time series with 400
sampling points. Until time 200, the time series was generated using A1 , which
was a two-variable one-order AR model parameter, and from time 201, generatedusing A2 , which was another AR model parameter. That is,x(t) = Ai y(t) + U,½
1 if 1 ≤ t ≤ 200,i=
2 if 201 ≤ t ≤ 400,

(18)
(19)

where U is 2-dimensional white noise with covariance matrix Σ. Then to artificially introduce an anomaly, we added 5 to the 1st variable at time 100. We
used following parameters:
µ
¶
µ
¶
µ
¶
0.2 0.3−0.5 1.3
1 0.2
A1 =, A2 =, Σ=,(20)−0.6 1.1
0.3 0.3
0.2 1
where Σ is time invariant. Figure 1 is an example of this time series.
Figure 2 shows the 3-dimensional visualization results of this time series using
our method with τ = 1 and K = 1, 2. In the case of K = 1 (Figure 2(a)), most
samples are clustered in the lower left. The anomalous sample at time 100 is
not located furthest from the cluster, that is, it is not visualized as the mostanomalous. In the case of K = 2 (Figure 2(b)), there are 2 clusters in the lower
left and the lower right, in which most samples are located. Unlike K = 1, theanomaly is located furthest from the others, so it is visualized as the anomaly.
Our method can detect anomalies with an appropriate K.

8

a

6
4
2
0
-2
-4
-6

0

50

100

150

200

250

300

350

400

350

400

A2

A1
8
6
4
2
0
-2
-4
-6
-8

0

50

100

150

200

250

300

Fig. 1. Two-variable artificial time series: (upper) the 1st variable, (lower) the 2nd
variable. ’a’ in the upper graph represents the time point of anomaly that we added.

Figure 2(b) also shows the efficiency of the mixture model and the characteristics of the given samples. In the lower right, mainly samples until 200(blackish points) are located, and in the lower left, mainly samples after 200(whitish points) are located. These result indicate that the estimation was successful on some level. However, there are many points in the middle of the 2clusters, implying that it is difficult for the model to classify some samples.
4.2

Japanese Economic Time Series

Next, we evaluated our method using real data from a monthly economic timeseries in Japan covering the period 1983-2003. There is a total of 240 points,
and the six measured variables are: monetary base, national bond interest rate,
wholesale price index, industrial produce index, machinery orders, and yen to
dollar exchange rate (Figure 3). Monetary base represents the average outstanding balance deposited in the Bank of Japan. The wholesale price index represents
an average of the selling prices of domestic products, excluding services. The industrial produce index measures the goods produced by the economy. The raw
data are non-stationary and transformed to a stationary series through several
operations, including adjustments for seasonality and trends [13].
We visualized this time series with τ = 1 and K = 1 ∼ 4 (Figure 4). One of
the biggest economic events during this period was the Russian economic crisis ofAugust 1998, which caused a large hedge fund company to fail, and the Bank ofJapan rapidly raised interest rates in September 1998. At K = 1 ∼ 4, the month(a) K=1

0

50

100

(b) K=2

150

200

250

300

350

400

(c)
Fig. 2. Results of 3D visualization of anomaly in artificial time series. 4 is artificial
anomaly. (a) K=1, (b) K=2, and (c) the color-map for time of points. The point changes
color from black to white going thorough time.

September 1998 (⊕) is visualized as an anomaly. However, in the case of K = 1,
some samples of other months are also located near September 1998. Also, in thecases of K = 2 and K = 3, other samples are visualized as anomalies. On theother hand, in the case of K = 4, there is only one anomaly, i.e. the sample ofSeptember 1998. If K = 1, the model is a linear model, and thus it is impossible
to appropriately model the data. As K increases, the model becomes flexible,
and then it becomes less likely to show false alarms.

5

Related Work

Anomaly detection methods can be broadly classified into 2 methods: profilingand discriminating [7]. A profiling method builds a model from normal samples,and samples which deviate significantly from the model are deemed anomalies.
A discriminating method builds a model from anomalous samples, and samples
which match the model are deemed anomalies. If a lot of anomalous samples
are given or the model of anomaly is known, a discriminating method is useful.
In general, however, it is difficult to obtain many anomalous samples, and themodel of the anomaly cannot be specified. So, we adopted a profiling methodthat builds a mixture model from data and samples which deviate significantly
from the model based on the 2-sigma value are assumed to be anomalies.

monetary base
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
-0.05

1985

1990

1995

2000

national bond interest rate
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

1985

1990

1995

2000

wholesale price index
0.02
0.01
0
-0.01
-0.02
-0.03
-0.04
-0.05
-0.06
-0.07

1985

1990

1995

2000

industrial produce index
0.15
0.1
0.05
0
-0.05
-0.1
-0.15

1985

1990

1995

2000

machinery orders
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3

1985

1990

1995

2000

yen to dollar exchange rate
0.15
0.1
0.05
0
-0.05
-0.1
-0.15
-0.2

1985

1990

1995

Fig. 3. Japanese economic time series

2000

(a) K=1(b) K=2(c) K=3(d) K=4

Fig. 4. Results of 3D visualization of anomalies in an economic time series. ⊕ is the
sample on 9/1998, °s are samples classified to component 1, 4s to component 2, ×s
to component 3, and ♦s to component 4.

The main characteristics of our method give the relationships among data
as well as the magnitude of anomalies using visualization. Davidson’s method[5] is similar to our method because it visualizes anomalies with the relationships among data. However, he focuses on anomalies in clustering problems, andconsiders samples located near cluster boundaries as anomalies. It is useful to detect vague samples that are difficult to classify when we execute clustering. Our
method can also detect such anomalies, although we do not explicitly consider
them anomalies. We consider samples that do not match any mixture components as anomalies, therefore, our method can apply usual anomaly detection
problems. In this point, our method is very different from Davidson.
We used a mixture of vector AR models for anomaly detection in a timeseries. The components of the mixture model can include more complex models
such as ARMA [3]. Next, P (k) can be estimated from inputs, using a gatednetwork [16] or a threshold [14]. If a mixture AR model is inadequate to modela given time series, our method can be modified easily for other mixture models.
There are also discriminating methods using AR models. Chen et al. considered 4 types of anomalies: (1) additive outliers, (2) innovational outliers, (3) temporary changes, and (4) level shifts. They detected them using 4 types anomaly
models [4]. Experiments using an artificial time series showed that our method
can detect simple additive outliers. We believe that our method can cope withother types of outliers because of the flexible expressive power of mixture models.
However, this claim must be confirmed by further extensive experiments.

6

Conclusions

We proposed a new method to visualize data for anomaly detection using a mixture model. This method visualizes samples based on the probability of belonging
to a component and not belonging to any components. Since our visualization
method is quite different from conventional approaches, it provides a new viewof data and helps to find new knowledge. We also showed the validity of our
method by applying it to an artificial time series and an actual economic timeseries.
Although our results have been encouraging to date, a number of directions
remain in which we must extend our approach before it can become a useful tool
for anomaly detection. First, we need to determine the number of components
K. In the application to an economic time series, we incremented K from K = 1
until only a few anomalies are visualized clearly. It is necessary to determine K
automatically by using a quantitative score such as AIC [1] or MDL [11]. Second,
we must avoid getting trapped in a local minimum. Visualization results of thesame data can be different. In the application we estimated parameters several
times and selected parameters whose likelihood were the highest. However, if
there are many local minima in the data, such a procedure is ineffective. We will
study these problems in the future.

References
1. Akaike, H. : A new look at the statistical model identification. IEEE Transactions
on Automatic Control, 19 (1974) 716–723.
2. Barnett, V., Lewis, T. : Outliers in statistical data. 2nd ed. John Wiley & Sons,
New York (1984).
3. BrockWell, P. J., Davis, R. A. : Introduction to Time Series and Forecasting.
Springer-Verlag (2002).
4. Chen, C., Liu, L. : Joint estimation of model parameters and outlier effects in timeseries. Journal of the American Statistical Association 88 (1993) 284–297.
5. Davidson, I., Ward, M. : A particle visualization framework for clustering andanomaly detection. Proceedings of KDD Workshop on Visual Data Mining (2001).
6. Dempster, A. P., Laird, N. M., Rubin, D. B. : Maximum likelihood from incomplete data via the EM algorithm (with discussion). Journal of the Royal StatisticalSociety, Series B 39 (1977) 1–38.
7. Fawcett, T., Provost, F. : Activity monitoring: Noticing interesting changes in behavior. Proceedings of Fifth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. (1999) 53–62.
8. Gorldenberg, A. : Early statistical detection of anthrax outbreaks by tracking overthe-counter medication sales. Proceedings of the National Academy of Sciences of
the United States of America, 99 (2002) 5237–5240.
9. Kamada, T., Kawai, S. : An algorithm for drawing general undirected graphs.
Information Processing Letters 31 (1989) 7–15.
10. Lutkepohl, H. : Introduction to multiple time series analysis. Spinger-Verlag (1993).
11. Rissanen, J. : Modeling by shortest data description. Automatica 14 (1978) 465–
471.
12. Rorvig, M. : Images of similarity: A visual exploration of optimal similarity metrics
and scaling properties of TREC document sets. Journal of the American Society
for Information Science 50 (1999) 639–651.
13. Sato, S. : Stepwise prediction for economic time series by using vector autoregressive model. Science of Modeling, AIC2003, ISM Report on Research and Education
17 (2003) 225–233.
14. Tong, H., Lim, K. S. : Threshold autoregression, limit cycles and cyclical data.
Journal of the Royal Statistical Society, Series B 42 (1980) 245–292.
15. Torgerson, S. : Theory and methods of scaling. Wiley, New York (1958).
16. Weigend, A., Mangeas, M., Srivastava, A. : Nonlinear gated experts for time series: discovering regimes and avoiding overfitting. International Journal of Neural
Systems 6 (1995) 373–399.
17. Wong, C. S., Li, W. K. : On a mixture auto regressive model. Journal of the Royal
Statistical Society, Series B 62 (2000) 95–115.
18. Yamada, T., Saito, K., Ueda, N. : Cross-entropy directed embedding of network
data. ICML2003 (2003) 832–839.
19. Ye, N., Chen, Q. : An anomaly detection technique based on a chi-square statistic
for detecting intrusions into information systems. Quality and Reliability Engineering International 17 (2001) 105-112.