Li

# Learning from Polar Representation: An Extreme-Adaptive Model for Long-Term
Time Series Forecasting

Yanhong Li  Department of Computer Science and Engineering, Santa Clara
University, Santa Clara, CA, USA  Jack Xu  Santa Clara Valley Water District,
San Jose, CA, USA  David C. Anastasiu Corresponding author: David C.
Anastasiu. E-mail<EMAIL_ADDRESS>Department of Computer Science and
Engineering, Santa Clara University, Santa Clara, CA, USA

###### Abstract

In the hydrology field, time series forecasting is crucial for efficient water
resource management, improving flood and drought control and increasing the
safety and quality of life for the general population. However, predicting
long-term streamflow is a complex task due to the presence of extreme events.
It requires the capture of long-range dependencies and the modeling of rare
but important extreme values. Existing approaches often struggle to tackle
these dual challenges simultaneously. In this paper, we specifically delve
into these issues and propose Distance-weighted Auto-regularized Neural
network (DAN), a novel extreme-adaptive model for long-range forecasting of
stremflow enhanced by polar representation learning. DAN utilizes a distance-
weighted multi-loss mechanism and stackable blocks to dynamically refine
indicator sequences from exogenous data, while also being able to handle uni-
variate time-series by employing Gaussian Mixture probability modeling to
improve robustness to severe events. We also introduce Kruskal-Wallis sampling
and gate control vectors to handle imbalanced extreme data. On four real-life
hydrologic streamflow datasets, we demonstrate that DAN significantly
outperforms both state-of-the-art hydrologic time series prediction methods
and general methods designed for long-term time series prediction.

###### keywords:

multivariate time series analysis | hydrologic prediction | streamflow prediction | representation learning | deep learning models

Figure 1: Streamflow over 3 years.

## Introduction

Time series forecasting has a critical role in diverse domains, enabling
effective resource management and informed policy decisions. However, certain
time series data pose a unique problem because they contain sporadic but
significant extreme events, such as unexpected flash floods or climate change-
induced droughts in the problem of streamflow prediction. The ability to
forecast time series that include these types of extreme occurrences is an
important research direction which has seen much attention in recent years (1,
2, 3, 4, 5, 6).

Traditionally, machine learning and statistics-based models were the basic
foundation for time series prediction (7, 8). However, techniques like
Autoregressive Integrated Moving Average (ARIMA) (9) seem to perform badly
when dealing with large variations in the streamflow values, while other
methods (10, 11, 12) are generally designed for short future horizon
forecasting.

A variety of neural network architectures have been investigated for
hydrologic forecasting, including recurrent neural networks (13, 14), hybrid
networks (15) and graph neural networks (16, 17). Some work employed Extreme
Value Theory (EVT) to enhance the hydrologic time series performance (6, 5).
However, these studies primarily concentrate on short-term forecasting and
their performance on longer time horizons is doubtful. While there has been a
surge in transformer-based forecasting models asserting their high-performance
capabilities for long-horizon time series tasks (18, 19, 20, 21, 22), recent
research has raised questions about their efficacy, indicating that simpler
linear models can outperform them (23, 24). Moreover, imbalanced data or
severe events might hurt all these state-of-art deep learning approaches when
it comes to long-term predictions.

We focus on these challenges and innovatively address them through
representation learning (25, 26, 27), a burgeoning field in unsupervised
learning. Our aim is to extract latent states containing the extreme features
in data for downstream tasks. To achieve this, we explore the potential of
multi-loss functions (28, 29) in shaping our training objective. The main
contributions of this work are as follows:

  * •

We propose a Distance-weighted Auto-regularized Neural network (DAN), which
uses expandable blocks to dynamically facilitate long-term prediction.

  * •

To improve the model’s robustness to severe events, DAN innovatively uses a
distance-weighted multi-loss method to extract the polar representations from
time series simultaneously.

  * •

We introduce a Kruskal-Wallis sampling policy to handle imbalanced extreme
data and gate control vectors to boost the discriminatory capacity of
indicators to accommodate imbalanced data.

  * •

We evaluate DAN and competing methods on four separate datasets and find that
DAN significantly outperforms state-of-the-art baselines. Additionally, we
carry out several ablation studies to comprehend the effects of specific
design decisions.

## Related work

Streamflow forecasting holds a pivotal role in enhancing water resource
allocation, management, flood warning, and mitigation of flood-related
damages. Traditional methods for streamflow forecasting included the
univariate Autoregressive (AR), Moving Average (MA), Simple Exponential
Smoothing (SES), and Extreme Learning Machine (ELM) algorithms, and most
famously the Autoregressive Integrated Moving Average (ARIMA) (9) method and
its several variants. Wang et al. (30) developed a hybrid model combining
Empirical Mode Decomposition (EMD), Ensemble Empirical Mode Decomposition
(EEMD) and ARIMA for long-term streamflow forecasting, but they did not
examine the effectiveness of their models on datasets with extreme values.

Recently, deep learning models have emerged as the preferred approach for
forecasting rich time series data (31), outperforming classical statistical
approaches such as ARIMA or GARCH (32). Boris et al. (15) proposed NBeats,
which shows good performance on general time series prediction (33). DeepAR
(33) learns a conditional distribution over the future values and uses the
shared RNN to predict future values and their confidence. To tackle the long-
term forecasting challenge, some recent transformer-based methods like
Autoformer (34) and Reformer (22) have been proposed to empower the
transformer with more sophisticated dependency discovery and representation
ability. Informer (20) proposed a ProbSparse self-attention mechanism and a
generative style decoder which drastically improves the inference speed of
long-sequence predictions. FEDFormer (21) represents time series by randomly
selecting a constant number of Fourier components to maintain the global
property and statistics of time series as a whole.

On the other hand, many generic time series prediction models can perform
poorly on data with high skewness and kurtosis scores. Conventional methods
often falter when confronted with extreme events, which, although infrequent,
hold considerable real-world implications—such as in specific instances of
streamflow forecasting. Singh et al. (35) proved machine learning approaches
suffer from the problem of imbalanced data distribution and noted that
balancing the dataset is an imperative sub-task. An and Cho (36) proposed an
anomaly detection method using the reconstruction probability, which is a
probabilistic measure that takes into account the variability of the data
distribution. Ding et al. (2) explored the central theme of improving the
ability of deep learning time series models to capture extreme events. Zhang
et al. (5) proposed a framework to integrate machine learning models with
anomaly detection algorithms. In an earlier work, we proposed NEC+ (6), a
model specifically designed to provide good prediction performance on
hydrological time series with extreme events. Additionally, concurrently with
this work, we recently developed SEED (37), a Segment-Expandable Encoder-
Decoder model for univariate streamflow prediction.

Few of these prior works have concentrated on addressing both prolonged
sequences and extreme occurrences. To bridge this gap, the proposed method,
DAN, employs joint learning of two polar hidden spaces within a single series,
guided by an associated distance-regularized loss function. This combined
approach facilitates accurate end-to-end structure prediction. Remarkably,
DAN’s performance surpasses that of state-of-the-art methods across four real-
life streamflow datasets.

Figure 2: DAN’s end-to-end extendable framework consists of two stages, named
RepGen and RepMerg, respectively. RepGen contains three parallel encoder-
decoder blocks, resulting in polar representations of ordinary series inputs
and refined indicators. These elements are further merged in the RepMerg
stack.

## Preliminaries

###  0.1 Problem Statement

Suppose we have a collection of $m$ ($m>=1$) related univariate time series,
with each row corresponding to a different time series. We are going to
predict the next $h$ time steps for the first time series $x_{1}$, given
historical data from multiple length-$t$ observed series. The problem can be
described as,

$\begin{bmatrix}x_{\mathbf{1},1}&\cdots&x_{\mathbf{1},t}\\\
x_{2,1}&\cdots&x_{2,t}\\\ \vdots&\ddots&\vdots\\\
x_{m,1}&\cdots&x_{m,t}\end{bmatrix}\in\mathbb{R}^{m\times
t}\rightarrow[x_{\mathbf{1},t+1},...,x_{\mathbf{1},t+h}]\in\mathbb{R}^{h}$

where $x_{i,j}$ denotes the value of time series $i$ at time $j$. The matrix
on the left are the inputs, and $x_{1,t+1}$ to $x_{1,t+h}$ are the outputs in
our method.

We first define this task by modeling the objective time series $x_{1}$ as the
ordinary series and the group of related time series $x_{2}$ to $x_{m}$ as
extraordinary indicators. When an extraordinary indicator series is not
available, our proposed model can generate a Gaussian Mixture Model (GMM)
indicator based on $x_{1}$, which becomes $x_{2}$. In such cases, the problem
can be reduced to that of univariate time series forecasting.

###  0.2 GMM Indicator

In our work, when there is no extraordinary indicator series provided, we use
a Gaussian mixture model (GMM) (38) to learn a group of distributions from the
input univariate time series. Then, we compute an indicator feature for each
value in the time series as the weighted sum of all component probabilities,
based on the weights learned during GMM model fitting. Due to lack of space,
we detail this step in Section 4 of our technical appendix.

###  0.3 Kruskal-Wallis Test in Time Series

To balance the sparse distribution of extreme events, we employ the Kruskal-
Wallis test (39) as a non-parametric method to evaluate the normality of a
training sample and guide our oversampling policy. The Kruskal-Wallis test
examines two or more groups of time series based on their medians, in which
the data are first ranked, and the sum of ranks is calculated for each group.
The $H$ value is then calculated from these rank sums, and compared to a
critical value to determine if there are significant differences between the
groups. Because the Kruskal-Wallis test does not assume a particular
distribution, it is sometimes referred to as a distribution-free test (40,
41). The $H$ value is computed as

$H=\frac{12}{n(n+1)}\sum_{j=1}^{k}\frac{R_{j}^{2}}{n_{j}}-3(n+1).$ (1)

In our work, we separate the sampled training sequence into $k$ sub sequence
groups of equal length. $H$ is the Kruskal-Wallis test statistic, $n$ is the
total number of samples across all groups, $R_{j}$ is the sum of ranks for the
$j$th group, and $n_{j}$ is the number of samples in the $j$th group. The
oversampling policy will be described in Section 0.7.

## Methods

An overview of our architecture has been presented in Figure 2.

###  0.4 Polar Representation Learning

The key innovations in DAN are new mechanisms to generate and exchange
information among polar (far and near) representations and the indicators
(ind) for direct improvement of the predictions. Polar representation learning
in RepGen allows for the separate encoding of extreme points, which are then
preserved in RepMerg. This ensures that, during training, these
representations are not compromised by the predominance of normal values, as
they often adhere to different distributions. As described in Figure 3, RepGen
consists of three encoder-decoder sub-networks. Please note that the CONV-LSTM
layers are specifically designed for learning to predict far points. Because
the hidden state for extreme events may be updated multiple times in repeated
blocks, we use convolution operations to shorten the input sequence, which
helps alleviate any potential exploding or vanishing gradient problem. The far
and near point predictions are extended into the RepMerg stack and further
refined as $\hat{y}_{f}$ and $\hat{y}_{n}$, respectively, after being added to
the prior output of RepGen. To reflect the change of predicted values,
$\hat{y}_{i}$ is expected to converge to the first order of y.

###  0.5 Architecture Items

Figure 3:  In the RepGen stack, “f-ED” is responsible for the representation
learning of those points that are far away from the mean of the series
($\hat{y}_{f}$), including the extreme values sparsely distributed in the
predicted zone. “n-ED” mainly focuses on learning hidden features of near
points ($\hat{y}_{n}$), which include most of the normal values. “i-ED” is
designed to learn the indicator representation ($\hat{y}_{i}$).

##### Representation Merge

In the RepMerg stack, to finish representation merging, the middle LSTM-FC
block in the red circle takes the far point representation $h_{far}$ as the
initial hidden state and $\hat{y}_{i}$ as input, which is then combined with
$\hat{y}_{n}$ to form a two-dimensional vector as input for another FC layer
that outputs $\hat{y}$.

##### Indicator Refine Layer

It should be highlighted that precise $\hat{y}_{i}$ prediction is important
for performance improvement as it can help predict extreme events. Therefore,
an additional refine layer made of double CNN components is intended to assist
in producing the expected indicator by first refining the value within a local
horizon. The indicator input can then be cyclically updated from the output of
the refine layer of the “i-ED” as the RepGen is designed to be an expandable
stack that can be repeated multiple times.

##### Gate Control Vector

Given the significance of $\hat{y}_{i}$, we provide another way to hone the
predicted indicator. As shown in Figure 3, we designed a gate control vector
$m_{far}$, whose values reflect precisely in which places $\hat{y}$ is closer
to $\hat{y}_{f}$. Similarly, $m_{near}$ is the complement of $m_{far}$. We
then compute $\hat{y}_{w}$ as
$m_{far}\odot\hat{y}_{f}+m_{near}\odot\hat{y}_{n}$, where $\odot$ is the
component-wise multiplication. We therefore use $\hat{y}_{w}$ to increase the
discriminatory ability of $\hat{y}_{i}$ when it is expected to reflect better
extreme values without sacrificing the overall normal values.

Figure 4: Gate control vector. $m_{far}$ is equal to
$sigmoid(\alpha\hat{y}_{i})$, where $\alpha$ is an amplifier and equals $4$ in
our experiments; $m_{near}$ is computed as $1-m_{far}$.

.

###  0.6 Auto-Regularized Loss Function

Different from the conventional usage of regularization loss as a penalty term
for preventing overfitting of the model to the training data, our approach
employs multiple distance-weighted loss functions when training the DAN model,
with the objective of compelling the model to learn more informative
representations. Moreover, it should be noted that our method can also serve
as an effective regularizer for preventing overfitting of the model to the
base normal values in the long-term time series prediction task.

We define $w_{f}=(\tanh(y))^{2}$ to be a weight that emphasizes the accuracy
of points further away from the mean value of 0 (since the series were
standardized to have 0 mean). In contrast, $w_{n}=(1-|\tanh(y)|)^{2}$ focuses
more on the points closer to zero. The square root of $w_{f}$, denoted by
$w_{p}$, is a more moderate way to maintain discriminatory output for
indication-related tasks. These weight definitions contribute to the accuracy
of the model in predicting extreme events in a long-term time series. Based on
these weights, we build our multi regulation loss as follows,

$\displaystyle\mathcal{L}_{1}$ $\displaystyle=RMSE((\hat{y}_{f}\odot
w_{f}),(y\odot w_{f})),$ $\displaystyle\mathcal{L}_{2}$
$\displaystyle=RMSE((\hat{y}_{n}\odot w_{n}),(y\odot w_{n})),$
$\displaystyle\mathcal{L}_{3}$ $\displaystyle=RMSE(\hat{y}_{w}\odot
w_{p},y\odot w_{p}),$ $\displaystyle\mathcal{L}_{4}$
$\displaystyle=RMSE(\hat{y}_{i}\odot w_{p},y_{i}\odot w_{p}),$

where $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ are used to regulate the bipolar
representation learning and $\mathcal{L}_{3}$ and $\mathcal{L}_{4}$ force the
predicted indicator to reflect the change of predicted values by setting
$y_{i}$ equal to the first order of $y$. Then, the overall loss is composed
as,

$\mathcal{L}=RMSE(\hat{y},y)+\lambda\times(\mathcal{L}_{1}+\mathcal{L}_{2}+\mathcal{L}_{3}+\mathcal{L}_{4}),$

where $\lambda$ is a multiplier ($\lambda=max(-1\cdot
e^{\frac{epoch}{45}}+2,0.2)$ in our experiments) applied on those regulation
items, decreasing with each epoch.

###  0.7 Kruskal-Wallis Sampling

Given that extreme events are rare within our data compared with normal ones,
we utilize Kruskal-Wallis sampling to over-sample regions with extreme events
in our training set that our model can learn appropriate patterns from.
Namely, for each random sample $x$ of size $t+h$ we draw from the input
sequence, we first split the sequence into $k$ consecutive sub-sequences of
equal size and compute the Kruskal-Wallis test statistic $H$ between the $k$
sub-sequences, using Equation 1. To avoid the $H$ statistic being affected by
minor differences in the sub-sequences, we round values in $x$ to the nearest
integer before computing $H$. We then include the sample in the training set
if $H>\epsilon$, where $\epsilon$ is a threshold, or otherwise include the
sample with probability $p<1$. The threshold $\epsilon$ allows us to set the
relative change in the sample that makes it more likely to contain an extreme
event, while the probability $p$ allows us to choose how many normal samples
should be included in the training set.

## Evaluation

Code and data for DAN can be found at https://github.com/davidanastasiu/dan.
In this section, we present empirical results from evaluating our proposed
framework. We are interested in answering the following research questions
with regards to prediction effectiveness: (1) How does DAN compare against
state-of-the-art baselines? (2) What is the effect of DAN’s extensible
framework? (3) What is the effect of the Kruskal-Wallis oversampling policy?
(4) How do the critical design elements of the framework affect performance?

###  0.8 Experimental Settings

##### Dataset

We used four groups of hydrologic datasets from Santa Clara County, CA, namely
Ross, Saratoga, UpperPen, and SFC, named after their respective locations.
Each group included a streamflow dataset and an associated rainfall dataset.
Statistics of our primary series are shown in Table 1. Our task was to
forecast the streamflow for wet seasons in a hydrologic year, excluding the
summer months, namely September 2021 to May 2022, in a rolling manner. The
training and validation sets were randomly sampled from series spanning from
January 1988 to August 2021. Inference involved predicting the streamflow
every 4 hours for the next 3 days. Since the sensors measure the streamflow
and precipitation every 15 minutes, we are attempting a lengthy forecasting
horizon ($h=288$), which is unquestionably an LSTF task based on the most
recent research (20, 21, 23, 24). Before training a model, all time series
were pre-processed by log transform ($x_{i}=\log(1+x_{i})\forall i$) and
standardization (subtract mean and divide by standard deviation). Inference
predictions were post-processed by inverting the standardization and log
transform operations.

Table 1: Streamflow Datasets Statistics

| Ross | Saratoga | UpperPen | SFC
---|---|---|---|---
min | 0.00 | 0.00 | 0.00 | 0.00
max | 1440.00 | 2210.00 | 830.00 | 7200.00
mean | 2.91 | 5.77 | 6.66 | 20.25
skewness | 19.84 | 19.50 | 13.42 | 18.05
kurtosis | 523.16 | 697.78 | 262.18 | 555.18

High skewness and kurtosis scores indicate that there is a significant
deviation from a normal distribution in our data.

##### Model Parameters

For all models, after testing different LSTM layer widths, we found that 512
node layers for ROSS and 384 layers for the other three sensors were the most
effective. We set $h=288$ (3 days) based on the problem definition and tested
different values of $t\in\\{288,672,1440\\}$ (3, 7, 15 days), with $t=1440$
producing the best results for all data streams. In the RepGen stage, the
three CNN layers produce 256 channels each. The kernel sizes used in these
layers are 11, 7, and 3, respectively. The stride, padding, and activation
function remain the same across all three layers, with a stride equal to the
kernel size, no padding, and a subsequent $tanh$ activation function. We use
two stacked CNN1d layers for indicator refinement, with the kernel size and
padding set to 7 for the first layer and 3 for the second.

###  0.9 Experimental Results

##### Baselines

We include seven state-of-the-art models for comparison, of which FEDFormer
(21), Informer (20), NLinear (23), and DLinear (23) focus on long-term time
series forecasting, while NEC+ (6) holds the best performance for hydrologic
time series prediction in the presence of extreme events. These five models
were used as baselines for both multivariate and univariate prediction. In
addition, Attention-LSTM (42) was used as a state-of-the-art hydrologic
multivariate model used to predict stream-flow using rainfall data. Finally,
N-BEATS (15), a state-of-the-art univariate baseline method that outperformed
all competitors on the standard M3 (43), M4 (44), and TOURISM (45) time series
datasets, was also used in the comparison.

Table 2: Multivariate/Univariate Long-Term ($h=288$) Series Forecasting
Results on Four Datasets

Methods | Metric | Ross | Saratoga | UpperPen | SFC
---|---|---|---|---|---
|  | Multi | Single | Multi | Single | Multi | Single | Multi | Single
FEDformer | RMSE | 6.01 | 6.49 | 6.01 | 6.85 | 3.05 | 2.38 | 23.54 | 24.10
| MAPE | 2.10 | 2.49 | 1.55 | 2.26 | 1.87 | 1.02 | 2.35 | 2.817
Informer | RMSE | 7.84 | 9.14 | 5.04 | 4.89 | 5.88 | 5.33 | 39.89 | 19.00
| MAPE | 4.05 | 5.45 | 1.43 | 0.73 | 4.10 | 4.21 | 8.64 | 0.54
Nlinear | RMSE | 6.10 | 5.84 | 5.23 | 4.98 | 1.57 | 1.74 | 18.47 | 18.43
| MAPE | 1.99 | 1.62 | 0.83 | 0.75 | 0.45 | 0.57 | 0.92 | 0.87
Dlinear | RMSE | 7.16 | 6.90 | 4.33 | 4.06 | 3.53 | 3.25 | 21.62 | 23.64
| MAPE | 3.10 | 2.79 | 1.40 | 1.31 | 2.35 | 2.04 | 2.74 | 4.02
NEC+ | RMSE | 9.44 | 9.33 | 1.88 | 1.95 | 2.22 | 1.94 | 17.00 | 16.39
| MAPE | 4.80 | 4.53 | 0.17 | 0.21 | 0.95 | 0.80 | 1.07 | 0.55
LSTM-Atten / | RMSE | 7.35 | 5.16 | 6.49 | 3.60 | 6.35 | 1.23 | 34.17 | 31.47
NBeats | MAPE | 3.74 | 1.25 | 1.80 | 0.70 | 4.76 | 0.25 | 9.90 | 3.24
DAN | RMSE | 4.25 | 4.24 | 1.80 | 1.84 | 1.10 | 1.31 | 15.23 | 15.20
| MAPE | 0.07 | 0.09 | 0.14 | 0.16 | 0.15 | 0.32 | 0.26 | 0.21

Over 1600 test points in the test set were inferenced on all datasets. The
best results are in bold and the second best results are underlined.

##### Multivariate and Univariate Results

Figure 5: Sampled multivariate inference for the Ross sensor. The right figure
emphasizes extreme events occurring during the example time period in the red
box of the left figure.

Table 2 shows the test root mean squared error (RMSE) and mean absolute
percentage error (MAPE) performance for the models that achieved the best
performance on our validation dataset. For these experiments, all DAN results
were achieved using the same random seed. Due to lack of space, we include the
definition of our performance metrics and multiple seed run performance
statistics in Section 5 of our technical appendix. In the multivariate
forecasting task, our proposed model DAN outperformed all baselines on all
four benchmark datasets. Compared to the second-best results, DAN achieved an
overall 19% relative RMSE reduction. Notably, the improvement was most
significant for the UpperPen dataset, where DAN achieved a 30% improvement.
For the univariate forecasting task, DAN outperformed other methods on three
out of four benchmark datasets. Although NBeats achieved a 6% relative RMSE
reduction compared to DAN for the UpperPen sensor, DAN surpassed NBeats with a
relative RMSE reduction of 18%, 49%, and 52% on the Ross, Saratoga, and SFC
datasets, respectively.

##### Inference Overall Analysis

Figure 5, in which we present rolling prediction results for the whole test
set (1600 time points) for the Ross sensor, helps explain DAN’s good
performance. To make it easier to visualize, the test set was sampled every 5
points (320 time points sampled) and a specific period including extreme
events is denoted by the red box and focused on in the right figure. We
observed that DAN performed better than other models on areas with extreme
events. While Informer and NLinear could follow the extreme events better than
other baselines, they predicted values much higher than the actual peaks. On
the other hand, DLinear and Attention-LSTM performed better than NEC+ and
FEDFormer, but they predicted values much lower than the ground truth.

###  0.10 Ablation Studies

The superior performance of our method comes from its extendable framework,
creative Kruskal-Wallis sampling policy, and three key architecture designs:
1) the use of a refine layer in RepGen that refines the indicator information
using double CNNs with moving kernel convolution operations, 2) the RepMerg
layer which combines the rich polar representations with the indicator learned
from RepGen, and 3) regularization loss components. We will examine their
effects respectively in this section.

Figure 6: Oversampling policy. Initially, a majority of the samples have $H$
values below 10. However, upon setting $p=50\%$ and $\epsilon=10$, the number
of samples with $H$ values below 10 decreases, while the count for higher $H$
values increases. The second row showcases the scenario for $p=25\%$.

####  0.10.1 Effects of DAN’s Extensible Framework

Our network comprises two main stages. The repeatable stacks are named “E”,
“D” and “R” in Figure 2, endowing DAN with the capability of updating the
indicator in an evolutional way by repeating “E”+“D” and including multiple
refinements by repeating “R”. The best stack configuration may vary for
different datasets, depending on the intrinsic relationships in the
multivariate series. We experimented with various combinations and identified
the best results as “EDEDRR”, “EDR”, “EDEDRR”, and “EDEDR” for Ross, Saratoga,
UpperPen, and SFC, respectively.

####  0.10.2 Effects of the Oversampling Policy

Figure 6 displays the distribution of the $H$ values before and after applying
the Kruskal-Wallis sampling algorithm. These observations highlight the impact
of adjusting $p$ and $\epsilon$ on the distribution of $H$ values in the
training set. If we maintain the $p$ value and increase the $\epsilon$ value,
the training set will contain more samples with $H$ values exceeding
$\epsilon$, as illustrated in the rightmost three figures in the first row of
Figure 6.

Table 3: $RMSE_{far}$ when Oversampling

Policy($\epsilon=10$) | Ross | SFC
---|---|---
$20\%<=p<=40\%$ | 23.5 | 118.1
$60\%<=p<=80\%$ | 25.0 | 122.0
$p=100\%$ | 28.7 | 129.0

We used $\epsilon$ = 10 and the $p$ values used were grouped into 3 cases.
$p=100\%$ means no Kruskal-Wallis oversampling was applied.

Our oversampling policy can help shift the focus of our model towards
improving “far” point prediction performance. To test this, we conducted
multiple runs of each model, averaging the RMSE values for points with values
greater than 1.5 standard deviations above the mean of the series, which we
call $RMSE_{far}$ in Table 3. The results show that the $RMSE_{far}$ can be
steadily decreased by decreasing $p$.

####  0.10.3 Effects of DAN Architecture Elements

Figure 7: Inference examples to show the effects of different architecture
elements on the Ross dataset.

In this section, we investigate the impact of different design elements on the
performance of our network. To isolate the effects of these design elements
and obtain a more comprehensive understanding of each one, we replace the
input indicator in our network with ground truth rain data and only use a
simple EDR architecture stack, which will remove the effect of varying
indicator performance from the comparison. By doing so, we can observe how the
network performs when certain layers are removed or added back.

To create a baseline structure, we first removed the refine and the RepMerge
layers, using only the encoder-decoder block to generate the indicator, and
set $\hat{y}$ equals to $\hat{y}_{w}$ to bypass the RepMerge structure. This
produced the first predicted sequence (named Vanilla_Ind in Figure 7), which
shows an example of inference on Ross sensor data. We then added the refine
layer back, which improved the results as shown in the second sub-figure
(named LCNN_Ind). Next, we added back the RepMerg layer. Concurrently, we also
removed the regularization loss items except $\mathcal{L}_{3}$, obtaining the
third figure (named Un-regulated). Finally, we added all regularization loss
items back, which gave the best result, as shown in the fourth figure.

These experiments demonstrate that the double CNNs with moving kernel
convolutional operations refines the indicator information, and RepMerg
produces better results as the Gate control vector mechanism increases the
discrimination of predicted values. Adding polar representation of the basic
series assists in identifying data distributions beyond the indicator
information and enhances the accuracy of data at corners of each fluctuation,
as denoted in the blue circles in Figure 7. Therefore, including the refine
layer, representation learning, and gate control vector resulted in the best
performance.

## Conclusion

In this work, we presented a novel end-to-end framework, DAN, designed to
better account for rare yet important extreme events in long single- and
multi-variate streamflow time series. Our framework learns polar
representations for predicting extreme and normal values, along with a
representation merging model that makes prediction in an expandable way. In
addition, to improve training performance, our framework uses Kruskal-Wallis
sampling policies to accommodate imbalanced extreme data and a distance-
weighted multi-loss regularization penalty. Extensive experiments using more
than 33 years of streamflow data from Santa Clara County, CA, showed that DAN
provides significantly better predictions than state-of-the-art baselines.

## Appendix

##  1 Network description

For all models, we used 10 bidirectional LSTMs with 512 hidden features for
ROSS and 384 for the others. Figure 8 shows a detailed representation of our
framework.

Figure 8: Data flow and architecture layers.

In the RepGen stage, the three CNN layers are labelled as C0, C1 and C2. The
stride, padding, and activation functions are listed in Table 4. Two stacked
CNN1d layers are used to finish the indicator refinement, which we noted as
C3. Each layer has 256 channels with a stride of 1. The kernel size and
padding for the first two layers are set to 7 and 3, respectively. All the
fully connected layers have $384$ (512 for ROSS)$\times 2$ input features and
1 output feature.

Table 4: Convolution Layers in Our Framework layer | in_dim | out_dim | kernel size | stride | padding | acti func
---|---|---|---|---|---|---
C0 | 1 | 256 | 11 | 11 | 0 | tanh
C1 | 1 | 256 | 7 | 7 | 0 | tanh
C2 | 1 | 256 | 3 | 3 | 0 | tanh
C3 | $512/384\times 2$ | 512/384 | 7 | 1 | 3 | relu
512/384 | 1 | 3 | 1 | 1
Table 5: LSTM Layers in Our Framework layer | in_dim | out_dim | layer | bidirectional
---|---|---|---|---
L1 | 256 | 512/384 | 1 | yes
L2 | 256 | 512/384 | 1 | yes
L3 | 1 | 512/384 | 1 | yes
L4 | 1 | 512/384 | 1 | yes
L5 | 2 | 512/384 | 1 | yes
L6 | 2 | 512/384 | 1 | yes
L7 | 2 | 512/384 | 1 | yes
L8 | 2 | 512/384 | 1 | yes
L9 | 1 | 512/384 | 1 | yes
L10 | 2 | 512/384 | 1 | yes

##  2 Dataset description

To analyze our data set, we computed the related statistical values of our
input time series. We use the measures of skewness and kurtosis to provide
valuable insights into the shape and distribution of the data. The high
skewness and kurtosis shows that the datasets deviate significantly from a
normal distribution. In our datasets, the positive skewness reflects the
asymmetry of the distribution. To be specific, it indicates a longer right
tail. On the other hand, the kurtosis values in the hundreds indicate heavier
tails and a higher peak.

Table 6: Statistics of our Stream Data Statistic / Stream | Ross | Saratoga | UpperPen | SFC
---|---|---|---|---
mean | 2.91 | 5.77 | 6.66 | 20.25
max | 1440.00 | 2210.00 | 830.00 | 7200.00
min | 0.00 | 0.00 | 0.00 | 0.00
median | 0.17 | 1.00 | 3.20 | 1.20
variance | 597.22 | 711.09 | 452.90 | 12108.14
std_deviation | 24.43 | 26.66 | 21.28 | 110.03
skewness | 19.84 | 19.50 | 13.42 | 18.05
kurtosis | 523.16 | 697.78 | 262.18 | 555.18

Table 6 shows the computed statistics of our input time series. From these
data, we can find these data sets have a pronounced departure from normality,
with a significant skew and heavy tails. This implies that the distribution of
the data is skewed towards one side and contains extreme values that deviate
from the average.

##  3 Related research

In this section, we discuss several studies that have employed traditional
machine learning methods for stream flow prediction. Nayak et al. (46)
utilized neural network models (NN) and adaptive neuro-fuzzy inference systems
(ANFIS) to forecast stream flow. Tabbussum et al (47) proposed a fuzzy
inference system termed as flood model using two fuzzy flood models. Gaussian
process regression (GPR) and quantile regression were used in some studies to
not only predict but also quantify forecast uncertainty. Tree-based models,
such as classification and regression trees (CARTs) and random forest (RF),
have been employed due to their computational efficiency and ability to handle
predictors without assuming any specific distribution. Nhu et al. (48)
demonstrated the use of CART and RF in solving hydrologic prediction problems.
While neural networks have been utilized in early hydrologic flow prediction
studies, they were often shallow networks unable to capture complex patterns
in the data, requiring extensive feature engineering and manual tuning based
on domain expertise to improve performance.

##  4 GMM Indicator Generation

A Gaussian mixture model (GMM) (49) can be described by the equation,

$p(x|\lambda)=\sum_{i=1}^{M}w_{i}\leavevmode\nobreak\
g(x|\mu_{i},\mathbf{\Sigma_{i}}),$

where $x$ is a $D$-dimensional continuous-valued vector, $w_{i}\ \forall
i=1,\ldots,M$ are the mixture weights, and $g(x|\mu_{i},\Sigma_{i})$, are the
component Gaussian densities. Each component density is a $D$-variate Gaussian
function, and the overall GMM model is a weighted sum of M learnable Gaussian
densities, i.e.,

$g(x|\mathbf{\mu_{i}},\mathbf{\Sigma_{i}})=\frac{1}{2\pi^{\frac{D}{2}}|\Sigma_{i}|^{\frac{1}{2}}}\exp{\left\\{-\frac{1}{2}{(x-\mu_{i})}^{T}\leavevmode\nobreak\
\Sigma_{i}^{-1}(x-\mu_{i})\right\\}},$

where $\mu_{i}$ is the mean vector and $\Sigma_{i}$ is the covariance matrix
of the $i$th component. The mixture weights are constrained such that
$\sum_{i=1}^{M}w_{i}=1$. The GMM’s capacity to produce smooth approximations
to arbitrarily shaped densities is one of its most impressive features.

In our work, we utilize the Expectation-Maximization (EM) algorithm to fit a
Gaussian Mixture Model (GMM) to the ordinary series training data. The GMM
consists of multiple components, and each component represents a Gaussian
distribution. The EM algorithm iteratively estimates the parameters of the GMM
by maximizing the likelihood of the observed data.

Once we have fitted the GMM, each model component can generate a probability
for each point in the time series. These probabilities represent the
likelihood that a particular component generated the observed data point. To
further analyze the time series, we compute an "indicator" feature for each
value in the time series. The indicator feature is computed as the weighted
sum of all component probabilities, given the weights learned during the
fitting process of the GMM model. This feature provides a measure of the
contribution of each component to the overall time series at a specific point.

In our framework, the number of components in the GMM, denoted as $M$, is a
hyper-parameter that can be tuned for each time series. We set it as 3 in all
our experiments.

##  5 Experiments

###  5.1 Evaluation metrics

We use the popular Root Mean Square Error (RMSE) and Mean Absolute Percentage
Error (MAPE) metrics to measure the effectiveness of our prediction methods.
RMSE is defined as the square root of the mean squared error between the
predicted and actual values, while MAPE is the mean of all absolute percentage
errors between the predicted and actual values.

###  5.2 Multiple seed performance statistics

Table 7: Multiple seed performance statistics RMSE | Ross | Saratoga | UpperPen | SFC
---|---|---|---|---
multi | single | multi | single | multi | single | multi | single
min | 4.24 | 4.23 | 1.71 | 1.82 | 1.03 | 1.12 | 14.99 | 15.20
max | 4.36 | 4.35 | 1.96 | 2.06 | 1.35 | 1.31 | 15.86 | 17.56
mean | 4.29 | 4.31 | 1.84 | 1.91 | 1.18 | 1.19 | 15.55 | 16.34
std_deviation | 0.04 | 0.02 | 0.08 | 0.09 | 0.12 | 0.06 | 0.26 | 0.55

In this part of the study, we aim to mitigate the impact of randomness and
obtain a more comprehensive evaluation by training our models multiple times
with different random seeds. Specifically, we performed 6 to 8 runs for the
multivariate models and the univariate models. The resulting minimum, maximum,
mean, and standard deviation values are summarized in the Table 7. The
oversampling parameters used in this experiment are $p=25\%,60\%,18\%,18\%$
and $\epsilon=30,150,80,80$ for the Ross, Saratoga, UpperPen, and SFC streams,
respectively.

Despite the average RMSE being slightly higher than the minimum values, it is
important to note that DAN consistently outperformed all baselines across all
four benchmark datasets.

###### Acknowledgements.

Funding for this project was made possible by the Santa Clara Valley Water
District.

## Bibliography

## References

  * Nguyen and Chan (2004) Hanh H Nguyen and Christine W Chan.  Multiple neural networks for a long term time series forecast.  _Neural Computing & Applications_, 13:90–98, 2004.
  * Ding et al. (2019) Daizong Ding, Mi Zhang, Xudong Pan, Min Yang, and Xiangnan He.  Modeling extreme events in time series prediction.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 1114–1122, 2019.
  * Qi and Majda (2020) Di Qi and Andrew J Majda.  Using machine learning to predict extreme events in complex systems.  _Proceedings of the National Academy of Sciences_ , 117(1):52–59, 2020.
  * Chen et al. (2020) Zechuan Chen, Haomin Yu, Yangli-ao Geng, Qingyong Li, and Yingjun Zhang.  Evanet: An extreme value attention network for long-term air quality prediction.  In _2020 IEEE International Conference on Big Data (Big Data)_ , pages 4545–4552. IEEE, 2020.
  * Zhang et al. (2021) Yifan Zhang, Jiahao Li, Ablan Carlo, Alex K Manda, Scott Hamshaw, Sergiu M. Dascalu, Frederick C. Harris, and Rui Wu.  Data regression framework for time series data with extreme events.  In _2021 IEEE International Conference on Big Data (Big Data)_ , pages 5327–5336, 2021.  10.1109/BigData52589.2021.9671387.
  * Li et al. (2023a) Yanhong Li, Jack Xu, and David C. Anastasiu.  An extreme-adaptive time series prediction model based on probability-enhanced lstm neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 2023a.
  * Box and Pierce (1970) George EP Box and David A Pierce.  Distribution of residual autocorrelations in autoregressive-integrated moving average time series models.  _Journal of the American statistical Association_ , 65(332):1509–1526, 1970.
  * Nielsen (2019) Aileen Nielsen.  _Practical time series analysis: Prediction with statistics and machine learning_.  O’Reilly Media, 2019.
  * Box and Jenkins (1976) George.E.P. Box and Gwilym M. Jenkins.  _Time Series Analysis: Forecasting and Control_.  Holden-Day, , 1976.
  * Shortridge et al. (2016) Julie E Shortridge, Seth D Guikema, and Benjamin F Zaitchik.  Machine learning methods for empirical streamflow simulation: a comparison of model accuracy, interpretability, and uncertainty in seasonal watersheds.  _Hydrology and Earth System Sciences_ , 20(7):2611–2628, 2016.
  * Papacharalampous and Tyralis (2022) Georgia Papacharalampous and Hristos Tyralis.  A review of machine learning concepts and methods for addressing challenges in probabilistic hydrological post-processing and forecasting.  _Frontiers in Water_ , 4:961954, 2022.
  * Cheng et al. (2020) M Cheng, F Fang, T Kinouchi, IM Navon, and CC Pain.  Long lead-time daily and monthly streamflow forecasting using machine learning methods.  _Journal of Hydrology_ , 590:125376, 2020.
  * Lai et al. (2018) Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu.  Modeling long-and short-term temporal patterns with deep neural networks.  In _The 41st international ACM SIGIR conference on research & development in information retrieval_, pages 95–104, 2018.
  * Siami-Namini et al. (2018) Sima Siami-Namini, Neda Tavakoli, and Akbar Siami Namin.  A comparison of arima and lstm in forecasting time series.  In _2018 17th IEEE international conference on machine learning and applications (ICMLA)_ , pages 1394–1401. IEEE, 2018.
  * Oreshkin et al. (2019) Boris N Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-beats: Neural basis expansion analysis for interpretable time series forecasting.  _arXiv preprint arXiv:1905.10437_ , 2019.
  * Wu et al. (2020) Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang.  Connecting the dots: Multivariate time series forecasting with graph neural networks.  In _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_, pages 753–763, 2020.
  * Cao et al. (2020) Defu Cao, Yujing Wang, Juanyong Duan, Ce Zhang, Xia Zhu, Congrui Huang, Yunhai Tong, Bixiong Xu, Jing Bai, Jie Tong, et al.  Spectral temporal graph neural network for multivariate time-series forecasting.  _Advances in neural information processing systems_ , 33:17766–17778, 2020.
  * Li et al. (2019) Shiyang Li, Xiaoyong Jin, Yao Xuan, Xiyou Zhou, Wenhu Chen, Yu-Xiang Wang, and Xifeng Yan.  Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting.  _Advances in neural information processing systems_ , 32, 2019.
  * Qin et al. (2017) Yao Qin, Dongjin Song, Haifeng Chen, Wei Cheng, Guofei Jiang, and Garrison Cottrell.  A dual-stage attention-based recurrent neural network for time series prediction.  _arXiv preprint arXiv:1704.02971_ , 2017.
  * Zhou et al. (2021) Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang.  Informer: Beyond efficient transformer for long sequence time-series forecasting.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 35, pages 11106–11115, 2021.
  * Zhou et al. (2022) Tian Zhou, Ziqing Ma, Qingsong Wen, Xue Wang, Liang Sun, and Rong Jin.  Fedformer: Frequency enhanced decomposed transformer for long-term series forecasting.  In _International Conference on Machine Learning_ , pages 27268–27286. PMLR, 2022.
  * Kitaev et al. (2020) Nikita Kitaev, Łukasz Kaiser, and Anselm Levskaya.  Reformer: The efficient transformer.  _arXiv preprint arXiv:2001.04451_ , 2020.
  * Zeng et al. (2022) Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu.  Are transformers effective for time series forecasting?  _arXiv preprint arXiv:2205.13504_ , 2022.
  * Das et al. (2023) Abhimanyu Das, Weihao Kong, Andrew Leach, Rajat Sen, and Rose Yu.  Long-term forecasting with tide: Time-series dense encoder.  _arXiv preprint arXiv:2304.08424_ , 2023.
  * Fortuin et al. (2018) Vincent Fortuin, Matthias Hüser, Francesco Locatello, Heiko Strathmann, and Gunnar Rätsch.  Som-vae: Interpretable discrete representation learning on time series.  _arXiv preprint arXiv:1806.02199_ , 2018.
  * Lei et al. (2017) Qi Lei, Jinfeng Yi, Roman Vaculin, Lingfei Wu, and Inderjit S Dhillon.  Similarity preserving representation learning for time series clustering.  _arXiv preprint arXiv:1702.03584_ , 2017.
  * Tonekaboni et al. (2022) Sana Tonekaboni, Chun-Liang Li, Sercan O Arik, Anna Goldenberg, and Tomas Pfister.  Decoupling local and global representations of time series.  In _International Conference on Artificial Intelligence and Statistics_ , pages 8700–8714. PMLR, 2022.
  * Ma et al. (2021) Haojie Ma, Zhijie Zhang, Wenzhong Li, and Sanglu Lu.  Unsupervised human activity representation learning with multi-task deep clustering.  _Proceedings of the ACM on Interactive, Mobile, Wearable and Ubiquitous Technologies_ , 5(1):1–25, 2021.
  * Ma et al. (2018) Jiaqi Ma, Zhe Zhao, Xinyang Yi, Jilin Chen, Lichan Hong, and Ed H Chi.  Modeling task relationships in multi-task learning with multi-gate mixture-of-experts.  In _Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining_, pages 1930–1939, 2018.
  * Wang et al. (2018) Zhi-Yu Wang, Jun Qiu, and Fang-Fang Li.  Hybrid models combining emd/eemd and arima for long-term streamflow forecasting.  _Water_ , 10(7), 2018.  ISSN 2073-4441.  10.3390/w10070853.
  * Sen et al. (2019) Rajat Sen, Hsiang-Fu Yu, and Inderjit S Dhillon.  Think globally, act locally: A deep neural network approach to high-dimensional time series forecasting.  _Advances in neural information processing systems_ , 32, 2019.
  * Box et al. (2015) George EP Box, Gwilym M Jenkins, Gregory C Reinsel, and Greta M Ljung.  _Time series analysis: forecasting and control_.  John Wiley & Sons, 2015.
  * Salinas et al. (2020) David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski.  Deepar: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ , 36(3):1181–1191, 2020.
  * Wu et al. (2021) Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long.  Autoformer: Decomposition transformers with auto-correlation for long-term series forecasting.  _Advances in Neural Information Processing Systems_ , 34:22419–22430, 2021.
  * Singh et al. (2022) Amit Singh, Ranjeet Kumar Ranjan, and Abhishek Tiwari.  Credit card fraud detection under extreme imbalanced data: a comparative study of data-level algorithms.  _Journal of Experimental & Theoretical Artificial Intelligence_, 34(4):571–598, 2022.
  * An and Cho (2015) Jinwon An and Sungzoon Cho.  Variational autoencoder based anomaly detection using reconstruction probability.  _Special lecture on IE_ , 2(1):1–18, 2015.
  * Li et al. (2023b) Yanhong Li, Jack Xu, and David C. Anastasiu.  Seed: An effective model for highly-skewed streamflow time series data forecasting.  In _2023 IEEE International Conference on Big Data (Big Data)_ , IEEE BigData 2023, Los Alamitos, CA, USA, Dec 2023b. IEEE Computer Society.
  * Day (1969a) N. E. Day.  Estimating the components of a mixture of normal distributions.  _Biometrika_ , 56(3):463–474, 12 1969a.  ISSN 0006-3444.  10.1093/biomet/56.3.463.
  * McKight and Najab (2010) Patrick E McKight and Julius Najab.  Kruskal-wallis test.  _The corsini encyclopedia of psychology_ , pages 1–1, 2010.
  * Ostertagova et al. (2014) Eva Ostertagova, Oskar Ostertag, and Jozef Kováč.  Methodology and application of the kruskal-wallis test.  In _Applied mechanics and materials_ , volume 611, pages 115–120. Trans Tech Publ, 2014.
  * Breslow (1970) Norman Breslow.  A generalized kruskal-wallis test for comparing k samples subject to unequal patterns of censorship.  _Biometrika_ , 57(3):579–594, 1970.
  * Le et al. (2021) Yan Le, Changwei Chen, Ting Hang, and Youchuan Hu.  A stream prediction model based on attention-lstm.  _Earth Science Informatics_ , 14:1–11, 06 2021.  10.1007/s12145-021-00571-z.
  * Makridakis and Hibon (2000) Spyros Makridakis and Michèle Hibon.  The m3-competition: results, conclusions and implications.  _International Journal of Forecasting_ , 16(4):451–476, 2000.  ISSN 0169-2070.  https://doi.org/10.1016/S0169-2070(00)00057-1.  The M3- Competition.
  * Makridakis et al. (2018) Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos.  The m4 competition: Results, findings, conclusion and way forward.  _International Journal of Forecasting_ , 34(4):802–808, 2018.
  * Athanasopoulos et al. (2011) George Athanasopoulos, Rob J Hyndman, Haiyan Song, and Doris C Wu.  The tourism forecasting competition.  _International Journal of Forecasting_ , 27(3):822–844, 2011.
  * Nayak et al. (2004) P.C Nayak, K.P Sudheer, D.M Rangan, and K.S Ramasastri.  A neuro-fuzzy computing technique for modeling hydrological time series.  _Journal of Hydrology_ , 291(1):52–66, 2004.  ISSN 0022-1694.  https://doi.org/10.1016/j.jhydrol.2003.12.010.
  * Tabbussum and Dar (2021) Ruhhee Tabbussum and Abdul Qayoom Dar.  Comparison of fuzzy inference algorithms for stream flow prediction.  _Neural Computing and Applications_ , 33:1643–1653, 2021\.
  * Nhu et al. (2020) Viet-Ha Nhu, Himan Shahabi, Ebrahim Nohani, Ataollah Shirzadi, Nadhir Al-Ansari, Sepideh Bahrami, Shaghayegh Miraki, Marten Geertsema, and Hoang Nguyen.  Daily water level prediction of zrebar lake (iran): A comparison between m5p, random forest, random tree and reduced error pruning trees algorithms.  _ISPRS International Journal of Geo-Information_ , 9(8), 2020.  ISSN 2220-9964.  10.3390/ijgi9080479.
  * Day (1969b) N. E. Day.  Estimating the components of a mixture of normal distributions.  _Biometrika_ , 56(3):463–474, 12 1969b.  ISSN 0006-3444.  10.1093/biomet/56.3.463.