# Auto-labelling of Bug Report using Natural Language Processing

Avinash Patil  Juniper Networks Inc.
Sunnyvale, USA
<EMAIL_ADDRESS>Aryan Jadon  San José State University
San Jose, USA
<EMAIL_ADDRESS>
###### Abstract

The exercise of detecting similar bug reports in bug tracking systems is known
as duplicate bug report detection. Having prior knowledge of a bug report’s
existence reduces efforts put into debugging problems and identifying the root
cause. Rule and Query-based solutions recommend a long list of potential
similar bug reports with no clear ranking. In addition, triage engineers are
less motivated to spend time going through an extensive list. Consequently,
this deters the use of duplicate bug report retrieval solutions. In this
paper, we have proposed a solution using a combination of NLP techniques. Our
approach considers unstructured and structured attributes of a bug report like
summary, description and severity, impacted products, platforms, categories,
etc. It uses a custom data transformer, a deep neural network, and a non-
generalizing machine learning method to retrieve existing identical bug
reports. We have performed numerous experiments with significant data sources
containing thousands of bug reports and showcased that the proposed solution
achieves a high retrieval accuracy of 70% for recall@5.

###### Index Terms:

Defect Reports, Bug Reports, Duplicate Detection, Similarity Search,
Information Retrieval, Machine Learning, Natural Language Processing, Text
Analysis.

##  I Introduction

A bug report, also known as a problem report, usually contains details that
help an engineer identify and reproduce issues in software or hardware. The
information can be structured or unstructured and the summary or description
of the issue constitute further detailed unstructured information. However,
impacted products, platforms, categories, etc., can be structured information.
Bug reports are managed using Bug Tracking Systems. These systems encounter
large volumes of duplicate reports constantly. A thorough root cause analysis
is needed before the problem report can be determined as a redundant report of
an already known issue. The larger the time gap between two reports, the
harder it is to identify the later filed report as duplicate. Typically, the
person triaging a new Bug Report determines if the similar issue at hand has
been reported in the past or not. To search for an existing bug report and
associated root cause analysis search query can be used to retrieve past
similar incidents. The query mechanism requires knowing the exact pattern,
string, or substring to look for. The difficulty with the traditional querying
approach is that every human describes the similar problem in different
manner. The summary and description of a problem vary when reported by
different users. This also applies when the same user is reporting the same
issue reencountered. The language and keywords used by the same person would
change over time. If the database query matches something, it will return all
the bug reports containing those patterns, strings, or substrings [3]. This is
not much help to the user as there is no way to determine how similar these
existing reports are to the new one unless done manually by going through each
result returned by the query. Moreover, identifying a duplicate bug report
before it enters the system results in but not limited to; preventing
repetitive root cause analysis, avoiding duplicate debug effort, enabling
auto-filing bug reports for automated workflows, keeping track of recurring
problems, and reducing the redundant data in the database.

Our proposed solution can flag a problem report as a duplicate before it
enters the system, which helps reduce operational expenses significantly. It
also allows customers to obtain a quicker resolution if the problem has
already been reported internally or by third-party stakeholders. Along with
preventing duplicates from entering the system [4], the proposed solution can
also help engineers debug or fix the issue efficiently by auto linking
existing identical bug reports to new bug reports.

###  I-A Terminology

A Bug Report can also be referred to as Problem Report or Defect Report in
this and other related studies. A Duplicate or Child bug report means that it
has been identified as a clone or duplicate of a previously submitted bug
report. The previously submitted report in the bug tracking system can be
referred to as a Master, Parent, or Original. Child and Parent reports are
mapped to each other. The rest of the reports that have never been linked as
Child or Parent to another bug report are called Unique Reports.

###  I-B Paper Outline

The rest of the paper is outlined as follows, Section 2 provides related work,
and Section 3 describes the methodology in detail. Section 4 reports the
evaluation results. Section 5 discusses the limitations of this approach.
Section 6 concludes the paper.

##  II Relevant Work

Runeson, Alexanderson and Nyholm conducted one of the earliest experiments
applying ML (Machine Learning) techniques to resolve the issue of duplicate
bug reports. For preprocessing, they used stemming, thesaurus, and
spellchecker and removed stop words. This resulted in a solution that can deal
with varying descriptions of the same bug report. After the preprocessing
steps, bug reports were transformed into vectors, and similarity computation
was done on these vectors to find duplicates. They reported a recall rate of
30% for k=5 and a search space of bug reports with age $<100$ days. [2].

In reference [5], Neysiani and Babamir contrasted Information Retrieval and
Machine Learning based approaches. The study observed that Information
Retrieval methods tend to return a list of potential duplicate bug reports
based on the similarity between the existing and new bug reports. While the ML
classification-based approach usually trains on pairs of identical and non-
identical bug reports. The results for retrieval-based strategies are
traditionally calculated using the Mean Average Precision and Recall rate.
Whereas the test results for ML classification are calculated using accuracy,
precision, recall, f1-score, etc. The study also reports that while the
accuracy reported for classification-based models is high, they are rarely
used in real-world bug tracking systems. This is because it requires two bug
reports to predict whether input reports are duplicates of each other. On the
other hand, information retrieval-based models are more useful as they report
a list of existing originals that are identical to the new bug report.

A. T. Nguyen et al. [6] used a hybrid approach combining information retrieval
and topic modeling. The study reported an accuracy improvement of 20%. It used
LDA for topic modeling bug reports. Topic modeling for bug reports considers
the probability of a word belonging to a particular topic. This works well
with bug reports using different terms describing the same issue. By using
BM25F and T-Model, Textual and Topic similarity is measured to detect existing
duplicate bug reports.

Haruna et al. Fine-tuned a pre-trained BERT triplet network that uses pairs.
Duplicate-Original bug reports made the positives class, and Original-
Unrelated bug reports made the negative class. The original bug report was the
same (anchor) in both positive and negative cases [7].

J. Lerch and M. Mezini proposed using stack traces available from execution to
run a search for existing bug reports with similar stack trace information in
them. Thus, saving the user the time to go through the hassle of filling out
the whole bug report only to find that it is a duplicate report [4].

Xiaoyin Wang et al. Used both the bug information like summary and description
along with the execution information to calculate similarities between two bug
reports. Using execution information is a interesting approach. They try to
replicate the bug using steps to reproduce and track the functions which get
called during the execution. They use this information to build a second set
of vectors and aggregate similarity scores from bug report information as well
as execution information. Their experiment achieved a 90% recall rate for k=5
on an extensively small dataset of only 232 bug reports with only 42 pairs of
duplicates [9].

Ashima Kukkar et al. Use CNN for feature extraction, a deep learning model
that captures the syntactic and semantic meaning of features used for
similarity computations. The study utilizes several openly available bug data
sources, for example, Eclipse, Open Office, Firefox, etc. The paper reports
89%,83%, and 78% recall rates for k=5 for Eclipse, Open Office, and Firefox,
respectively. The formula for the recall is based on the classification
approach. It also accounts for the correct classification of non-duplicate
reports, which comprise 80% to 90% of the data [10].

Alipour, Hindle and Stroulia, Incorporated more features by inference from bug
reports apart from the information already available. They feature engineered
contextual features to improve duplicate bug retrieval. They created 6
contexts, each linked to a list of most used words to describe the context.
When computing similarity between two bug reports, they use vectors from these
new features along with the typically available details from a bug report
[11].

B. Kucuk and E. Tuzun characterized Master, Duplicate, and Unique bug reports.
Also identified traits of resolved master-duplicate bug reports and unresolved
master-duplicate bug reports and proposed actions that may be taken for each
category [12]. The proposed actions seem pragmatic as they suggest re-opening
the resolved existing bug report or creating a new scope if the bug tracking
system allows it. For the Master bug report that is still open, the new bug
report can be prevented from being submitted. They also identified differences
in Parent and Child bug reports in the above two categories based on
characteristics such as severity, users involved, bug surface time, and life
cycle time.

Rocha and Carvalho, proposed SiameseQAT, which considers reports for
individual bugs as well as information from clusters of Parent-Child Reports
where each parent can have multiple children. Thus, the Parent Report becomes
the centroid of the collection. SiameseQAT combines syntactical and semantical
learning on structured, unstructured, and topic extraction-based features as
well. They reported results on Eclipse, Netbeans, and Open Office datasets and
reported a recall@25 mean of 85% [13].

H. Mahfoodh and M. Hammad used Word2Vec and Tensorflow to build a neural
network-based machine learning model to identify duplicate bug reports. They
use a word matching algorithm to determine the similarities between two
sentences [14].

In reference [15] Xiao, Du, Sui and Yue, introduced a novel heterogeneous
information network is used to incorporate structured details in addition to
unstructured characteristics of a bug report, for example, version, priority,
severity, etc. The HIN combines representation learning to develop a deep
neural network, hence HINDBR. They report a high 98% classification accuracy
on a small dataset made of pairs of duplicate and non-duplicate bug reports
[15].

##  III Our Approach

In this research, we have built a similarity search search mechanism by
leveraging NLP approaches that interact with each other to justify the extent
to which a Bug Report is similar to an existing one. Our system is made of
four components; Preprocessor, Data Transformer, Nominator, and Evaluator. As
outlined in Fig 1.

####  III-1 Preprocessor

As part of data preprocessing techniques, we first applies basic preprocessing
techniques to clean the data by tokenizing the unstructured text-based
information followed by some data domain relevant feature engineering steps
such as in our usecase, the alpha or alphanumeric are considered features. We
use the alphanumeric terms as these features can refer to vital information
that helps bug identification. For example, an alphanumeric script name is
provided as a part of a stack trace in the description of a bug report.
However, stack traces can also have process identification numbers or
hexadecimal memory addresses, which do not contribute specific information
that helps triage a bug report. Furthermore, it will try to find execution
failure traces in the description of the problem [4][9]. If found, the
preprocessor adds the extracted stack trace to the bug report summary again to
increase its weightage in learning [2]. Identical bug reports can have the
same values for specific characteristics like symptoms, impact, and priority
[12], considering these data is useful for the model [15]. If present these
characteristic data are concatenated to Summary. Not using them as independent
features helps reduce mismatches for reports where the values for these
attributes are not provided. The processed corpus or document is then fed to
the Data Transformer.

Figure 1: Proposed Solution Flowchart

####  III-2 Data Transformer

The Data Transformer component employs two methods to transform the text
documents into vectors. The first one is Scikit Learn’s Column Transformer
which acts as a wrapper for distinct types of data transformers. It allows the
aggregation of heterogeneous data as a single vector. For example, the column
transformer used here is made of TFIDF Vectorizer, Dictionary Vectorizer, and
a One Hot Encoder. Bug report attributes like Summary and Description are
transformed using TFIDF, while Component can be vectorized using OHE. The
Dictionary Vectorizer can be used for bug reports with structured data in a
dictionary format, for example, Platforms. The column transformer also
supports weight assignment. The Summary, Description, Component, and Platform
data transformers have weights of 45%, 25% [8], 25%, and 5%, respectively.
These transformations are used to build the Nearest Neighbors model. The
second method employed by Data Transformer simply concatenates the earlier
mentioned characteristics of the bug report and uses Input Example from the
Sentence Transformer library. These transformations are used in pairs to train
the BERT model.

####  III-3 Nominator

The Nominator utilizes unsupervised nearest neighbors learning. It can select
from three nearest neighbor algorithms: Brute Force, KD Tree, or Ball Tree.
The nominator algorithm attempts to determine the best approach from the
training data. The choice depends on factors like the number of samples, the
number of features, if the data is sparse or not, the number of neighbors
requested for each query, and the number of queries. If the number of samples
and searches is large, brute force would be a wrong choice, but if the amount
of data and inquiries is small, we can avoid the overhead of constructing the
KD or Ball Trees. The input to The Nominator component is the vectors from the
Column Transformer. For training, the Nominator would just determine the
nearest neighbor algorithm it needs to use and construct the Trees if
required. Which does not take more than a few minutes. We also maintain a
cache of recently submitted unique bug reports. This helps us detect identical
submissions in case a duplicate report is filed within a few minutes. Thus, we
ensure that model knowledge is current up to the latest bug report. For
Information Retrieval, the Nominator will return a list of nearest neighbors
along with the Euclidean distance from them. The similarity is then computed
by subtracting distance from the unity. This list of nearest neighbors
contains indices that can be used to get the Bug Report ID. We use these IDs
to load, preprocess and transform the nominated bug reports using the Input
Example transformer mentioned in the data transformer component. These vectors
are sent to the Evaluator for further assessment.

####  III-4 Evaluator

The evaluator uses BERT, a deep neural network, as a filter. BERT is trained
differently compared to unsupervised nearest neighbors learning used in the
Nominator. To train BERT we use pairs of known original-duplicate and
original-original bug reports, also known as the Siamese network. The
original-duplicate couple makes the positive case, while the original-original
couple is a negative case. For training BERT, we considered three times the
data samples of the positives as compared to the negatives. This helps reflect
the nature of data and the fact that fewer duplicate bug reports are filed
compared to unique bug reports. Furthermore, this allows BERT to learn a
distinct set of patterns from the Nearest Neighbors model. Internally, BERT
makes use of Transformer, an attention mechanism that learns contextual
relations between words (or sub-words) in a text. Based on this, the BERT
model evaluates predictions from unsupervised nearest neighbors. A New bug
report is paired with each recommendation from the model to assess whether
they are identical or not. Suppose the prediction is identified as a Negative
case in BERT, a new bug report paired with a similar bug report being an
original-original pair. In that case, it will suppress that recommendation
from the Nominator. This suppression applied by BERT eliminates False
Positives resulting in a shorter precise list of existing identical bug
reports.

##  IV Evaluation & Results

###  IV-A Dataset

  1. 1.

Private Dataset: A collection of bug reports from a private enterprise, It has
100+ components and 60K+ bug reports.

  2. 2.

Firefox Dataset: A collection of bug reports for the open-source web browser
Firefox. It has 50+ components and 100K+ bug reports.

  3. 3.

Eclipse Dataset: A collection of bug reports for the open-source Integrated
Development Environment. It has 20+ components and 80K+ bug reports.

###  IV-B Evaluation Metrics

For Information Retrieval systems, recall rate accuracy is typically important
performance evaluation metrics [3]. Recall rate is the percentage of
duplicates for which the original (also known as parent or master) bug report
is found in a list of size $n$ [1][2] expressed as below:

$Recall(n)={\sum_{i=1\ldots\\#duplicates}listed(n)\over\\#actual\ duplicates}$
(1)

The list presents identical bug reports that exist in the bug tracking system.
For a search of duplicates, where the size of the list of potentially
identical bug reports is minimal, i.e., k=5, we call it a successful search if
the Original or Parent bug report is present in the list. The recall rate [1]
thus helps us measure the efficiency of similarity search tools. We would also
like to see the parent report at or near the top of the list. The total
percentage of bug report from test data for which correct known parent report
is found in the list measures the accuracy of the model.

We use the Euclidean Distance to measure the similarity between two bug
reports. If the distance between two data samples is zero, the samples are
identical. Greater distance between two documents lowers the similarity.
Euclidean Distance calculation can be expressed as following:

$Distance(p,q)=\sqrt{\sum_{i=1}^{d}\left(q_{i}-p_{i}\right)^{2}}$ (2)

Where, $p,q$ are two points in Euclidean d-space, $q_{i},p_{i}$ are euclidean
vectors, starting from the origin of the space (initial point), and $d$ = d
dimensional space. Similarity can be calculated by subtracting the distance
from the unity for the data samples for which space is zero; the resemblance
is 1 or 100% in percentages expressed by (3) below:

$Similarity(p,q)=1-\sqrt{\sum_{i=1}^{d}\left(q_{i}-p_{i}\right)^{2}}$ (3)

We consider three significant data sources (private as well as public)
containing thousands of bug reports, one private and two public. Table I is
the count of bug reports used in training and testing the model. Experiments
indicate that the proposed solution achieves a high average accuracy of 70%
for k=5, i.e., recall@5 across all data. 0This means that the original is
found in a list recommending 5 similar bug reports that already exist in the
bug tracking system for 70% of duplicate bug reports. The reported results
could be conservative as the recall rate can be higher. This happens because
there are duplicate bug reports in the system which are allowed to exist as
Unique reports.

TABLE I: Training & Testing data count Dataset | Train Count | Test Count
---|---|---
Firefox | 90963 | 9332
Eclipse | 71131 | 6543
Private | 61396 | 4680

###  IV-C Results

The bug reports used are the ones that are known duplicates and each
corresponding Parent bug report. For each known Child bug report, we search
for a Parent bug report. We also keep track of a couple of things. First, if
the parent is present in the list or not. Second, the index of the parent if
found in the list, and lastly, the similarity percentage between the Child and
the Parent bug report. Following are the results per data set.

####  IV-C1 For Private Dataset

We observed that the Parent bug report was the top recommendation for 2332
Child reports. The distribution of position and count of bug reports,
similarity, and accuracy is represented by the charts Fig 2, Fig 3, Fig 4.

Figure 2: Identical Bug Report rang in Top 5 Recommendations, Private.

Figure 3: Distribution of Similarity for internal bug reports data.

Figure 4: 71.58% Parent BRs in Top 5 for internal data.

####  IV-C2 For Firefox Dataset

The Parent bug report was the top recommendation for this dataset for 3675
Child reports. The distribution of position and count of bug reports,
similarity, and accuracy is represented by the charts Fig 5, Fig 6, Fig 7.

Figure 5: Identical Bug Report rang in Top 5 Recommendations, Firefox.

Figure 6: Distribution of Similarity for Firefox bug reports data.

Figure 7: 67% Parent BRs in Top 5 for Firefox data.

####  IV-C3 For Eclipse Dataset

The Parent bug report was the top recommendation for 2466 Child reports. The
distribution of position and count of bug reports, similarity, and accuracy is
represented by the charts Fig 8, Fig 9, Fig 10.

Figure 8: Identical Bug Report rang in Top 5 Recommendations, Eclipse.

Figure 9: Distribution of Similarity for Eclipse bug reports data.

Figure 10: 71.14% Parent BRs in Top 5 for Eclipse data.

The overall performance of the proposed solution is depicted in the Table II.
The left column represents the position where the parent bug report was found.
The numbers in each row show the count of bug reports per data set found in
the corresponding place. -1 is the count of bug reports which had no parent in
the predicted list from model. Fig 11 represents Cumulative Count vs. Position
of Identified Parent Report.

TABLE II: Position and Count of Identified Parent Bug Reports Position | Private | Firefox | Eclipse
---|---|---|---
1 | 2332 | 3675 | 2466
2 | 482 | 901 | 535
3 | 219 | 469 | 295
4 | 157 | 328 | 209
5 | 158 | 277 | 159
-1 | 1329 | 2774 | 1486
Total | 4677 | 8424 | 5150

Figure 11: Cumulative Count vs. Position of Identified Parent Report

##  V Limitations

The results rely on the correctness of the data evaluated. We know that a bug
tracking system may not have all the duplicate bug reports identified and
marked, which restricts the performance calculation of the model. Also, some
child reports do not directly link to a parent report because they map to a
different child of the same parent report [13]. However, these discrepancies
only affect the score of the model, not the model itself. The model is able to
weed out the duplicate reports that existed in the system as originals. In
addition, the model links child reports to the correct parents as the Nearest
Neighbor’s model does not contain duplicate bug reports in its knowledge base.

##  VI Conclusion

Using deep neural network with BERT as the evaluator of Nearest Neighbors
Unsupervised learning, and an enhanced weighted heterogeneous data transformer
is a novel idea that helps us get high recall rate accuracy. Independently,
the nearest neighbor learning, which uses Term Frequency-based vector, fails
to achieve equivalent results [2][6]. It cannot learn the relationship among
the words in the Summary and Description of the bug report. This is because
TF-IDF does not consider the semantic meaning or context of the words[17],
whereas BERT does. On the other hand, BERT does not consider the all the
tokens in the document like TFIDF does, but BERT can derive the contexts of
the terms and understand their importance [16][7]. Since BERT uses deep neural
networks as part of its architecture [16], it can be computationally expensive
compared to Nearest Neighbors. BERT also limits the number of tokens it
considers to 512, which Nearest Neighbors does not. Put together, we overcome
the deficiencies of individual models. We get the best from each, Weighted
frequencies from TFIDF, fast search of Nearest Neighbors, and Contextual
understanding from BERT.

## References

  * [1] J. Natt och Dag, V. Gervasi, S. Brinkkemper and B. Regnell, ”Speeding up requirements management in a product software company: linking customer wishes to product requirements through linguistic engineering,” Proceedings. 12th IEEE International Requirements Engineering Conference, 2004., 2004, pp. 283-294, doi: 10.1109/ICRE.2004.1335685.
  * [2] P. Runeson, M. Alexandersson and O. Nyholm, ”Detection of Duplicate Defect Reports Using Natural Language Processing,” 29th International Conference on Software Engineering (ICSE’07), 2007, pp. 499-510, doi: 10.1109/ICSE.2007.32.
  * [3] Manning, C. D. and Schütze, H., Foundations of Statistical Natural Language Processing. Cambridge, USA: MIT Press, 1999.
  * [4] J. Lerch and M. Mezini, ”Finding Duplicates of Your Yet Unwritten Bug Report,” 2013 17th European Conference on Software Maintenance and Reengineering, 2013, pp. 69-78, doi: 10.1109/CSMR.2013.17.
  * [5] B. S. Neysiani and S. Morteza Babamir, ”Automatic Duplicate Bug Report Detection using Information Retrieval-based versus Machine Learning-based Approaches,” 2020 6th International Conference on Web Research (ICWR), 2020, pp. 288-293, doi: 10.1109/ICWR49608.2020.9122288.
  * [6] A. T. Nguyen, T. T. Nguyen, T. N. Nguyen, D. Lo and C. Sun, ”Duplicate bug report detection with a combination of information retrieval and topic modeling,” 2012 Proceedings of the 27th IEEE/ACM International Conference on Automated Software Engineering, 2012, pp. 70-79, doi: 10.1145/2351676.2351687.
  * [7] H. Isotani, H. Washizaki, Y. Fukazawa, T. Nomoto, S. Ouji and S. Saito, ”Duplicate Bug Report Detection by Using Sentence Embedding and Fine-tuning,” 2021 IEEE International Conference on Software Maintenance and Evolution (ICSME), 2021, pp. 535-544, doi: 10.1109/ICSME52107.2021.00054.
  * [8] T. Hirsch and B. Hofer, ”Identifying non-natural language artifacts in bug reports,” 2021 36th IEEE/ACM International Conference on Automated Software Engineering Workshops (ASEW), 2021, pp. 191-197, doi: 10.1109/ASEW52652.2021.00046.
  * [9] X. Wang, L. Zhang, T. Xie, J. Anvik and J. Sun, ”An approach to detecting duplicate bug reports using natural language and execution information”, Proc. 13th Int. Conf. Softw. Eng., pp. 461-470, 2008.
  * [10] A. Kukkar, R. Mohana, Y. Kumar, A. Nayyar, M. Bilal and K. -S. Kwak, ”Duplicate Bug Report Detection and Classification System Based on Deep Learning Technique,” in IEEE Access, vol. 8, pp. 200749-200763, 2020, doi: 10.1109/ACCESS.2020.3033045.
  * [11] A. Alipour, A. Hindle and E. Stroulia, ”A contextual approach towards more accurate duplicate bug report detection,” 2013 10th Working Conference on Mining Software Repositories (MSR), 2013, pp. 183-192, doi: 10.1109/MSR.2013.6624026.
  * [12] B. Kucuk and E. Tuzun, ”Characterizing Duplicate Bugs: An Empirical Analysis,” 2021 IEEE International Conference on Software Analysis, Evolution and Reengineering (SANER), 2021, pp. 661-668, doi: 10.1109/SANER50967.2021.00084.
  * [13] T. M. Rocha and A. L. D. C. Carvalho, ”SiameseQAT: A Semantic Context-Based Duplicate Bug Report Detection Using Replicated Cluster Information,” in IEEE Access, vol. 9, pp. 44610-44630, 2021, doi: 10.1109/ACCESS.2021.3066283.
  * [14] H. Mahfoodh and M. Hammad, ”Word2Vec Duplicate Bug Records Identification Prediction Using Tensorflow,” 2020 International Conference on Innovation and Intelligence for Informatics, Computing and Technologies (3ICT), 2020, pp. 1-6, doi: 10.1109/3ICT51146.2020.9311954.
  * [15] G. Xiao, X. Du, Y. Sui and T. Yue, ”HINDBR: Heterogeneous Information Network Based Duplicate Bug Report Prediction,” 2020 IEEE 31st International Symposium on Software Reliability Engineering (ISSRE), 2020, pp. 195-206, doi: 10.1109/ISSRE5003.2020.00027.
  * [16] Devlin, Jacob, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. “Bert: Pre-Training of Deep Bidirectional Transformers for Language Understanding.” arXiv.org, May 24, 2019.
  * [17] M. Li and B. -B. Yin, ”ARB-BERT: An Automatic Aging-Related Bug Report Classification Method based on BERT,” 2021 8th International Conference on Dependable Systems and Their Applications (DSA), 2021, pp. 474-483, doi: 10.1109/DSA52907.2021.00071.
  * [18] L. Wang, L. Zhang and J. Jiang, ”Detecting Duplicate Questions in Stack Overflow via Deep Learning Approaches,” 2019 26th Asia-Pacific Software Engineering Conference (APSEC), 2019, pp. 506-513, doi: 10.1109/APSEC48747.2019.00074.
  * [19] Z. Chen, X. Ju, Y. Shen and X. Chen, ”Improving Blocking Bug Pair Prediction via Hybrid Deep Learning,” 2021 IEEE 21st International Conference on Software Quality, Reliability and Security Companion (QRS-C), 2021, pp. 727-732, doi: 10.1109/QRS-C55045.2021.00110.
  * [20] S. K. Chauhan, A. Goel, P. Goel, A. Chauhan and M. K. Gurve, ”Research on product review analysis and spam review detection,” 2017 4th International Conference on Signal Processing and Integrated Networks (SPIN), 2017, pp. 390-393, doi: 10.1109/SPIN.2017.8049980.


# Edgeworth expansion for Bernoulli weighted mean

Pierre-Louis Cauvin
<EMAIL_ADDRESS>Criteo AI Lab.Grenoble INP - UGA,
Ensimag, 38400 Saint-Martin-d’Hères, France.

###### Abstract

In this work, we derive an Edgeworth expansion for the Bernoulli weighted mean
$\hat{\mu}=\frac{\sum_{i=1}^{n}Y_{i}T_{i}}{\sum_{i=1}^{n}T_{i}}$ in the case
where $Y_{1},\dots,Y_{n}$ are i.i.d. non semi-lattice random variables and
$T_{1},\dots,T_{n}$ are Bernoulli distributed random variables with parameter
$p$. We also define the notion of a semi-lattice distribution, which gives a
more geometrical equivalence to the classical Cramér’s condition in dimensions
bigger than 1. Our result provides a first step into the generalization of
classical Edgeworth expansion theorems for random vectors that contain both
semi-lattice and non semi-lattice variables, in order to prove consistency of
bootstrap methods in more realistic setups, for instance in the use case of
online AB testing.

##  1 Introduction and Main Result

In probability and statistic applied problems, one often have to deal with the
limit of sample distribution when its size gets bigger, in order to construct
asymptotically consistent estimators or confidence intervals for instance. The
most widely used theorem in this case is the classical Central Limit theorem
(CLT) and its variation the Delta method, which guarantee that a lot of
statistics converge in law to a Gaussian distribution. However, the CLT only
gives a convergence result but not its speed of convergence nor the error
terms of the approximation when we only deal with low finite sample size. The
first point can be tackme with Berry-Esseen bounds, stating that under some
regularity and finite moments conditions the convergence is of order
$\mathcal{O}(1/\sqrt{n})$ where $n$ is the sample size. Whereas for the second
point we need a precise analysis of error terms of the approximation, what
motivates the theory of Edgeworth expansions [7].

We say that a statistics $T(Y_{1},\dots,Y_{n})$ of some sample
$(Y_{1},\dots,Y_{n})$ admits an Edgeworth expansion of order $q$ if there
exist polynomials $p_{1},\dots,p_{q}$ whose coefficients depend only on the
moments of $(Y_{1},\dots,Y_{n})$ such that

$\displaystyle\mathbb{P}\left(T\leq
x\right)=\Phi(x)+\sum_{j=1}^{q}n^{-j/2}p_{j}(x)\phi(x)+\mathcal{O}(n^{-(q+1)/2})$

uniformly in $x$, where $\Phi$ and $\phi$ are the cumulative distribution
function and the probability distribution function of a standard Gaussian
distribution respectively.

The theory of Edgeworth expansions was vastly explored since the first formal
asymptotic expansions were proposed in a paper by Tchebycheff in 1890 [10],
and later on at the end of the $20^{th}$ century when contemporary
statisticians found its usefulness to study properties of empirical
statistical methods. One pioneer of this revival of interest was Peter Hall,
who single-handedly proved consistency, orders of convergence and coverage
error of confidence intervals for the B. Efron’s bootstrap [6], one of the
most used statistical approximation method, based on general existence
properties of Edgeworth expansions [7].

However, most of these general existence theorems hold only for function of
empirical means of i.i.d. samples and under what is called the Cramér’s
condition, which denies the use of semi-lattice (see Definition 2.6) random
variables, including many of the classical discrete distributions. Some
progress have been made since the start of the years 2000s in order to
generalize and weaken such conditions, for instance [8] shows the existence of
Edgeworth expansions for a broader class of statistics which contains $U$\-
and $L$-statistics, but still need the Cramér’s condition, whereas some other
authors have proposed some Edgeworth expansions for purely lattice triangular
array [2], bounded integer valued random variables [5], or have given a weaker
Cramér’s condition to allow discrete but non semi-lattice distribution [1].

Therefore, in this work we propose a first step into the generalization of
Edgeworth expansions for some semi-lattice random variables, which are of
vastly practical interests. More precisely, we prove the following theorem
which provides an Edgeworth expansion for the Bernoulli weighted mean of non
semi-lattice random variables, which is purely semi-lattice due to the mixture
of Bernoulli and non semi-lattice random variables as we will see in Section
2.

######  Theorem 1.1 (Edgeworth expansion for Bernoulli weighted mean).

Suppose that $Y_{1},\dots,Y_{n}$ are i.i.d. non semi-lattice random variables
with mean $\mu$ and variance $\sigma^{2}$. Let $T_{1},\dots,T_{n}$ be i.i.d.
Bernoulli distributed random variables with parameter $p>0$, and write
$N=\sum_{i=1}^{n}T_{i}$. Then, if $Y$ has at least $q$ finite moments, we have
the Edgeworth expansion

$\displaystyle\mathbb{P}\left(\sqrt{N}\frac{\hat{\mu}-\mu}{\sigma}\leq
x\right)=\Phi(x)+\sum_{j=1}^{q-2}n^{-j/2}p^{*}_{j}(x)\phi(x)+\mathcal{O}\left(n^{-(q-1)/2}\right)$
(1)

uniformly in $x$, where
$\hat{\mu}=\begin{cases}\frac{\sum_{i=1}^{n}Y_{i}T_{i}}{\sum_{i=1}^{n}T_{i}}&\text{if
$\sum_{i=1}^{n}T_{i}\neq 0$}\\\ \hfil 0&\text{otherwise}\end{cases}$ is the
Bernoulli weighted mean, and $p^{*}_{j}$’s are polynomials of degree at most
$3j$ depending only on the first $j+2$ moments of $Y_{1}$ and on moments of
$T$.

Bernoulli weighted means are the basis of many statistical online estimators,
such as uplift modeling in AB testing with various applications in health or
online advertising, where we want to estimate the different behaviors between
a treatment group and a control group, where each individual is assigned at
random based on a Bernoulli random variable. As such, we obtain estimators
which are difference or ratios of Bernoulli weighted means, motivating our
work. Another application of Bernoulli weighted means appears when we have a
big sample, e.g. the whole population, and we want to estimate statistical
properties of smaller samples taken at random with a Bernoulli variable
flipped for each individual.

This paper is organized as follow : in Section 2, we first state preliminary
definitions and results that will be useful later on. More precisely, we give
the classical Edgeworth expansion theorem [7], and we provide an equivalence
for the multivariate Cramér’s condition as the definition of a semi-lattice
random variable. We further use this equivalence to show that the classical
theorem doesn’t hold for Bernoulli weighted means. Then, in Section 3 we
provide a proof for Theorem 1.1, based on intermediate results using a
conditioning with respect to the Bernoulli random variables
$T_{1},\dots,T_{n}$ and asymptotic expansions for inverse binomial moments. We
also give an easier proof for the bound of the expectation of arbitrary powers
of a binomial random variable than the ones that can be found in [4] and [11].

##  2 Preliminary results

We will extensively use the notion of uniformly bounded equality, which we
recall its definition below, and we prove that it is stable with respect to
addition.

######  Definition 2.1 (Uniformly bounded).

For two sequences of functions
$(f_{n})_{n\in\mathbb{N}},(g_{n})_{n\in\mathbb{N}}:\mathbb{R}\to\mathbb{R}$
and a sequence $(u_{n})_{n\in\mathbb{N}}$ of $n$, we say that
$f_{n}(x)=g_{n}(x)+\mathcal{O}(u_{n})$ uniformly in $x$ if there exists an
integer $N>0$ such that

$\displaystyle\underset{x\in\mathbb{R}}{\sup}\,|f_{n}(x)-g_{n}(x)|\leq Cu_{n}$

for all $n\geq N$, where $C>0$ is a constant which does not depend on $x$.

######  Proposition 2.2 (Linearity of uniform bounds).

Let
$(f^{1}_{n})_{n\in\mathbb{N}},(f^{2}_{n})_{n\in\mathbb{N}},(g^{1}_{n})_{n\in\mathbb{N}},(g^{2}_{n})_{n\in\mathbb{N}}:\mathbb{R}\to\mathbb{R}$
be sequences of functions, and let
$(u^{1}_{n})_{n\in\mathbb{N}},(u^{2}_{n})_{n\in\mathbb{N}}$ be sequences of
$n$. Suppose that $f^{1}_{n}(x)=g^{1}_{n}(x)+\mathcal{O}(u^{1}_{n})$ and
$f^{2}_{n}(x)=g^{2}_{n}(x)+O(u^{2}_{n})$ both uniformly in $x$. Then

$\displaystyle
f^{1}_{n}(x)+f^{2}_{n}(x)=g^{1}_{n}(x)+g^{2}_{n}(x)+\mathcal{O}(u^{1}_{n}+u^{2}_{n})$

uniformly in $x$.

###### Proof.

By triangular inequality, we have

$\displaystyle\underset{x\in\mathbb{R}}{\sup}\,|f^{1}_{n}(x)+f^{2}_{n}(x)-g^{1}_{n}(x)-g^{2}_{n}(x)|$
$\displaystyle\leq\underset{x\in\mathbb{R}}{\sup}\,|f^{1}_{n}(x)-g^{1}_{n}(x)|+\underset{x\in\mathbb{R}}{\sup}\,|f^{2}_{n}(x)-g^{2}_{n}(x)|$
$\displaystyle\leq C^{1}u^{1}_{n}+C^{2}u^{2}_{n}$ $\displaystyle\leq
C(u^{1}_{n}+u^{2}_{n})$

where $C=\max(C^{1},C^{2})$. ∎

The following classical theorem on Edgeworth expansion for functions of sample
mean is due to [7] and can also be found in [9].

######  Theorem 2.3.

Let $(X_{n})$ be a sequence of i.i.d. random $d$-vectors and let $g$ be a
Borel measurable function on $\mathbb{R}^{d}$ taking real values. Suppose that
the following conditions hold

  1. 1.

$X_{1}$ has at least $q\geq 3$ finite moments,

  2. 2.

$g$ is $\mathcal{C}^{p-1}$ in a neighborhood of
$\mu=\mathbb{E}\left[\mspace{1.0mu}X_{1}\mspace{1.0mu}\right]$,

  3. 3.

$\nabla g(\mu)\neq 0$ where $\nabla$ denotes the gradient,

  4. 4.

(Cramér’s condition)
$\underset{||\mathbf{t}||\to\infty}{\limsup}~{}|\psi_{X_{1}}(\mathbf{t})|<1$
where
$\psi_{X_{1}}(\mathbf{t})=\mathbb{E}\left[\mspace{1.0mu}\exp(i\mathbf{t}\cdot
X_{1})\mspace{1.0mu}\right]$ is the characteristic function of $X_{1}$.

Then, we have the Edgeworth expansion

$\displaystyle\mathbb{P}\left(\sqrt{n}\,\dfrac{g\left(\bar{X}_{n}\right)-g(\mu)}{\sigma}\leq
x\right)=\Phi(x)+\sum_{j=1}^{q-2}n^{-j/2}p_{j}(x)\phi(x)+\mathcal{O}\left(n^{-(q-1)/2}\right)$

uniformly in $x$, where $\sigma^{2}=\nabla g(\mu)^{T}Var(X)\nabla g(\mu)$ and
$p_{j}$ is a polynomial of degree at most $3j$ whose coefficients depend only
on the first $j+2$ moments of $X_{1}$ and the partial derivatives of $g$ at
$\mu$.

This theorem is a powerful tool to prove existence of Edgeworth expansions for
a variety of different statistics. However, one of its hypothesis if quite
restrictive in many applied problems : Cramér’s condition. Cramér’s condition
is strongly linked with the notion of lattice distributions, that we define
below.

######  Definition 2.4 (Lattice distribution).

A random $d$-vector $X$ is said to have a lattice distribution if it takes
values almost surely in a lattice
$\\{\mathbf{x_{0}}+k\boldsymbol{\delta},~{}k\in\mathbb{Z}\\}$ spanned by
$\mathbf{x_{0}},\boldsymbol{\delta}\in\mathbb{R}^{d}$.

######  Example 2.5.

  1. 1.

The Bernoulli distribution in $d=1$ is a lattice distribution, because it
takes values in $\\{0,1\\}$ which is included in the lattice $\\{0+k\times
1,k\in\mathbb{Z}\\}$.

  2. 2.

If we take the uniform distribution on the Boolean hypercube $\\{0,1\\}^{d}$,
then it is a lattice distribution of dimension $d$.

  3. 3.

All absolutely continuous distributions with respect to Lebesgue measure are
not lattice.

For $d=1$, the Cramér’s condition is equivalent to the fact that $X$ has a non
lattice distribution. For instance, Cramér’s condition is always true when $X$
is absolutely continuous with respect to Lebesque measure. However, in the
case where $d>1$ we no longer have the equivalence between Cramér’s condition
and classical lattice distributions.

This is what motivates the following definition of semi-lattice distribution,
which is equal to lattice distribution for $d=1$, and that we will show is
also equivalent to Cramér’s condition in any finite dimension. This definition
also allows for a more geometrical view of multivariate Cramér’s condition :
there is no direction such that the distribution is lattice.

######  Definition 2.6 (Semi-lattice distribution).

A random $d$-vector $X$ is said to have a semi-lattice distribution if there
exists a vector $\mathbf{t^{*}}\in\mathbb{R}^{d}\setminus\\{0\\}$ such that
$\mathbf{t^{*}}\cdot X$ takes values almost surely in a lattice
$\\{x_{0}+k\delta,~{}k\in\mathbb{Z}\\}$ spanned by
$x_{0},\delta\in\mathbb{R}$.

######  Remark 2.7.

  1. 1.

When $d=1$, all semi-lattice distributions are discrete, but not all discrete
distributions are semi-lattice (e.g. $X\in\\{e,3,\pi\\}$ uniformly).

  2. 2.

When $d>1$, semi-lattice distributions are not necessarily discrete. For
example, the random vector $(N,4N)$ where $N\sim\mathcal{N}(0,1)$ is semi-
lattice as $(-4,1)\cdot(N,4N)=0$.

Now we prove the equivalence between semi-lattice behavior and Cramér’s
condition. More precisely, Proposition 2.8 shows that all semi-lattice
distributions do not verify Cramér’s condition, and on the other hand a
consequence of Proposition 2.9 is that if the Cramér’s condition doesn’t hold,
then the distribution is necessary semi-lattice.

######  Proposition 2.8.

If a random $d$-vector $X$ has a semi-lattice distribution, then there exists
a vector $\mathbf{t^{*}}\in\mathbb{R}^{d}\setminus\\{0\\}$ such that
$r\mapsto|\psi_{X}(r\mathbf{t^{*}})|$ is periodic. In particular, it means
that we have
$\underset{||\mathbf{t}||\to\infty}{\limsup}~{}|\psi_{X}(\mathbf{t})|=1$.

###### Proof.

Suppose that we have a $\mathbf{t^{*}}\neq 0$ such that $\mathbf{t^{*}}\cdot
X\in\\{x_{0}+k\delta,~{}k\in\mathbb{Z}\\}$ almost surely. Then, denoting
$\mu_{X}$ the law of $X$, we can write the function
$r\mapsto|\psi_{X}(r\mathbf{t^{*}})|$ as

$\displaystyle|\psi_{X}(r\mathbf{t^{*}})|=\left|\int_{\mathbb{R}^{d}}e^{ir\mathbf{t^{*}}\cdot\mathbf{x}}d\mu_{X}(\mathbf{x})\right|$
$\displaystyle=\left|\sum_{k\in\mathbb{Z}}e^{ir(x_{0}+k\delta)}\mathbb{P}\left(\mathbf{t^{*}}\cdot
X=x_{0}+k\delta\right)\right|$
$\displaystyle=\left|\sum_{k\in\mathbb{Z}}e^{irk\delta}\mathbb{P}\left(\mathbf{t^{*}}\cdot
X=x_{0}+k\delta\right)\right|$

Let us fix $r^{*}\in\mathbb{R}$, and consider the particular case
$r_{n}=\frac{2n\pi}{\delta}+r^{*}$. Using the equality above, we find that

$\displaystyle|\psi_{X}(r_{n}\mathbf{t^{*}})|=\left|\sum_{k\in\mathbb{Z}}e^{2ink}e^{ir^{*}k\delta}\mathbb{P}\left(\mathbf{t^{*}}\cdot
X=x_{0}+k\delta\right)\right|$
$\displaystyle=\left|\sum_{k\in\mathbb{Z}}e^{ir^{*}k\delta}\mathbb{P}\left(\mathbf{t^{*}}\cdot
X=x_{0}+k\delta\right)\right|$ $\displaystyle=|\psi_{X}(r^{*}\mathbf{t^{*}})|$

Therefore the function $r\mapsto|\psi_{X}(r\mathbf{t^{*}})|$ is
$\frac{2\pi}{\delta}$-periodic. In particular, setting
$\mathbf{t_{n}}=\frac{2n\pi}{\delta}\mathbf{t^{*}}$, it means that
$\underset{||\mathbf{t}||\to\infty}{\limsup}~{}|\psi_{X}(\mathbf{t})|\geq\underset{n\to\infty}{\lim}~{}|\psi_{X}(\mathbf{t_{n}})|=1$
and so
$\underset{||\mathbf{t}||\to\infty}{\limsup}~{}|\psi_{X}(\mathbf{t})|=1$. ∎

######  Proposition 2.9.

Let $X$ be a random $d$-vector admitting a characteristic function $\phi_{X}$.
Suppose that there exists a $t^{*}\neq 0$ such that $|\phi_{X}(t^{*})|=1$.
Then the distribution of $X$ is semi-lattice.

###### Proof.

If $|\phi_{X}(t^{*})|=1$, it means that there exists a $\theta\in\mathbb{R}$
such that $\phi_{X}(t^{*})=e^{i\theta}$. Rewriting $\phi_{X}$ by its
expression and dividing both sides by $e^{i\theta}$, it implies that
$\mathbb{E}\left[\mspace{1.0mu}e^{i(t^{*}\cdot
X-\theta)}\mspace{1.0mu}\right]=1$. Therefore, if we use the decomposition of
$e^{it}$ as cosinus and sinus, we obtain easily that
$\mathbb{E}\left[\mspace{1.0mu}\cos(t^{*}\cdot
X-\theta)\mspace{1.0mu}\right]=1$. But $\cos(t^{*}\cdot x-\theta)-1\geq 0$ for
all vector $x\in\mathbb{R}^{d}$, thus this equality is possible if and only if
$\cos(t^{*}\cdot X-\theta)=1$ almost surely. And so we have necessarily
$t^{*}\cdot X=\theta+k2\pi$ for $k\in\mathbb{Z}$, meaning that $X$ is semi-
lattice. ∎

This equivalence between semi-lattice and Cramér’s condition gives a
convenient way to verify if we can apply Theorem 2.3 on some random vector
without having to manipulate its characteristic function. For instance,
suppose that we would like to use Theorem 2.3 to show existence of an
Edgeworth expansion for the Bernoulli weighted mean
$\hat{\mu}=\frac{\sum_{i=1}^{n}Y_{i}T_{i}}{\sum_{i=1}^{n}T_{i}}$. This
statistic can indeed be expressed as a function of samples means by noticing
that $\hat{\mu}=\overline{YT}/\overline{T}=g(\overline{YT},\overline{T})$
where $g:x,y\mapsto x/y$, so we are in the framework of the theorem. However,
according to the following Lemma 2.10, the distribution of $(YT,T)$ doesn’t
follow Cramér’s condition, hence we unfortunately cannot apply Theorem 2.3 for
Bernoulli weighted means.

######  Lemma 2.10.

Let $Y$ be an arbitrary random variable and $T$ a Bernoulli distributed random
variable with parameter $p$. Then
$\underset{||\mathbf{t}||\to\infty}{\limsup}~{}|\psi_{(YT,T)}(\mathbf{t})|=1$
where $\psi_{(YT,T)}$ is the characteristic function of the vector $(YT,T)$.

###### Proof.

Setting $t^{*}=(0,1)$, we obtain $t^{*}\cdot(YT,T)=T$ which takes values
almost surely in the lattice $\\{0,1\\}$, meaning that the vector $(YT,T)$ has
a semi-lattice distribution. As such, we can use Prop 2.8 to conclude that
$\underset{||\mathbf{t}||\to\infty}{\limsup}~{}|\psi_{(YT,T)}(\mathbf{t})|=1$.
∎

##  3 Proof of Theorem 1.1

We recall that we define the Bernoulli weighted average $\hat{\mu}$ as

$\displaystyle\hat{\mu}=\begin{cases}\dfrac{\sum_{i=1}^{n}Y_{i}T_{i}}{\sum_{i=1}^{n}T_{i}}&\text{if
$\sum_{i=1}^{n}T_{i}\neq 0$}\\\ \hfil 0&\text{otherwise}\end{cases}$

where $Y_{1},\dots,Y_{n}$ are i.i.d non semi-lattice random variables, and
$T_{1},\dots,T_{n}$ are Bernoulli random variables with parameter $p>0$.
Furthermore, we let $N=\sum T_{i}$ and
$Z=\sqrt{N}\frac{\hat{\mu}-\mu}{\sigma}$ where $\mu$ and $\sigma^{2}$ are the
mean and variance of $Y$ respectively. For any $k=1,\dots,n$, we also define
$\overline{Y}_{k}=\frac{1}{k}\sum_{i=1}^{k}$ the empirical mean of sample size
$k$.

Before proceeding through the proof of Theorem 1.1, we first need to state and
prove some intermediate useful results. The following Lemma allows us to
reduce the problem of expanding the whole c.d.f. of $Z$, which depend both on
the $Y_{i}$’s and the $T_{i}$’s, to a problem of approximation for each
$k=1,\dots,n$ the easier c.d.f.
$F_{Z_{k}}(x):=\mathbb{P}\left(\sqrt{k}\frac{\overline{Y}_{k}-\mu}{\sigma}\leq
x\right)$ involving only the $Y_{i}$’s.

######  Lemma 3.1.

For every real $\alpha$, the cumulative distribution function $F_{Z}$ of $Z$
can be written as

$\displaystyle
F_{Z}(x)=\sum_{k=1}^{n}\binom{n}{k}p^{k}(1-p)^{n-k}\mathbb{P}\left(\sqrt{k}\dfrac{\overline{Y}_{k}-\mu}{\sigma}\leq
x\right)+\mathcal{O}(n^{\alpha})$

uniformly in $x$.

###### Proof.

Let $q_{1},\dots,q_{n}\in\\{0,1\\}$ where not all $q_{i}$’s are $0$. Then, by
conditioning we have

$\displaystyle\mathbb{P}\left(Z\leq
x~{}|~{}T_{1}=q_{1},\dots,T_{n}=q_{n}\right)=\mathbb{P}\left(\dfrac{\sqrt{\sum
q_{i}}}{\sigma}\left[\frac{1}{\sum q_{i}}\sum_{q_{i}=1}Y_{i}-\mu\right]\leq
x\right)$ (2)

Now, the idea is to notice that $\sum_{i=1}^{\sum q_{i}}Y_{i}$ has the same
distribution as $\sum_{q_{i}=1}Y_{i}$ due to the fact that the $Y_{i}$’s are
independent and of the same law. Therefore, if we define $k=\sum q_{i}$, we
can rewrite (2) as

$\displaystyle\mathbb{P}\left(Z\leq
x~{}|~{}T_{1}=q_{1},\dots,T_{n}=q_{n}\right)=\mathbb{P}\left(\sqrt{k}\dfrac{\overline{Y}_{k}-\mu}{\sigma}\leq
x\right)$ (3)

Now, we observe that the right-hand side of (3) only depends on $k$, meaning
that the conditional c.d.f. is invariant by permutations of the $q_{i}$’s.
Thus, by defining $N=\sum T_{i}$, we obtain

$\displaystyle\mathbb{P}\left(Z\leq
x~{}|~{}N=k\right)=\mathbb{P}\left(\sqrt{k}\dfrac{\overline{Y}_{k}-\mu}{\sigma}\leq
x\right)$

In the case where $q_{1}=\dots=q_{n}=0$, we have instead
$\mathbb{P}\left(Z\leq x~{}|~{}N=0\right)=\mathbb{P}\left(0\leq
x\right)=\mathbb{I}\\{x\geq 0\\}$. Hence, by the law of total expectation, the
c.d.f. of $Z$ can be expressed as

$\displaystyle F_{Z}(x)$
$\displaystyle=\sum_{k=0}^{n}\mathbb{P}\left(N=k\right)\mathbb{P}\left(Z\leq
x~{}|~{}N=k\right)$
$\displaystyle=\sum_{k=1}^{n}\binom{n}{k}p^{k}(1-p)^{n-k}\mathbb{P}\left(\sqrt{k}\dfrac{\overline{Y}_{k}-\mu}{\sigma}\leq
x\right)+\mathbb{I}\\{x\geq 0\\}(1-p)^{n}$

where we have used that $N$ follows a binomial distribution with parameters
($n$, $p$). To obtain the result it only remains to show that for all real
$\alpha$, $\mathbb{I}\\{x\geq 0\\}(1-p)^{n}=\mathcal{O}(n^{\alpha})$ uniformly
in $x$, which is immediate because $|\mathbb{I}\\{x\geq 0\\}|\leq 1$ and
$(1-p)^{n}$ decays to $0$ faster than $n^{\alpha}$. ∎

Hence, now we can focus ourselves on the existence of Edgeworth expansions for
each c.d.f. $F_{Z_{k}}$, which is what we prove in the following Lemma.

######  Lemma 3.2.

Suppose that $Y$ has at least $q$ finite moments. Then, for all positive
integer $k$, $\overline{Y}_{k}$ admits the Edgeworth expansion

$\displaystyle\mathbb{P}\left(\sqrt{k}\dfrac{\overline{Y}_{k}-\mu}{\sigma}\leq
x\right)=\Phi(x)+\sum_{j=1}^{q-2}k^{-j/2}p_{j}(x)\phi(x)+\mathcal{O}\left(k^{-(q-1)/2}\right)$

uniformly in $x$, where $p_{j}$ is a polynomial of degree at most $3j$ and
depending only on the first $j-2$ moments of $Y$.

###### Proof.

This a special case of Theorem 2.3 where $g=x\mapsto x$. Another direct proof
of this result is conducted in [7, Chapter 2.2]. ∎

Putting together this Edgeworth expansion result with the expression of Lemma
3.1, we will find ourselves with sums of the form
$\sum_{k=1}^{n}k^{\alpha}\binom{n}{k}p^{k}(1-p)^{k}$ with $\alpha$ real, that
we call the Bernoulli sum of $k^{\alpha}$. In what follows, we will provide
bounds for these Bernoulli sums in order to control exactly the error terms of
the asymptotic expansions.

First, we prove a simple bound on the expectation of powers of binomial random
variables, result that was already shown in [4] and [11], but we provide an
easier and shorter proof based on the Kullback-Leibler upper-bound for
binomial random variables.

######  Proposition 3.3.

Let $N$ be a binomial random variable with parameters $(n,p)$ where
$p\in(0,1/2]$. Then, for all $\alpha\in\mathbb{R}$ we have

$\displaystyle\mathbb{E}\left[\mspace{1.0mu}(N+1)^{\alpha}\mspace{1.0mu}\right]=\mathcal{O}(n^{\alpha})$

###### Proof.

If $\alpha\geq 0$ We have $N\leq n$, so
$\mathbb{E}\left[\mspace{1.0mu}(N+1)^{\alpha}\mspace{1.0mu}\right]\leq(n+1)^{\alpha}$
and we obtain immediately the result by noticing that
$(n+1)^{\alpha}=\mathcal{O}(n^{\alpha})$.
If $\alpha=-\beta<0$ In order to prove the result, we will use the following
classical inequality obtained by optimizing the Chernoff bound for a binomial
distribution : for all $\delta\in(0,p)$, we have

$\displaystyle\mathbb{P}\left(N\leq\delta n\right)\leq
e^{-n\mathcal{D}(\delta||p)}$ (4)

where $\mathcal{D}$ is the Kullback-Leibler divergence between two Bernoulli
distributions of parameters $\delta$ and $p$ respectively. Turning back to our
former expectation, by conditioning on the event $\\{N\leq\delta n\\}$ we
obtain

$\displaystyle\mathbb{E}\left[\mspace{1.0mu}(N+1)^{\alpha}\mspace{1.0mu}\right]=\mathbb{P}\left(N\leq\delta
n\right)\mathbb{E}\left[\mspace{1.0mu}\dfrac{1}{(N+1)^{\beta}}~{}|~{}N\leq\delta
n\mspace{1.0mu}\right]+\mathbb{P}\left(N>\delta
n\right)\mathbb{E}\left[\mspace{1.0mu}\dfrac{1}{(N+1)^{\beta}}~{}|~{}N>\delta
n\mspace{1.0mu}\right]$

Now, we can use (4) to bound $\mathbb{P}\left(N\leq\delta n\right)$, and we
also have
$\mathbb{E}\left[\mspace{1.0mu}\frac{1}{(N+1)^{\beta}}~{}|~{}N\leq\delta
n\mspace{1.0mu}\right]\leq 1$ due to $N+1\geq 1$, $\mathbb{P}\left(N>\delta
n\right)\leq 1$ because it’s a probability, and
$\mathbb{E}\left[\mspace{1.0mu}\frac{1}{(N+1)^{\beta}}~{}|~{}N>\delta
n\mspace{1.0mu}\right]\leq\frac{1}{(\delta n+1)^{\beta}}$. All in all, we
finaly obtain

$\displaystyle\mathbb{E}\left[\mspace{1.0mu}(N+1)^{\alpha}\mspace{1.0mu}\right]\leq
e^{-n\mathcal{D}(\delta||p)}+(\delta n+1)^{\alpha}$

where the right-hand side is a $\mathcal{O}(n^{\alpha})$ because the Kullback-
Leibler divergence is strictly positive for $\delta<p$, proving the
proposition. ∎

An immediate consequence of this bound is the following Lemma, controlling the
speed of convergence of the Bernoulli sum of $k^{\alpha}$.

######  Lemma 3.4.

Let $p\in(0,1/2]$ and $n\geq 1$. Then, for all $\alpha\in\mathbb{R}$

$\displaystyle\sum_{k=1}^{n}k^{\alpha}\binom{n}{k}p^{k}(1-p)^{n-k}\,\overset{(1)}{=}\,np\mathbb{E}\left[\mspace{1.0mu}(N^{*}+1)^{\alpha-1}\mspace{1.0mu}\right]\,\overset{(2)}{=}\,\mathcal{O}(n^{\alpha})$

where $N^{*}$ is a binomial random variable with parameters $(n-1,p)$.

###### Proof.

Expanding the expectation and using a change of index, we obtain

$\displaystyle\mathbb{E}\left[\mspace{1.0mu}(N^{*}+1)^{\alpha-1}\mspace{1.0mu}\right]$
$\displaystyle=\sum_{k=0}^{n-1}(k+1)^{\alpha-1}\binom{n-1}{k}p^{k}(1-p)^{n-1-k}$
$\displaystyle=\sum_{k=1}^{n}k^{\alpha-1}\binom{n-1}{k-1}p^{k-1}(1-p)^{n-k}$
$\displaystyle=\sum_{k=1}^{n}k^{\alpha-1}\dfrac{k}{n}\binom{n}{k}\dfrac{1}{p}p^{k}(1-p)^{n-k}$
$\displaystyle=\dfrac{1}{np}\sum_{k=1}^{n}k^{\alpha}\binom{n}{k}p^{k}(1-p)^{n-k}$

which proves equality $(1)$. Now, by Proposition 3.3 it is straightforward
that
$\mathbb{E}\left[\mspace{1.0mu}(N^{*}+1)^{\alpha-1}\mspace{1.0mu}\right]=\mathcal{O}\left((n-1)^{\alpha-1}\right)=\mathcal{O}(n^{\alpha-1})$,
and so equality $(2)$ holds immediately by multiplying both side by $np$. ∎

However, in order to prove the Edgeworth expansion of Theorem 1.1 we will also
need an even stronger result on binomial random variables, namely an
asymptotic expansion for the inverse binomial moments.

######  Definition 3.5 (Inverse binomial moment).

For any real $\alpha>0$, we define the inverse binomial moment $f_{\alpha}(n)$
of order $\alpha$ as

$\displaystyle
f_{\alpha}(n)=\sum_{k=1}^{n}\dfrac{1}{k^{\alpha}}\binom{n}{k}p^{k}(1-p)^{n-k}$

The following result was first stated and proved in [11] with a general closed
form for the constants of the expansion.

######  Proposition 3.6.

The inverse binomial moment $f_{\alpha}(n)$ can be expanded in terms of
inverse powers of $(np)$, meaning that for any integer $q$

$\displaystyle
f_{\alpha}(n)=\dfrac{1}{(np)^{\alpha}}\left(\sum_{k=0}^{q-1}C_{\alpha,k}\dfrac{1}{(np)^{k}}\right)+\mathcal{O}\left(\dfrac{1}{(np)^{\alpha+q}}\right)$

where $C_{\alpha,k}$ are constants depending only on $\alpha$ and $p$.

###### Proof.

See [11] for the original proof using expansion of the generating function of
central moments, or [3] for another generalized approach based on Bernstein
polynomials. ∎

Finally, we need this last Lemma in order to construct asymptotic expansions
for sums of inverse binomial moments.

######  Lemma 3.7.

For every $j=1,\dots,q-2$, let $p_{j}$ be a polynomial of degree at most $3j$
and depending only on the first $j+2$ moments of some r.v. $Y$. Furthermore,
let $\phi$ be the p.d.f. of the standard Gaussian distribution. Then, we have

$\displaystyle\sum_{j=1}^{q-2}f_{j/2}(n)p_{j}(x)\phi(x)=\sum_{j=1}^{q-2}n^{-j/2}p^{*}_{j}(x)\phi(x)+\mathcal{O}\left(n^{-(q-1)/2}\right)$

uniformly in $x$, where $p^{*}_{j}$’s are polynomials of degree at most $3j$
and depending only on the first $j+2$ moments of $Y$ and on moments of $T$.

###### Proof.

For each $j=1,\dots,q-2$, we can use Proposition 3.6 to obtain

$\displaystyle\left|f_{j/2}(n)-\sum_{i=0}^{q-2+\left\lceil
j/2\right\rceil}C_{j,i}(np)^{-j/2+i}\right|\leq K_{j}n^{-(q-1)/2}$

where $K_{j}>0$ is a constant. This comes from the fact that
$-j/2+q-2+\left\lceil j/2\right\rceil<q-1$ and $-j/2+q-1+\left\lceil
j/2\right\rceil\geq q-1$. Therefore, for all $x\in\mathbb{R}$, we have also
that

$\displaystyle\left|f_{j/2}(n)p_{j}(x)\phi(x)-\sum_{i=0}^{q-2+\left\lceil
j/2\right\rceil}C_{j,i}(np)^{-j/2+i}p_{j}(x)\phi(x)\right|\leq\left|p_{j}(x)\phi(x)\right|K_{j}n^{-(q-1)/2}.$

However $p_{j}(x)\phi(x)$ is bounded uniformly because $\phi(x)$ has an
exponential decay, and so we can put the bound inside the constant $K_{j}$.
Hence, for each $j=1,\dots,q-2$ we conclude that

$\displaystyle f_{j/2}(n)p_{j}(x)\phi(x)=\sum_{i=0}^{q-2+\left\lceil
j/2\right\rceil}C_{j,i}(np)^{-j/2+i}p_{j}(x)\phi(x)+\mathcal{O}\left(n^{-(q-1)/2}\right).$
(5)

uniformly in $x$. Now we can sum equality (5) over every $j$, giving

$\displaystyle\sum_{j=1}^{q-2}f_{j/2}(n)p_{j}(x)\phi(x)=\sum_{j=1}^{q-2}~{}\sum_{i=0}^{q-2+\left\lceil
j/2\right\rceil}C_{j,i}(np)^{-j/2+i}p_{j}(x)\phi(x)+\mathcal{O}\left(n^{-(q-1)/2}\right)$

still uniformly in $x$ due to Proposition 2.2. Noticing that for every couple
$(i,j)$ we have $-j/2+i=-k/2$ where $1\leq k\leq q-2$ and summing over each of
such $k$, we can therefore rewrite the equality as

$\displaystyle\sum_{j=1}^{q-2}f_{j/2}(n)p_{j}(x)\phi(x)=\sum_{j=1}^{q_{2}}n^{-j/2}p^{*}_{j}(x)\phi(x)+\mathcal{O}\left(n^{-(q-1)/2}\right)$

unformly in $x$, where $p^{*}_{j}$ is a polynomial given by

$\displaystyle
p^{*}_{j}(x)=\begin{cases}p^{-j/2}\sum_{i=0}^{\frac{j-1}{2}}C_{2i+1,\frac{j-1}{2}-i}~{}p_{2i+1}(x)&\text{if
$j$ odd}\\\
p^{-j/2}\sum_{i=1}^{\frac{j}{2}}C_{2,\frac{j}{2}-i}~{}p_{2i}(x)&\text{if $j$
even}\end{cases}$

With this expression, it’s easy to see that $p^{*}_{j}$ has the same degree as
$p_{j}$, that it only depends on the first $j+2$ moments of $Y$ because it is
a linear combination of all polynomials $p_{k}$ where $k\leq j$ is of same
parity as $j$, and that it also depends on the moments of $T$ because of the
constants $C_{j,i}$ depending on $p$. Thus proving the result. ∎

After all these intermediate results, we are finally ready to prove Theorem
1.1.

###### Proof of Theorem 1.1.

Using the Edgeworth expansions of Lemma 3.2 in each term of the decomposition
of the c.d.f. $F_{Z}$ provided by Lemma 3.1, for $\alpha=-(q-1)/2$ we obtain

$\displaystyle F_{Z}(x)$
$\displaystyle=\left[1-(1-p)^{n}\right]\Phi(x)+\sum_{j=1}^{q-2}\left(\sum_{k=1}^{n}k^{-\frac{j}{2}}\binom{n}{k}p^{k}(1-p)^{n-k}\right)p_{j}(x)\phi(x)$
$\displaystyle+\mathcal{O}\left(\sum_{k=1}^{n}k^{-\frac{q-1}{2}}\binom{n}{k}p^{k}(1-p)^{n-k}\right)+\mathcal{O}(n^{-(q-1)/2})$
(6)

uniformly in $x$ due to Proposition 2.2. Here the
$\left[1-(1-p)^{n}\right]\Phi(x)$ term comes from the fact that
$\sum_{k=1}^{n}\binom{n}{k}p^{k}(1-p)^{n-k}=1-(1-p)^{n}$ since the sum starts
at $1$ and not at $0$. Now, the idea is to bound every term involving binomial
sums in order to find an Edgeworth expansion as powers of $n^{-1/2}$ instead
of powers of $k^{-1/2}$.
First, we have that $|(1-p)^{n}\Phi(x)|\leq(1-p)^{n}$ which is uniformly
bounded in $x$ by $n^{-(q-1)/2}$, and so we can put it inside the already
present $\mathcal{O}(n^{-(q-1)/2})$.
Next, according to Lemma 3.4 we have
$\sum_{k=1}^{n}k^{-\frac{q-1}{2}}\binom{n}{k}p^{k}(1-p)^{n-k}=\mathcal{O}(n^{-(q-1)/2})$,
hence allowing us to rewrite (3) as

$\displaystyle F_{Z}(x)$
$\displaystyle=\Phi(x)+\sum_{j=1}^{q-2}\left(\sum_{k=1}^{n}k^{-\frac{j}{2}}\binom{n}{k}p^{k}(1-p)^{n-k}\right)p_{j}(x)\phi(x)+\mathcal{O}(n^{-(q-1)/2})$
(7)

uniformly in $x$. Finally, it only remains to apply Lemma 3.7 in order to
bound uniformly the last sum term of (7), giving

$\displaystyle F_{Z}(x)$
$\displaystyle=\Phi(x)+\sum_{j=1}^{q-2}n^{-j/2}p^{*}_{j}(x)\phi(x)+\mathcal{O}(n^{-(q-1)/2})$

where $p^{*}_{j}$ are polynomials verifying the required conditions, thus
proving our main theorem. ∎

## Acknowledgments

This work was done as part of a 3-month internship in the Causality team of
Criteo AI Lab. The author would like to thank all members of this team for
their hospitality, especially Matthieu Martin for having supervised this
internship and reviewed parts of this paper several times.

## References

  * [1] J. Angst and G. Poly, “A weak Cramér condition and application to Edgeworth expansions,” Electronic Journal of Probability 22 no. 59, (2017) 1 – 24.
  * [2] A. Bock, “Edgeworth expansions for lattice triangular arrays.” 2014.
  * [3] J. Cichoń and Z. Golebiewski, “On bernoulli sums and bernstein polynomials,” Discrete Mathematics & Theoretical Computer Science DMTCS Proceedings vol. AQ, 23rd Intern. Meeting on Probabilistic, Combinatorial, and Asymptotic Methods for the Analysis of Algorithms (AofA’12) (Jan., 2012) 179–190.
  * [4] F. Cribari-Neto, N. L. Garcia, and K. L. P. Vasconcellos, “A note on inverse moments of binomial variates,” Brazilian Review of Econometrics 20 no. 2, (2000) .
  * [5] D. Dolgopyat and Y. Hafouta, “Edgeworth expansions for independent bounded integer valued random variables,” 2020.
  * [6] B. Efron, “Bootstrap methods: Another look at the jackknife,” The Annals of Statistics 7 no. 1, (1979) 1–26.
  * [7] P. Hall, The Bootstrap and Edgeworth Expansion.  Springer Series in Statistics. Springer New York, 1997.
  * [8] B.-Y. Jing and Q. Wang, “A unified approach to edgeworth expansions for a general class of statistics,” Statistica Sinica 20 no. 2, (2010) 613–636.
  * [9] J. Shao and D. Tu, The Jackknife and Bootstrap.  Springer Series in Statistics. Springer New York, 1995.
  * [10] P. Tchebycheff, “Sur deux théorèmes relatifs aux probabilités,” Acta Mathematica 14 no. none, (1890) 305 – 315.
  * [11] M. Znidaric, “Asymptotic expansion for inverse moments of binomial and poisson distributions,” The Open Statistics & Probability Journal 1 no. 1, (Jan, 2009) 7–10.


# Time Series Forecasting Based on Augmented Long Short-Term Memory

Daniel Hsu

###### Abstract

In this paper, we use recurrent autoencoder model to predict the time series
in single and multiple steps ahead. Previous prediction methods, such as
recurrent neural network (RNN) and deep belief network (DBN) models, cannot
learn long term dependencies. And conventional long short-term memory (LSTM)
model doesn’t remember recent inputs. Combining LSTM and autoencoder (AE), the
proposed model can capture long-term dependencies across data points and uses
features extracted from recent observations for augmenting LSTM at the same
time. Based on comprehensive experiments, we show that the proposed methods
significantly improves the state-of-art performance on chaotic time series
benchmark and also has better performance on real-world data. Both single-
output and multiple-output predictions are investigated.

##  1 Introduction

Time series forecasting and modeling is an important interdisciplinary field
of research, involving among others Computer Sciences, Statistics, and
Econometrics. Made popular by Box and Jenkins [1] in the 1970s, traditional
modeling procedures combine linear autoregression (AR) and moving average.
But, since data are nowadays abundantly available, often complex patterns that
are not linear can be extracted. So, the need for nonlinear forecasting
procedures arises.

Recently, neural networks with deep architectures have proven to be very
successful in image, video, audio and language leaning tasks [6]. In time
series forecasting area, though traditionally shallow neural networks are
generally adopted, the deep neural networks have also aroused enormous
interests among researchers. Deep belief networks (DBN) are frequently
employed in current short-term traffic forecasting [7][8], and pre-training
strategies with unsupervised learning algorithms such as Restricted Boltzmann
machine (RBM) [9] and Stacked AutoEncoder (SAE) [11] are also used. However,
these deep architectures can not capture the long dependencies across data
points which are beyond input observations.

RNNs are particularly suitable for modeling dynamical systems as they operate
on input information as well as a trace of previously acquired information
(due to recurrent connections) allowing for direct processing of temporal
dependencies. RNNs can be employed for a wide range of tasks as they inherit
their flexibility from plain neural networks. Among all RNN architectures, the
most successful one to characterize long-term memory is the long short-term
memory network (LSTM) [3], which learns both short-term and long-term memory
by enforcing constant error flow through the designed cell state.

Recent advances in deep learning, especially recurrent neural network (RNN)
and long short-term memory (LSTM) models [5][4][13], have been proposed to
tackle this problem. However, these models still have some disadvantages.
Especially LSTM model cannot work well on cases where the prediction is
primarily based on recent past observations [5].

In this work, we investigate time series forecasting by combining AutoEncoder
(AE) and LSTM. In this model, the input observations are first encoded to
latent variables, which is equivalent to feature extraction. In the LSTM cell,
the input observations and latent variables are incorporated into the state
transition. The decoder is a multi-layer neural network which maps hidden
state of LSTM and latent variables of input into the predicted values. We term
the proposed model as augmented long short-term memory (A-LSTM). The LSTM
cell, encoder and decoder are trained together by ADAM [20]. This model can
not only capture long-term dependencies, but also augment the prediction by
the extracted features from recent observations at the same time.

The experiments show that this model can not only improve the performance of
benchmarks on chaotic time series prediction, but also performs better on
real-world data than previous methods such as DBN pre-trained with RBM [9],
SAE [11] and LSTM models [13]. Moreover, with 5-step and 10-step experiments,
the performance on multi-output prediction task doesn’t degrade as the number
of output increases, showing the proposed method has strong potential for
scalability.

##  2 Preliminary

###  2.1 Autoencoder

The Autoencoder (AE) was first introduced as a dimension-reduction model [19].
An autoencoder takes an input vector $\bm{x}$ and transforms it into a latent
representation $\bm{z}$. The transformation, typically referred as encoder,
follows the following equation:

$\bm{z}=\sigma(\bm{W}\bm{x}+\bm{b})$

where $\bm{W}$ and $\bm{b}$ correspond to the weighs and bias in the encoder,
and $\sigma$ is the sigmoid function.

The resulting latent representation $\bm{z}$ is then mapped back into the
reconstructed feature space $\bm{y}$ in the decoder as follows:

$\bm{y}=\sigma(\bm{W}^{\prime}\bm{z}+\bm{b}^{\prime})$

where $\bm{W}^{\prime}$ and $\bm{b}^{\prime}$ correspond to the weighs and
bias in the decoder. The autoencoder is trained by minimizing the
reconstruction error $\|\bm{y}-\bm{x}\|$. Multiple autoencoders can be
connected to construct stacked autoencoder [11], which can be used to learn
multiple levels of non-linear features. In this work, we utilize AE to
extractor features from input observations, which can make LSTM learn more
efficiently.

###  2.2 Recurrent Neural Network

RNNs are discrete-time state–space models trainable by specialized weight
adaptation algorithms. The input to RNN is a variable-length sequence
$\bm{x}=(\bm{x}_{1},\ldots,\bm{x}_{T})$ which can be recursively processed.
And when processing each symbol, RNN maintains its internal hidden state
$\bm{h}$. The operation of RNN at each timestep $t$ can be formulated as

$\bm{h}_{t}=f_{\theta}(\bm{x}_{t},\bm{h}_{t-1})$

where $f$ is the deterministic state transition function and $\theta$ is the
parameter of $f$. The output of RNN is computed using the following equation:

$\bm{y}_{t}=g_{\varphi}(\bm{h}_{t})$ (1)

where function $g$ can be modeled as a neural network with weights $\varphi$.
In implementation, the function $f$ can be realized by long short-term memory
[3].

Although traditional RNN exhibits a superior capability of modeling nonlinear
time series problems in an effective fashion, there are still several issues
to be addressed [2]:

  * •

Traditional RNNs cannot train the time series with very long time lags, which
is commonly seen in real-world datasets.

  * •

Traditional RNNs rely on the predetermined time lags to learn the temporal
sequence processing, but it is difficult to find the optimal window size in an
automatic way.

###  2.3 Long Short-Term Memory

To overcome the aforementioned disadvantages of traditional RNNs, Long Short-
Term Memory (LSTM) neural network is proposed in this study to predict time
series in single-step and multi-step ahead. LSTM was initially introduced in
[3] with the objective of modeling long-term dependencies and determining the
optimal time lag for time series problems. A LSTM is composed of one input
layer, one recurrent hidden layer, and one output layer. The basic unit in the
hidden layer is memory block, containing memory cells with self-connections
memorizing the temporal state, and a pair of adaptive, multiplicative gating
units to control information flow in the block. It also has input gate and
output gate controlling the input and output activations into the block.

The memory cell is primarily a recurrently self-connected linear unit, called
Constant Error Carousel (CEC), and the cell state is represented by the
activation of the CEC. Because of CEC, the multiplicative gates can learn when
to open and close. Then by keeping the network error constant, the vanishing
gradient problem can be solved in LSTM. Moreover, a forget gate was added to
the memory cell, which can prevent the gradient from exploding when learning
long time series. The operation and structure of LSTM can be described as
below:

$\displaystyle\bm{i}_{t}$ $\displaystyle=$
$\displaystyle\sigma(\bm{W}_{i}\bm{x}_{t}+\bm{U}_{i}\bm{m}_{t-1}+\bm{b}_{i})$
$\displaystyle\bm{s}_{t}$ $\displaystyle=$
$\displaystyle\text{tanh}(\bm{W}_{s}\bm{x}_{t}+\bm{U}_{s}\bm{m}_{t-1}+\bm{b}_{s})$
$\displaystyle\bm{f}_{t}$ $\displaystyle=$
$\displaystyle\sigma(\bm{W}_{f}\bm{x}_{t}+\bm{U}_{f}\bm{m}_{t-1}+\bm{b}_{f})$
$\displaystyle\bm{y}_{t}$ $\displaystyle=$
$\displaystyle\sigma(\bm{W}_{o}\bm{x}_{t}+\bm{U}_{o}\bm{m}_{t-1}+\bm{b}_{o})$
$\displaystyle\bm{c}_{t}$ $\displaystyle=$
$\displaystyle\bm{c}_{t-1}\circ\bm{f}_{t}+\bm{s}_{t}\circ\bm{i}_{t}$
$\displaystyle\bm{m}_{t}$ $\displaystyle=$
$\displaystyle\bm{s}_{t}\circ\bm{o}_{t}$ (2)

where $\bm{i}_{t},\bm{f}_{t}$ and $\bm{y}_{t}$ are denoted as input gate,
forget gate and output gate at time $t$ respectively, and $\bm{m}_{t}$ and
$\bm{c}_{t}$ represent the hidden state and cell state of the memory cell at
time $t$.

##  3 Related Work and Motivation

Driven by the recent success of deep learning [6], several different Deep
Learning approaches can be found in the literature for performing time series
predictions. For example, deep belief networks are used in the work of [9]
along with restricted Boltzman machine (RBM). [10] also compares the
performance of Deep Belief Networks with that of Stacked Denoising
Autoencoders. This last type of network is also employed by [11] to predict
the temperature of an indoor environment. Another application of Time-Series
forecasting can be found in [8], which uses Stacked Autoencoders (SAE) to
predict the flow of traffic from a big data dataset. However, as compared in
[12], deep learning models such as RBM and SAE perform worse than LSTM because
they cannot capture long-term dependencies across data points, and fixed size
window size is another factor for suboptimal performance.

LSTM [3] is another learning structure often used recently in time series
prediction. [5] first used LSTM to predict chaotic time series. In [4], an
LSTM sequence-to-sequence model was used to predict next values. A survey [18]
reviews many applications of LSTM and other RNN architectures to short-term
load forecasting problem. However, the sequence-to-sequence RNN models [4],
aiming at predicting a sequence of future values on the basis of a sequence of
past values, can not handle very long sequences and model periods in time
series [15]. And every input sequence was padded to the same length. As
indicated in [5], LSTM model such as [13] cannot utilize recent observations
effectively, since it spends too much resources on long-term memory.

The time series forecasting with multi-output and multi-step is an active
research topic. [25] investigated this problem by employing multiple-output
support vector regression (M-SVR) with multiple-input multiple-output (MIMO)
prediction strategy. In [26], a cooperative neural evolution method for multi-
task learning was adopted, which enables neural networks to be trained with
shared knowledge representation.

Motivated by disadvantages of previous models above, we propose a LSTM model,
where the hidden state transition is dependent on both input feature extracted
by the autoencoder and the previous hidden state. The prediction output of the
proposed model is decoded from both hidden state and latent representation of
the input. So, our model can utilize long-term dependencies and recent
observations at same time. And it can work on long sequences with variable
lengths. We conduct performance comparison on both single-output and multi-
output forecasting experiments.

We notice that there are some recent work [16][17] which stack AE and RNN or
LSTM together to predict time series. Different from our work, the AE and LSTM
are trained separately in their models, and the decoder of AE is directly cut
off, which may yield suboptimal weights for AE and LSTM cell. And they don’t
have results on multi-output prediction. The proposed method can outperform
these models on chaotic and real-world datasets.

##  4 Augmented Long Short Term Memory Model

In this work, we incorporate autoencoder (AE) into LSTM model to improve time
series forecasting performance. Here, $\bm{x}_{t},\bm{z}_{t}$ and $\bm{h}_{t}$
denote the input observations, latent variables of the input and the hidden
state of LSTM at time $t$ respectively. The hidden state
$\bm{h}_{t}\triangleq(\bm{c}_{t},\bm{m}_{t})$ includes both the cell state
$\bm{c}_{t}$ and hidden state $\bm{m}_{t}$ of the memory cell as described in
(2). In the proposed model, the LSTM is augmented by the latent variables
learned from AE and corresponding feature extractors.

###  4.1 Encoder Model

Different from conventional autoencoder in [17], the latent variable extracted
at time $t$ is dependent on both the input observations and the previous
hidden state of LSTM cell. It can be described as below:

$\bm{z}_{t}=\varphi^{\text{enc}}_{\tau}(\varphi^{\bm{x}}(\bm{x}_{t}),\bm{h}_{t-1})$
(3)

where the decoder $\varphi^{\text{enc}}_{\tau}$ and feature extractor
$\varphi^{\bm{x}}(\bm{x}_{t})$ are both implemented as a two-layer neural
network. Here, the input $\bm{x}_{t}$ is not directly processed by the encoder
together with hidden state. The encoder only read the extracted feature of
input $\bm{x}_{t}$, which can reduce some redundant information in the input
observations. It is especially important when processing time series with low
sampling rate.

###  4.2 Decoder Model

This model applies AE in a recurrent setting. At each timestep, different from
conventional AE, the output of the decoder, i.e. predicted value, is dependent
on both the current latent variables of input and the previous hidden state of
LSTM. This process can be depicted as:

$\hat{\bm{y}}_{t}=\varphi^{\text{dec}}_{\tau}(\varphi^{\bm{z}}_{\tau}(\bm{z}_{t}),\bm{h}_{t-1})$
(4)

where the decoder $\varphi^{\text{dec}}_{\tau}$ and the feature extractor
$\varphi^{\bm{z}}_{\tau}$ can be modeled by two-layer neural networks. In
order to improve expressive capability, we also add feature extractor
$\varphi^{\bm{z}}_{\tau}$ to latent variables, which can learn multi-scale and
hierarchical features embedded in the time series. It can significantly
improve the prediction performance on complex and highly dynamic time series.

Inside the LSTM cell, the hidden state is updated as follows:

$\bm{h}_{t}=f_{\theta}(\varphi^{\bm{x}}_{\tau}(\bm{x}_{t}),\varphi^{\bm{z}}_{\tau}(\bm{z}_{t}),\bm{h}_{t-1})$
(5)

where transition function $f$ is essentially the operation described in (2).
We implement function $f$ by multi-layer neural networks with weights
$\theta$. Different from previous LSTM work [13][17][18], the state transition
incorporates the feature of input and feature of latent variables and previous
state, rather than the original input and latent variables.

###  4.3 Objective Function

The target function is the prediction error regularized by the L2 norm of all
trainable weights. We can describe this target as below:

$\sum_{t=1}^{T}\|\hat{\bm{y}}_{t}-\bm{y}\|+\beta\|\bm{W}\|_{2}$ (6)

where $\hat{\bm{y}}_{t}$ and $\bm{y}$ are the predicted values and ground
truth at time $t$, $\bm{W}$ denotes the weights of all neural netwokrs in the
proposed model, and $\beta$ is the trade-off coefficient for the regularizor.
With all functions implemented by neural networks, the accumulative ELBO (6)
is minimized by ADAM [20] algorithm. $\beta$ is chosen to be specific to the
experiment.

##  5 Simulation Results

This section presents a comprehensive experimental study that evaluates the
performance of augmented LSTM method (A-LSTM) in terms of predictive ability.
The performance is compared with cutting-edge time-series prediction models.
For single-step ahead prediction task, we compare the proposed method with the
deep belief network trained with restricted Boltzmann Machine (RBM) [9],
stacked autoencoder model (SAE) [11], and autoencoder stacked on LSTM (Auto-
LSTM) [17]. For multi-step ahead prediction task, the performance comparison
is conducted with multi-output support vector regression (MSVR) [25] and co-
evolutionary multi-task learning (CMTL) [26]. We conduct 5-step and 10-step
ahead prediction experiments, and the scalability is also discussed. All these
methods above have been introduced in Section 3.

The experiments are conducted using three simulated and two real-world time
series datasets. The simulated time series are chaotic nonlinear systems such
as Mackey and Glass [21], Lorenz [22], and Rossler [23]. One real world
dataset is individual household electric power consumption dataset [24]. All
the respective time series data sets are scaled in the range $[0,1]$ in order
to be used for sigmoid units in the feedforward neural network. The
performance is measured by the root-mean-squared error (RMSE):

$\text{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}}$

where $y_{i},\hat{y}_{i}$ are observed and predicted data.

###  5.1 Chaotic Processes Prediction

Across four simulations in this section, the length of training sequence is
1000, and the length of testing sequence is 500. The optimization is performed
by ADAM with learning rate $0.0005$.

####  5.1.1 Mackey and Glass Time Series

We first experiment on Mackey-Glass time series [21], which is derived from
differential equation as below:

$\frac{dx}{dt}=\beta x(t)+\frac{\alpha x(t-\delta)}{1+x(t-\delta)^{10}}$ (7)

The parameters are chosen to be $\alpha=0.2,\beta=-0.1$ and $\delta=17$. The
dataset is constructed by second-order Runge-Kutta method with the step size
of 0.1, with plot as below.

Figure 1: Plot of Mackey and Glass Time Series.

For one-step ahead prediction, the input is a sequence of 5 samples with
stride of 6. For 5-step ahead prediction, the input is a sequence of 15
samples with stride of 6. For 10-step ahead prediction, the input is a
sequence of 10 samples with stride of 6. The RMSE and its variance of the
proposed method (ALSTM), RBM, SAE and Auto-LSTM are shown in tables in the
following.

| ALSTM | RBM | SAE | Auto-LSTM
---|---|---|---|---
RMSE | 0.004778 | 0.010218 | 0.010233 | 0.008231
Table 1: Performance Comparison for One-step Ahead Prediction

Figure 2: RMSE for 5-step Ahead Prediction

Figure 3: RMSE for 10-step Ahead Prediction

| ALSTM | MSVR | CMTL
---|---|---|---
1 Step | 0.010839 | 0.022342 | 0.021562
2 Step | 0.011292 | 0.026897 | 0.019023
3 Step | 0.012740 | 0.035123 | 0.015021
4 Step | 0.013451 | 0.045213 | 0.010228
5 Step | 0.013305 | 0.054234 | 0.009342
Table 2: Variance for 5-step Ahead Prediction

The variance of ALSTM for 10-step ahead prediction is 0.011232, which is
comparable with CMTL and MSVR The complete variance table is omitted.

####  5.1.2 Lorenz Time Series

The Lorenz attractor [22] is a nonlinear 3-D system that provides a simplified
model of convection in the atmosphere. The Lorenz equations are given by:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle\sigma(y-x)$
$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle\rho x-y-xz$
$\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle xy-\beta z$

In this experiment, the parameters are set to be $\sigma=10,\rho=28$ and
$\beta=8/3$. The 3-D attractor and $x,y$ components are plotted as below.

Figure 4: Plots of Lorenz Attractor

The performance of one-step ahead prediction is shown in the following table.

| ALSTM | RBM | SAE | Auto-LSTM
---|---|---|---|---
RMSE | 0.006752 | 0.043278 | 0.020632 | 0.018992
Table 3: Performance Comparison for One-step Ahead Prediction

Figure 5: One-step prediction performance on test data. Left: plots of
confidence interval and mean of generating distribution. Right: plot of error,
which is the difference between the mean of generating distribution and ground
truth.

For 5-step ahead, we use 6 samples with stride of 6. For 10-step ahead, we use
10 samples with stride of 4. The results for 5-step and 10-step ahead
prediction are shown below.

Figure 6: RMSE for 5-step and 10-step Ahead Prediction on Lorenz | ALSTM | MSVR | CMTL
---|---|---|---
1 Step | 0.011467 | 0.015583 | 0.014023
2 Step | 0.010589 | 0.018234 | 0.012882
3 Step | 0.012034 | 0.018199 | 0.010568
4 Step | 0.012589 | 0.019202 | 0.010228
5 Step | 0.014892 | 0.022021 | 0.011772
Table 4: Variance for 5-step Ahead Prediction

####  5.1.3 Rossler Time Series

The Rossler attractor [23] is generated by a system of three differential
equations as below:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle-z-y$
$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle x+ay$
$\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle b+z(x-c)$

where coefficients are set to be $a=0.2,b=0.2$ and $c=4.6$. The figure below
plots the 3-D Rossler Attractor and its $x,y$ components.

Figure 7: Plots of Rossler Attractor. Right: blue line represents $x$ and red
line shows $y$.

For one-step ahead prediction, we use 5 samples with stride of 6. The plots of
prediction results are shown in the following figure.

Figure 8: One-step prediction performance on test data. Left: plots of
confidence interval, mean of generating distribution and ground truth. Right:
prediction error.

The performance comparison for one-step ahead prediction is shown as below:

| ALSTM | RBM | SAE | Auto-LSTM
---|---|---|---|---
RMSE | 0.006824 | 0.013278 | 0.010632 | 0.012992
Variance | 0.010558 | 0.009832 | 0.011342 | 0.012421

The performance comparison for 5-step and 10-step ahead predictions are shown
below.

Figure 9: RMSE for 5-step and 10-step Ahead Prediction on Rossler

The variance of ALSTM is 0.103214, and the complete table for variance is
omitted.

###  5.2 Power Consumption Data

This dataset contains measurements of electric power consumption in one
household with a one-minute sampling rate over a period of almost 4 years
[24], which is equivalent to 7000 data points. The first 1000 data points are
used for training and the remaining for testing.

The experiment shows that ALSTM model can also work well on real-world data.
The following plots show the estimated value and its error, compared with
ground truth.

Figure 10: One-step prediction performance on test data. Left: plots of predicted data (green) and ground truth (blue). Right: plot of prediction error. | ALSTM | RBM | SAE | Auto-LSTM
---|---|---|---|---
RMSE | 0.011562 | 0.035586 | 0.030211 | 0.022520
Table 5: One-step Ahead Prediction Preformance Comparison

Figure 11: RMSE for 5-step and 10-step Ahead Prediction on Power Data

Comparing the performance of experiments with different number of output, we
can see that performance of ALSTM model doesn’t degrade in multi-output cases,
which proves the scalability in practical applications.

##  6 Conclusion

This paper shows the advantage of ALSTM model on time series prediction. With
comprehensive experiments, we prove that the ALSTM model can outperform
cutting-edge algorithm on time series prediction and improve the state-of-art
performance on both chaotic benchmarks and real-world datasets. Moreover, the
predictions from ALSTM model have comparable variance with previous methods,
resulting in small confidence intervals. Comparing performance of different
number of prediction outputs, we can see the performance of 10-step prediction
doesn’t degrade too much, showing that the proposed method has strong ability
of scaling to multi-output prediction. In the future, we are going to apply
ALSTM to predict the financial data which has more variance and uncertainty.

## References

  * [1] G. Box and G. Jenkins, Time Series Analysis: Forecasting and Control. San Francisco, CA: Holden-Day, 1970.
  * [2] F.-A. Gers, J. Schmidhuber, F. Cummins. Learning to Forget: Continual Prediction with LSTM. Technical Report IDSIA-01-99, 1999.
  * [3] S. Hochreiter and J. Schmidhuber. Long short-term memory. Neural Computation, 9(8):1735-1780, 1997.
  * [4] A. Graves. Generating sequences with recurrent neural networks. arXiv preprint arXiv:1308.0850, 2013.
  * [5] F.A. Gers, D. Eck, and J. Schmidhuber, Applying LSTM to time series predictable through time-window approaches. Neural Nets WIRN Vietri-01. Springer, London, 2002. 193-200.
  * [6] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning . Nature, 521(7553):436-444, 2015.
  * [7] W. Huang, G. Song, H. Hong and K. Xie. Deep architecture for traffic flow prediction: deep belief networks with multi-task learning. IEEE Transactions on Intelligent Transportation Systems 15(5) 2191-2201, 2014.
  * [8] Y. Lv, Y. Duan, W. Kang, Z. Li, and F.-Y. Wang. Traffic flow prediction with big data: a deep learning approach. IEEE Transactions on Intelligent Transportation Systems,16(2):865-873, 2015.
  * [9] T. Kuremoto, S. Kimura, K. Kobayashi, and M. Obayashi. Time series forecasting using a deep belief network with restricted Boltzmann machines. Neurocomputing 137 (2014): 47-56.
  * [10] J.-T. Turner. Time series analysis using deep feed forward neural networks. Ph.D. thesis, University of Maryland, Baltimore County, 2014.
  * [11] P. Romeu, F. Zamora-Martı́nez, P. Botella-Rocamora, J. Pardo. Time-series forecasting of indoor temperature using pre-trained deep neural networks. In: Artificial Neural Networks and Machine Learning–ICANN 2013, pp. 451-458, Springer ,2013.
  * [12] P. Kanestrom. Traffic flow forecasting with deep learning. MS thesis. NTNU, 2017.
  * [13] X. Ma, et al. Long short-term memory neural network for traffic speed prediction using remote microwave sensor data. Transportation Research Part C: Emerging Technologies 54 : 187-197, 2015.
  * [14] G. Hinton, and R. Salakhutdinov. Reducing the dimensionality of data with neural networks. Science,313.5786:504-507, 2006.
  * [15] Y. Cinar, H. Mirisaee, P. Goswami, etc. Time Series Forecasting using RNNs: an Extended Attention Mechanism to Model Periods and Handle Missing Values. arXiv preprint arXiv:1703.10089.
  * [16] H.-M. Nguyen, S. Woo, J. Im, T. Jun, and D. Kim. A Workload Prediction Approach Using Models Stacking Based on Recurrent Neural Network and Autoencoder. High Performance Computing and Communications; IEEE 14th International Conference on Smart City, 2016.
  * [17] A. Gensler, J. Henze, B. Sick, and N. Raabe. Deep Learning for solar power forecasting—An approach using AutoEncoder and LSTM Neural Networks. In Systems, Man, and Cybernetics (SMC), 2016 IEEE International Conference on (pp. 002858-002865), Oct. 2016.
  * [18] F.-M. Bianchi, et al. An overview and comparative analysis of Recurrent Neural Networks for Short Term Load Forecasting. arXiv preprint arXiv:1705.04378, 2017.
  * [19] Bengio, Y., 2007. Learning deep architectures for AI. Technical Report 1312. Dept. IRO, Universite de Montreal.
  * [20] D. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  * [21] M. Mackey and L. Glass. Oscillation and chaos in physiological control systems, Science 197 (4300), pages 287-289, 1977.
  * [22] E. Lorenz. Deterministic non-periodic flows, J. Atmos. Sci., vol. 20, pages 267-285, 1963.
  * [23] O. E. Rossler. An equation for continuous chaos. Physics Letters A, 57(5), 397-398, 1976.
  * [24] UCI Machine Learning Repository. Available at http://archive.ics.uci.edu/ml/datasets
  * [25] Y. Bao, T. Xiong, and Z. Hu. Multi-step-ahead time series prediction using multiple-output support vector regression. Neurocomputing 129 (2014): 482-493.
  * [26] R. Chandra, Y.-S. Ong, and C.-K. Goh. Co-evolutionary multi-task learning for dynamic time series prediction.” arXiv preprint arXiv:1703.01887, 2017.