# Classification of events from $\alpha$-induced reactions in the MUSIC
detector via statistical and ML methods

K. Raghavan  M. L. Avila  P. Balaprakash  H. Jayatissa  D. Santiago-Gonzalez
<EMAIL_ADDRESS>Mathematics and Computer Science Division, Argonne National
Laboratory, Lemont, IL 60439, USA  Physics Division, Argonne National
Laboratory, Lemont, IL 60439, USA

###### Abstract

The Multi-Sampling Ionization Chamber (MUSIC) detector is typically used to
measure nuclear reaction cross sections relevant for nuclear astrophysics,
fusion studies, and other applications. From the MUSIC data produced in one
experiment scientists carefully extract an order of $10^{3}$ events of
interest from about $10^{9}$ total events, where each event can be represented
by an 18-dimensional vector. However, the standard data classification process
is based on expert-driven, manually intensive data analysis techniques that
require several months to identify patterns and classify the relevant events
from the collected data. To address this issue, we present a method for the
classification of events originating from specific $\alpha$-induced reactions
by combining statistical and machine learning methods that require
significantly less input from the domain scientist, relative to the standard
technique. We applied the new method to two experimental data sets and
compared our results with those obtained using traditional methods. With few
exceptions, the number of events classified by our method agrees within $\pm
20\%$ with the results obtained using traditional methods. With the present
method, which is the first of its kind for the MUSIC data, we have established
the foundation for the automated extraction of physical events of interest
from experiments using the MUSIC detector.

###### keywords:

AI/ML methods, $\alpha$-induced reactions , active target, experimental
nuclear physics

††journal: Nuclear Instruments and Methods in Physics Research Section A:
Accelerators, Spectrometers, Detectors and Associated Equipment

##  1 Introduction

The quest for answers to some of the most fundamental questions in nuclear
physics has led to the design and construction of increasingly complex
experimental equipment and facilities. Through careful analysis of the
experimental data generated by detector systems, nuclear physicists make
discoveries and measure properties of atomic nuclei for basic research and
diverse applications. The scale of the experiments, detector systems, and
collected data demands advances in current analysis techniques in order to
efficiently extract the relevant information from such large and complex data
sets. At present, expert-driven, manually intensive data analysis techniques
require several months to identify patterns and classify the relevant events
from the collected data. Artificial intelligence (AI) techniques, in
particular machine learning (ML) methods, have the potential to automate,
accelerate, and improve the data analysis techniques, thereby reducing the
effort to obtain the physical data, enhancing the scientific throughput, and
potentially enabling further discoveries.

From the AI/ML perspective, the analysis of data produced by the detector
systems poses a number of unique challenges that need the development of
domain-specific AI/ML methods. For example, in the case of the Multi-Sampling
Ionization Chamber (MUSIC) detector data, the number of background events is
typically several orders of magnitude more than the events of interest,
rendering out-of-the-box anomaly detection methods ineffective. A few recent
examples of domain-specific AI/ML methods can be found in the works of Kuchera
et al. [1] (experiment) and Raghavan et al. [2] (theory). For a recent review,
see [3].

In this paper we present a method for the classification of data obtained from
MUSIC, an active-target detector specializing in cross-section measurements of
various nuclear reactions such as fusion [4] and $\alpha$-induced reactions
[5].

Our proposed approach comprises four phases: (1) filter easy-to-identify
background events; (2) detect and remove background events that are similar to
the events of interest; (3) get input from the domain scientists and design a
classifier; and (4) detect the strip location at which the event took place.
We leverage statistical and machine learning methods in phases 1–3. In
addition to experimental data, we use simulated data in phases 1–3 to assist
in the event classification.

We applied our method to two data sets obtained under different experimental
conditions. In the first data set, $\alpha$-induced reactions on Fluorine-17
(17F) nuclei were observed and classified with standard analysis techniques.
In the second data set, $\alpha$-induced reactions on Oxygen-17 (17O) nuclei
were studied. The 17O data was previously classified with standard techniques,
and those results are published in [5].

We show that our new method classifies the MUSIC data—that is, it separates
relevant events for a particular nuclear reaction from other events in a given
MUSIC strip—in a way such that the percent differences with respect to human
expert classification are typically within $\pm 20\%$.

This paper is organized as follows. In Sec. 2 we briefly describe the MUSIC
detector and its data structure. In Sec. 3 we give a brief description of the
experiments that produced the data analyzed by the methods presented in this
paper. In Sec. 4 we present basic information on the simulation of the
detector response. In Sec. 5 we describe our data classification method, which
is based on statistical and ML methods. In Sec. 6 and 7 we present results and
conclusions, respectively.

##  2 The MUSIC detector and its data structure

MUSIC is an active-target gas-filled ionization chamber, where the gas serves
as both target and detection material. The MUSIC detector has been used to
measure excitation functions of heavy-ion fusion reactions [6, 4] and to
measure $(\alpha,p)$ and $(\alpha,n)$ reaction cross sections in inverse
kinematics [5, 7, 8, 9], involving stable and radioactive beams. For the
$(\alpha,p)$ and $(\alpha,n)$ reactions, cross sections down to 1 mb have been
measured. If the $(\alpha,p)$ and $(\alpha,n)$ channels are both energetically
allowed at the energies studied, then the MUSIC detector is able to probe them
simultaneously [7]. The high detection efficiency of MUSIC (close to 100%)
makes the most of heavy-ion beams with relatively low intensity (less than
$10^{6}$ particles per second). In addition, MUSIC can measure averaged cross-
section points at different beam energies with a single incident beam energy.
The reason is that the incident beam particles lose energy through
interactions with the gas molecules as they move deeper into the detector.

The anode of MUSIC has 18 segments, with the middle 16 segments subdivided in
left and right as shown in Fig. 1. With this segmentation, the total number of
electronic channels in the anode is 34. The data from the left and right
segmentation has proven useful when studying heavy-ion fusion reactions, since
the outgoing particles from elastic scattering events deposit enough energy in
both the left and right segments of a single strip (multiplicity 2) and can be
clearly separated from fusion events that have multiplicity 1; see [4]. For
$\alpha$-induced reactions, however, separation of elastic scattering events
from the events of interest based on “multiplicity” is not as effective since
the energy deposited by the elastically scattered helium ions in the gas is
such that it is normally comparable to the electronics noise. Therefore, only
the sum of the left and right segments was used for the ML analysis of the
$\alpha$-induced reactions presented in this paper and not the individual
contributions for left and right. This reduces the number of dimensions of an
event from 34 to 18. For a more detailed description of the detector and the
standard data analysis technique, see [6, 5].

Figure 1:  Schematic of the MUSIC detector anode structure. There are 18
strips composed by 34 electronics channels.

While the data collection of an experiment performed with the MUSIC detector
is carried out in a couple of days, it usually takes a human expert several
months to both find the relevant information from the $\sim 10^{9}$ recorded
events and to classify the $\sim 10^{3}$ events of interest, in order to then
extract the reaction cross sections. In order to separate the background
events from the events of interest, henceforth referred to as golden events,
vectors containing the particle’s energy loss in each detector segment are
studied. The energy loss or energy deposited by different particles in a
segment is described by the Bethe formula [10] that, in the present energy
regime, is proportional to the square of the number of protons, $Z$.
Therefore, under appropriate experimental conditions, particles in MUSIC can
be separated if they have different $Z$ number. To distinguish between
different reactions, the energy loss in several MUSIC segments, after the
reaction occurs, is summed and plotted against the total energy deposited in
the detector—the so-called $\Delta E$–$E$ plots. Depending on the energy and
type of reaction, sometimes several permutations of the plots have to be
performed in order to optimize the separation of different events. Then, two-
dimensional graphical “cuts” are used to select promising events, which are
later classified by studying the characteristic shape of their energy loss
profiles across the segmented anode strips, known as event traces. Event
traces are further discussed in Sec. 4 (examples are shown in Fig. 2).

Table 1 shows basic information about the data structure of an event in the
MUSIC detector. Currently, the data from each of the 18 strips of the MUSIC
detector can be acquired by using analog or digital electronics. In the analog
case, 10-bit Philips analog to digital converters, controlled by a Computer-
Aided Measurement And Control (CAMAC) crate, are used. In the digital case,
14-bit CAEN digitizers, controlled by a Versa Module Eurocard (VME) crate, are
used.

Dimension | 18
---|---
Raw data type | integer
Raw data range | 0 to 4,096, or 0 to 16,384
Calibrated data type | real (floating point)
Calibrated data range | depends on experimental conditions
Table 1: Data structure of a MUSIC event (point). Raw data range depends on
the data acquisition system (10- or 14-bit) used, but it is constant during a
single experiment.

##  3 Experimental data

The experimental data shown in this paper was obtained with the MUSIC detector
from experiments aimed to probe the $\alpha$-induced reactions
17F$(\alpha,p)^{20}$Ne and 17O$(\alpha,n)^{20}$Ne. These experiments were
carried out at the Argonne Tandem Linac Accelerator System (ATLAS) in inverse
kinematics. Both of these reactions are relevant for nuclear astrophysics and
were measured at low energies. The beam energies for the 17F and 17O beams
were about 50 MeV and 35 MeV, respectively. The 17F beam is radioactive and
was produced by using the in-flight technique via the 16O$(d,n)^{17}$F
reaction. The maximum beam intensity of the 17F beam was about 3000 pps. The
MUSIC detector was filled with pure 4He gas at room temperature and a pressure
of 400 Torr for the 17F$+\alpha$ case and 206 Torr for the 17O$+\alpha$ case.
In both cases, the data acquisition system was based on analog electronics.
The physics motivation, final cross section, and more details of the
17F$(\alpha,p)^{20}$Ne reaction will be presented in a future publication. The
cross section of the 17O$(\alpha,n)^{20}$Ne reaction and more details of the
experiment have been published in [5].

##  4 Simulated data

The computer program musicsim [11] was written to simulate the detector
response based on the kinematics of the reactions occurring in the detector,
the geometry of the detector, and the energy loss of the particles as they
travel through the gas (detector medium). The simulation takes as inputs the
kinetic energy of the beam particles as they enter the gas volume, the beam
energy spread, the random noise for all the detector segments, the ion
species—e.g., 12C for a Carbon-12 nucleus—for initial (beam and target), and
final (i.e., after a nuclear reaction) particles. To simulate a nuclear
reaction in a specific detector strip, the program randomly selects a point of
interaction along the trajectory of the beam particle (within the selected
strip) and calculates the energy loss of the beam as it travels from the
entrance window to the interaction point. To calculate the energy loss, one
can use stopping powers from different tables, such as the Ziegler [12] and
ATIMA [13]. After the energy is calculated, the simulation program randomly
selects an angle in the center of mass and calculates the kinematics of the
outgoing particles, then boosts the 4-momenta of the outgoing particles from
the center of mass to the laboratory frame through a Lorentz transformation.
Using the angles and energies of the different particles (in the laboratory
frame), the simulation program uses energy loss tables to calculate the energy
deposited per strip, as it is measured in the MUSIC detector.

Figure 2 shows examples of simulated (histogram on the background) and
experimental (blue lines on the foreground) data for the 17O$+\alpha$ case. In
this data visualization, the total energy deposited in each strip is plotted
as a function of the strip number. For the strips with left/right
segmentation, the y-axis is composed of the sum of the energy deposited in the
left and right segments. Three cases are shown in the figure: golden events
(top), which in this case correspond to events from 17O$(\alpha,n)^{20}$Ne
reactions; not golden events, in this case 17O$(\alpha,\alpha)^{17}$O
reactions; and background (bottom), for example from the unreacted 17O beam,
which are the most abundant background events but experimentally could include
other spurious events and noise. For the histograms (from simulation) the bins
with dark color indicate a higher frequency than those with lighter colors.
For the 17O$(\alpha,n)$ and 17O$(\alpha,\alpha)$ simulations, an isotropic
angular distribution (in the center of mass) has been assumed; however, the
actual physical distribution is unlikely to be isotropic.

---
Figure 2:  Comparison of simulated (histograms) and experimental (blue lines)
data for the 17O$+\alpha$ case. Several individual experimental traces are
shown with blue lines and are displayed over the histograms for (bottom)
unreacted 17O beam and for events occurring in strip 5 classified as “not
golden events” (middle), and “golden events” (top).

Empirically, separation between golden events (GE) and not golden events (not-
GE) becomes increasingly more challenging as we probe strips located deeper
into the detector, where the energy loss differences are smaller. The
separation between GE and not-GE also is expected to diminish when the
outgoing heavy ion differs by only one proton from the beam, as in the case of
17F$(\alpha,p)^{20}$Ne.

##  5 AI/ML methods for MUSIC

The experimental data obtained with the MUSIC detector can be denoted as a
matrix $\mathcal{X}\in\mathbb{R}^{n\times p}$ that is a collection of events
from the detector, where $n$ is the total number of events data from the
detector and $p$ is the dimension of the each event vector containing the
particle’s energy loss in each detector strip. As indicated in Sec. 2, the
nominal value of $p$ is $18$, and the typical value of $n$ is $\sim 10^{9}$.
We divide $\mathcal{X}$ into background ($\mathcal{X}_{\mathrm{back}}),$
golden events (events of interest, $\mathcal{X}_{GE}$), and other events
($\mathcal{X}_{OE}$ with
$\mathcal{X}=\mathcal{X}_{\mathrm{back}}\cup\mathcal{X}_{GE}\cup\mathcal{X}_{OE}$).
The exact proportion of the background, golden, and other events is unknown.
The goal of the analysis is to extract $\mathcal{X}_{GE}$ from $\mathcal{X}$
under the assumption that the total number of background events is
significantly larger than for golden events. For example, in the case of
$\alpha$-induced reactions on 17F, the golden events $\mathcal{X}_{GE}$ are
$(\alpha,p)$ events. Similarly, $(\alpha,n)$ events are the golden events for
$\alpha$-induced reactions on 17O.

The simulation data is denoted as $\mathcal{S}$, which can be divided into
three groups: $\mathcal{S}_{\mathrm{back}}$, typical background events (mostly
from unreacted beam); $\mathcal{S}_{GE}$, golden events (events of interest);
and $\mathcal{S}_{OE}$, the other events. These samples are obtained with the
simulation code described in Sec. 4.

The methodology to extract $\mathcal{X}_{GE}$ comprises four phases. See Fig.
3 for a high-level overview.

Figure 3: Schematic of our four-phase approach to golden event detection in
MUSIC detector.

Phase 1:

Remove easy-to-identify background events from $\mathcal{X}$ using statistical
techniques.

Phase 2:

Use a neural network transformation to remove background events that are
similar to golden events.

Phase 3:

Design a classifier using input from domain scientists to extract three
clusters that indicate events that are surely golden (“Surely GE”); are surely
other events (“Surely OE”); and may or may not be golden events (“Maybe GE”).

Phase 4:

Using the rate of change of energies at each strip, determine the strips
corresponding to the location (detector segment, signified by the strip
number) where the event took place.

Below, each phase is discussed in detail.

###  5.1 Phase 1: Removing background events

The total number of events in $\mathcal{X}$ is around $10^{9}$. The majority
of these (roughly $10^{5}$) are background events with traces similar to those
shown in the bottom panel of Fig. 2. The easy-to-identify background events
are removed by using a distance – a measure-based detection.

The Mahalanobis distance (MD) is a measure of the distance between a point and
a distribution [14]. Typically, one sample is denoted as reference
($\mathcal{R}$ with mean $\bar{r}$ and covariance matrix $\Sigma$) and one
sample as test (in this case $\mathcal{X}$). MD then quantifies the difference
between reference sample $\mathcal{R}$ and a point $x\in\mathcal{X}$ such that
the MD, $d(x,\mathcal{R})$, is given as

$\displaystyle d(x,\mathcal{R})=\|\Sigma^{-1/2}(x-\bar{r})\|_{2}^{2},$ (1)

where $\Sigma^{-1/2}$ denotes the Cholesky decomposition [15] of the
covariance matrix and $\|.\|_{2}^{2}$ is the squared 2-norm. In practice,
$\Sigma$ is replaced by its sample estimate $\hat{\Sigma}$. Typically, for any
normally distributed reference and a test point $x$, $d(x,\mathcal{R})<1$
would imply that $x$ belongs to the same distribution as the reference.
Therefore, with the goal of removing background events, one may obtain a
reference sample that is representative of background events and designate any
event $x\in\mathcal{X}$ as background, if $d(x,\mathcal{R})$ is less than a
threshold $d_{B}.$

Formally, the background events can be removed by testing the null hypothesis

$H_{0}:x\sim p(\mathcal{R}),$

where $\mathcal{R}$ is the reference sample and $p(\mathcal{R})$ is a
distribution over it. In order to test $H_{0}$ based on MD, a distance
$d(x,\mathcal{R})$ is calculated for each $x\in\mathcal{X}.$ Therefore, for
rows $x_{i}$ in $\mathcal{X}$, where $i=1,\cdots,n$, the distances are denoted
as $\mathcal{D}=\\{d_{1},d_{2},\cdots d_{N}\\}$. With sample $\mathcal{D},$ a
sampling distribution $p(\mathcal{D})$ and a distance threshold $(d_{B})$,
$\forall d\in\mathcal{D},$ the null and alternate hypothesis can be described
as

$\displaystyle H_{0}:d\leq d_{B}\quad\quad H_{1}:d>d_{B}.$ (2)

Satisfying the null hypothesis indicates that $x\sim p(\mathcal{R})$ (i.e.,
$x$ is a background event), and satisfying the alternate hypothesis indicates
the $x$ may be a golden event. To perform this hypothesis testing, one must
obtain an accurate $\mathcal{R}$ reference sample and an accurate $d_{B}$
distance threshold.

Reference Sample: With the assumption that the total number of background
events in $\mathcal{X}$ far exceeds the total number of golden events, a
uniform sample from $\mathcal{X}$ denoted as $\tilde{X}_{b}$ provides a
reasonable representation of the background events. Furthermore, a simulation
sample $\mathcal{S}_{\mathrm{back}}$ corresponding to the background events is
also provided, to be used along with $\tilde{X}_{b}$ to constitute the
reference sample—$\mathcal{R}=\mathcal{S}_{\mathrm{back}}\cup\tilde{X}_{b}$.

Distance Threshold: The distance threshold is obtained by designating
$d_{B}=P_{r}$, where $P_{r}$ denotes the $r$th percentile of $p(\mathcal{D}).$
A common choice of $r$ is $99$ such that $P_{r}$ denotes the $99$th
percentile. We note that a low value of $r$ yields a large type one error
(i.e., more background events are misclassified as potential golden events).
All the $x\in\mathcal{X}$ that satisfies the null hypothesis are designated as
background events; the rest of the events are collectively denoted as
$\mathcal{X}_{PE}$.

###  5.2 Phase 2- Outlier detection

As previously mentioned, a large number of events in $\mathcal{X}$ are
background events originating from unreacted beam particles. Since the
reference sample is obtained by uniformly sampling $\mathcal{X},$ it is a good
representation of an unreacted beam. MD-based detection is limited to removing
such background events. Thus, $\mathcal{X}_{PE}$ comprises golden events,
other events and many background events that are similar in their trace to
golden events or other events. Therefore, to precisely extract golden events,
one must remove other events and leftover background events from
$\mathcal{X}_{PE}$. We will first utilize DBSCAN (density-based spatial
clustering of applications with noise) to remove other events and develop a
neural-network-based transformation to reduce similarities between golden
events and background events.

####  5.2.1 DBSCAN – Density-Based Spatial Clustering of Applications with
Noise

DBSCAN [16] is a nonparametric algorithm where, given a set of points, DBSCAN
finds points that are closely packed (points with many nearby neighbors—the
high-density region) and points whose nearest neighbors are too far away (the
low-density region). Under the assumption that all the “OE” will cluster
together in the high-density regions, DBSCAN is used to extract points that
are closely packed, and these points are designated as $\mathcal{X}_{OE}.$
Similarly, the events that lie in the low-density region are designated as
$\mathcal{X}_{PGE}$, which is referred as the set of all potential golden
events such that $\mathcal{X}_{PGE}$ is a collection of golden events and
background events that are similar in their trace to golden events. In
particular, we would like to design a classifier that can take as input
$\mathcal{X}_{PGE}$ and label them as other events and golden events. To
design this classifier, however, we need to obtain labeled data. The next two
steps are therefore employed to extract that core subset of data that can be
labeled by the domain scientists.

####  5.2.2 Neural Network-Based Transformation

It has been empirically observed that in cases such as the one for
17F$(\alpha,p)$ reactions, several background events are similar to golden
events. In such scenarios one can perform a transformation that reduces the
similarity between golden and hard-to-identify background events. In order to
construct this transformation, an arbitrary neural network (NN) can be
employed. NNs comprise layers where each layer is built by using a linear
transformation employing weights and a nonlinear function. Typically, the
weights of the NN can be trained by minimizing a loss function. The weights
are collectively denoted as $\theta$ such that the NN can be denoted as
$f_{\theta}.$ In this work, the network is trained to expand the differences
between golden and background events with the goal of minimizing a cost
$J(\theta)$ such that

$\theta^{*}=\mathrm{argmin}_{\theta}J(\theta).$

The cost $J(\theta)$ is designed to quantify the differences between golden
and background events. To this end, three sample sets are defined: the
background
$\mathcal{T}_{\mathrm{back}}=\mathcal{S}_{\mathrm{back}}\cup\mathcal{R}$;
golden events $\mathcal{T}_{GE}=\mathcal{S}_{GE}$; and other events
$\mathcal{T}_{OE}=\mathcal{S}_{OE}\cup\mathcal{R}_{OE},$ where
$\mathcal{R}_{OE}$ is a sample from $\mathcal{X}_{OE}$ obtained at the end of
Subsection 5.2.1. Equipped with these samples, we define $J(\theta)$ as
follows.

$\displaystyle
J(\theta)=\mathbb{E}_{x\in\mathcal{T}_{\mathrm{back}},y\in\mathcal{T}_{GE},z\in\mathcal{T}_{OE}}\bigg{[}\|z-f_{\theta}(x)\|_{2}^{2}-\gamma\big{[}\|x-f_{\theta}(y)\|_{2}^{2}+\|z-f_{\theta}(y)\|_{2}^{2}\big{]}\bigg{]}$
(3)

where $\mathbb{E}$ is the expected value operator. In this cost function, the
distance between a golden event and reference/other events is increased by
maximizing the term
$[\|x-f_{\theta}(y)\|_{2}^{2}+\|z-f_{\theta}(y)\|_{2}^{2}]$, and the distance
between the other-events and background sample is reduced by minimizing the
term $\|z-f_{\theta}(x)\|_{2}^{2}$. The $\theta^{*}$ that satisfies the
objective Eq. (3) increases the distance between golden events and the rest
(other events and background). This transformation can be estimated by using
standard learning techniques such as stochastic gradient descent [17]. Once
this transformation is learned, $f({\theta^{*}})$ is used on all
$\mathcal{X}_{PGE}$ such that
$\mathcal{X}^{T}_{PGE}=f_{\theta^{*}}(\mathcal{X}_{PGE}).$ From here on, the
application of the NN transformation is denoted by a superscript T. Once we
transform all $\mathcal{X}_{PGE}$, a nearest-neighbor-based scheme is utilized
to remove most of the leftover background events and other events.

####  5.2.3 Nearest-Neighbor-Based Outlier Detection

The events from $\mathcal{X}^{T}_{PGE}$ that look similar to background and
other events are removed by using a nearest-neighbor-driven search. In order
to efficiently perform this search, a tree-based approach [18] is employed,
and a k-d tree (hereafter denoted as $Kdt$) is defined on $\mathcal{R}^{T},$
where $\mathcal{R}^{T}$ is the reference sample transformed using the NN
(refer to Section 5.2.2). With this tree, a nearest-neighbor search algorithm
[19] is executed to obtain k-neighbors from $\mathcal{R}^{T}$ that are closest
to a given $x\in\mathcal{X}^{T}_{PGE}.$ Based on these k neighbors, one
determines whether an event is golden or not according to the average distance
between $x$ and different neighbors.

Formally, for every $x\in\mathcal{X}^{T}_{PGE},$ $k$ neighbors are first
determined by using $Kdt.$ Next, the average distance between $x$ and the $k$
neighbors ($d_{\mathrm{avg}}$) is evaluated, where
$d_{\mathrm{avg}}=(\sum_{i=1}^{k}d_{i})/k)$ with $d_{i}$ being the distance
between each neighbor and $x.$ Based on a threshold $d_{\mathrm{thre}},$ the
following hypothesis is tested for each $x$.

$\displaystyle H_{0}:d\leq d_{\mathrm{thre}}\implies x\sim
p(\mathcal{R}^{T})\quad H_{1}:d>d_{\mathrm{thre}}\implies x\sim
p(\mathcal{T}_{GE})$ (4)

For any $x,$ if the null hypothesis is accepted, a background event is
detected; if the alternate hypothesis is accepted, a golden event is detected.
The strategy utilized in phase 1 can be used to determine $d_{\mathrm{thre}}.$
Based on the results of this hypothesis testing, outliers are extracted and
designated as $O.$ The set of outliers $(O)$ are then clustered into multiple
groups by using the Gaussian mixture model and given to the domain scientists
for labeling. The clusters obtained at the end of last phase are labeled by
the domain scientists as either an other events cluster (a collection of other
events clusters is denoted as $\tilde{\mathcal{X}}_{OE}$) or a golden events
cluster (a collection of golden events cluster is denoted as
$\tilde{\mathcal{X}}_{GE}$). These two class samples can now be employed to
classify $\mathcal{X}_{PGE}$ and obtain the final clusters.

###  5.3 Phase 3: Classification

With these two class samples, $\tilde{\mathcal{X}}_{GE}$ and
$\tilde{\mathcal{X}}_{OE}$, a two-step classification method is designed to
determine the three event groups “Surely GE,” “Maybe GE,” and “Surely OE.”
Almost any classification method [19] could be used to determine the three
event groups since the training data sample corresponding to GE and Not GE is
now available. In our preliminary analysis, the results recorded in Table 4 in
Sec. 6 show that the most GE’s were detected by One-class SVM [20] albeit with
many false positives. Since most of the background events (the source of false
positives) have been removed and a good training sample is now available, one-
class SVM is an ideal choice as a classifier in the two-step classification
procedure.

In the first step, the “Surely GE” class is extracted by using an SVM
classifier [19] that is designed to separate “Surely GE” $(\mathcal{L}_{GE})$
from $\mathcal{X}_{PGE}$. Specifically, the GE class for the classifier is
defined as $\tilde{\mathcal{X}}_{GE}\cup\mathcal{S}_{GE}$ and the OE class as
$\tilde{\mathcal{X}}_{OE}\cup\mathcal{S}_{OE}\cup\mathcal{S}_{\mathrm{back}}.$

Second, the “Maybe GE” class is extracted from
$\mathcal{X}_{PGE}|\mathcal{L}_{GE}$, that is, $\mathcal{X}_{PGE}$ with
$\mathcal{L}_{GE}$ removed. The training data for the GE class is redefined as
$\mathcal{S}_{GE}\cup\mathcal{L}_{GE},$ and the OE class is provided as
$\tilde{\mathcal{X}}_{OE}\cup\mathcal{S}_{OE}\cup\mathcal{S}_{\mathrm{back}}.$
Using this new training data, a classifier is used to extract all the
remaining golden events from $\mathcal{X}_{PGE}$, where $\mathcal{L}_{GE}$ are
removed. The result of this classifier is designated as “Maybe GE,” and all
the remaining events are designated as the “Surely OE” class.

###  5.4 Phase 4: Strip Detection

The final step is to determine at which detector segment a golden event
occurred. Experimentally, when a nuclear reaction occurs within one MUSIC
strip, the measured energy deposited in that strip, and the following strips
change in a way that depends on the atomic number of the resulting ions and on
their kinematics. In the $(\alpha,p)$ and $(\alpha,n)$ cases, an event is said
to happen at strip $i$ if there occurs a significantly large increment in the
energy deposited in the strip. This characteristic “jump” can be detected by
using the gradient of the energy level with respect to the strip number (that
is, the gradient is maximum at the strip corresponding to the characteristic
jump). Therefore, for any data point $x,$ the gradient at each strip for
$x=\\{x_{i},i=1\cdots p\\}$ is given as $G^{x}_{i}=\frac{\partial x}{\partial
i}\approx|x_{i}-x_{i-1}|.$ The strip number $s_{x}$ for $x$ then is given as
$s_{x}=\mathrm{argmax}\;G_{x}$ such that $G_{x}=\\{G^{x}_{i},i=1,\cdots p\\}.$
Equipped with this algorithm, we determine the strip numbers for each event
categorized at the end of the last phase, in other words, an event in
the“Surely GE” class, “Maybe GE” class, and “Surely OE” class.

##  6 Results and discussion

The new method described above was developed for the strip-wise classification
of $(\alpha,p)$ and $(\alpha,n)$ reactions taking place in the MUSIC detector.
We have applied the method to two experimental data sets: the 17F$+\alpha$
data set, where 17F$(\alpha,p)^{20}$Ne reactions are the “golden events,” and
the 17O$+\alpha$ data set for the classification of 17O$(\alpha,n)^{20}$Ne
events. In the former case, no published experimental data exists;
nevertheless, a preliminary analysis using traditional classification methods
to tag the 17F$+\alpha$ events has been performed so that the performance of
the new method can be assessed. For the latter case, the results are compared
with the experimental data published in [5]. Based on the physics of the
detector response (see Sec. 2), typically the principal differences in the
MUSIC traces are enhanced when more protons are removed from or added to the
beam particles ($\Delta Z$). Therefore, it is expected that for the
17F$(\alpha,p)^{20}$Ne case ($\Delta Z=1$) the classification of GE would be
more challenging than for the 17O$(\alpha,n)^{20}$Ne case ($\Delta Z=2$). In
this section we discuss the results obtained by applying the new method to
these two MUSIC data sets.

###  6.1 Classification of 17F$(\alpha,p)^{20}$Ne events

In the top panel of Fig. 4, we show the number of events classified as “Surely
GE,” in this case 17F$(\alpha,p)$ events, plotted as a function of the MUSIC
strip number. Such event classification is a key step in the extraction of the
17F$(\alpha,p)$ reaction cross section, which is the quantity of interest for
the domain scientists. In Fig. 4, the black points are obtained by using the
standard technique (traditional method), and the red open circles are obtained
with the new (present) method. In the top panel the error bars of black points
include statistical and systematic uncertainties while the error bars for the
red points relate to the aleatoric and epistemic uncertainties of the new
method (see Sec. 6.3). In the bottom panel we show the percent change in the
number of events classified as “Surely GE” using the new method with respect
to the number of 17F$(\alpha,p)$ events classified using the traditional
method as a function of the strip number. The shaded area encompasses the
region of $\pm 20\%$. Error bars in the bottom graphs are the experimental
ones. The up arrows in strip 1 indicate that the corresponding % change values
are off the scale.

---
Figure 4:  Comparison of results for the classification of 17F$(\alpha,p)$
events with the traditional technique (black points) or with the present
method (red open circles). See text.

The largest disagreement is found in the classification of GE for strip 1.
Such a dramatic increment in the number of counts, going from strip 2 to strip
1, is not expected from the experimental point of view, and it is likely not
physical. This suggests that, for the new method, refinements are needed to
distinguish strip-1 GE from not-GE or from background signals. Intuitively,
this is not too surprising since there is not a robust baseline—in this case
the baseline is just formed by the strip 0 data–upon which to signature “jump”
in the MUSIC traces would be apparent. Empirically, we also find challenging
the classification of GE in the strips near the edges of the detector (in this
case, strips 1 and 2, and beyond strip 10) when using traditional methods.
Therefore, larger systematic uncertainties are expected near the edges of the
detector.

17F$(\alpha,p)$ | GE pres. meth. | Not-GE pres. meth.
---|---|---
GE trad. meth. | 1825(32) | 186(32)
Not-GE trad. meth. | 7968(1311) | 569882(1311)
Table 2: Confusion matrix for the 17F$(\alpha,p)$ case showing the mean number
of events with standard deviation in parentheses.

Table 2 shows the confusion matrix for two types of events classified by the
new method and traditional methods. Again, in this case the golden events (GE)
correspond to 17F$(\alpha,p)$ reactions, and not-GE are events that are not
considered background but rather other types of reactions such as
$(\alpha,\alpha)$ events (background is mostly unreacted beam events).
Physically, the not-GE events likely contain some $\alpha$-scattering events
but might also contain spurious events originating from electronics noise in
some of the detector channels. In the table we show the mean number of events
with standard deviation in parentheses. In order to evaluate these quantities,
the four phases of the new method were executed for 50 repetitions. The values
are rounded off to the nearest integer.

As mentioned above, the large number of false positives, namely, the 7968
events in Table 2, originate from the disagreement in classification of events
in strip 1. When strip 1 events are removed from the analysis, the number of
false positives is dramatically reduced, from 7968 to 279 events, while the
number of true positives remains high, changing from 1825 to 1736 events. The
number of false negatives, namely, the 186 events in Table 2, is less than 10%
of the total number of GE classified using traditional methods. We compare the
accuracy of the present method with other outlier detection methods in Sec.
6.4.

###  6.2 Classification of 17O$(\alpha,n)^{20}$Ne events

In Fig. 5 we compare the results for the classification of 17O$(\alpha,n)$
events obtained from the traditional method and the new method. As previously
mentioned, the experimental data was published in [5]. In the top panel, we
show the number of events classified as 17O$(\alpha,n)$ by the standard
technique (black circles) or by the present method (red open circles) plotted
as a function of the MUSIC strip number. As in the preceding section, the
error bars of the black points are only statistical, while the error bars for
the red points relate to the aleatoric and epistemic uncertainties of the
present model (see Sec. 6.3). In the bottom panel, we show the percent change
of number of events classified as “Surely GE,” using the new method, with
respect to the number of 17O$(\alpha,n)$ events classified using traditional
technique as a function of the strip number. The shaded area encompasses the
region of $\pm 20\%$.

---
Figure 5:  Similar to Fig. 4 but for a different data set, where the “golden
events” correspond to 17O$(\alpha,n)$ reactions. The experimental data is from
[5].

In this case the largest disagreement is found in the classification of GE for
strips 8 and 13. As mentioned before, the classification of GE in the strips
near the edges of the detector is more challenging than away from the edges.
Therefore, the disagreement in strip 13 is not surprising. On the other hand,
we currently do not fully understand the disagreement in strip 8. However, the
cross section for that strip reported in [5] is smaller than the one obtained
in an earlier experimental work [21], as seen in Fig. 5 of Ref. [5]. A
separate investigation to see whether the prior data analysis misclassified
some events in this region could help us understand this difference, but such
effort is beyond the scope of this paper.

Table 3 shows the confusion matrix for GE and not-GE classified by the new
method and traditional methods for the 17O$(\alpha,n)$ case. In the table, we
show the mean number of events with standard deviation in parentheses, which
are evaluated by executing the new method for 50 repetitions with different
seeds. The values are rounded to the nearest integer.

As expected, the classification of GE in this case is less difficult than in
the 17F$(\alpha,p)$ case, because of the larger difference in proton number
between the heavy ion (20Ne) and the beam particles (17O). This can be seen by
the relatively large diagonal elements and relatively small off-diagonal
elements of the confusion matrix.

17O$(\alpha,n)$ | GE pres. meth. | Not-GE pres. meth.
---|---|---
GE trad. meth. | 2860(7) | 53(7)
Not-GE trad. meth. | 330(46) | 81004(46)
Table 3: Confusion matrix for the 17O$(\alpha,n)$ case showing the mean number
of events with standard deviation in parentheses.

###  6.3 Sources of uncertainty

In Tables 2 and 3 the standard deviation provides an estimate of how much
uncertainty is incurred by the proposed approach in extracting golden events.
There are multiple sources of this uncertainty in the four phases of the
approach. In phases 2 and 4, the models (the neural network transformation and
the classifier) comprise estimated parameters. Typically, these parameters are
estimated by using repetitive updates (using stochastic gradient descent)
starting from a random initialization point. The stochastic nature of the
repetitive updates contributes to parametric uncertainty [22].

All machine learning models are plagued with the choice of training data that
is used to estimate the parameters. Typically and in our methodology, the
training data is uniformly sampled from a larger set, therefore contributing
aleatoric uncertainty [22]. Furthermore, the choice of thresholds in phase 1,
2, and 3 also contributes to uncertainty where different choices of thresholds
governs the number of false positives (number of not-golden events classified
as golden) in the detection. In this work the parametric and aleatoric
uncertainty of the proposed approach is summarized by the standard deviations
obtained through repetition of the detection experiment.

###  6.4 Other ML methods for GE detection (17F$(\alpha,p)$)

Methods | Total GE | Correctly identified | Wrongly identified
---|---|---|---
One-class SVM [20] | 51852 | 1540 | 50312
Isolation Forest [23] | 53916 | 1194 | 52722
Local Outlier Factor [23] | 53766 | 1233 | 52533
Elliptical Envelope [24] | 51782 | 776 | 51006
Present | 9793 | 1825 | 7968
Table 4: Results with different outlier detection (GE are the outliers to be
detected) approaches for the 17F$(\alpha,p)$ data.

To establish the efficiency of the new methodology, we compare our approach
with several out-of-the-box outlier detection methods with the assumption that
GE’s are the outliers. We use one-class SVM [20], isolation forest [23], local
outlier factor [23], and elliptical envelope [24] for this purpose. First, we
remark that all of these out-of-the-box methods are computationally
prohibitive for detecting GE from $10^{9}$ (1 billion) points. Therefore, we
apply these methods to the result of phase 1 where most of the easy-to-
identify background events have been removed. The first four rows in Table 4
correspond to results of these outlier detection methods, and the last row
corresponds to the methodology presented in this paper. In the table we record
the total number of events identified by all these methods as GE the total
number of events that were correctly identified as GE, and the total number of
events that were wrongly identified. In this context, “correctly” or “wrongly”
identified GE means that the GE agree or disagree with the ones found using
the traditional method, respectively. We remark that, while each of the
methods is able to correctly identify several GE’s, the number of GE’s
identified by our approach is much greater than that identified by all the
out-of-the-box methods. For example, our approach correctly detects a total of
$1825$ GE’s whereas one-class SVM detects only $1540$—which is the most number
of events correctly detected by any out-of-the-box outlier detection method.
While $1540$ is a reasonable number, the main drawback of one-class SVM and
other outlier detection methods is the number of false positives. The one-
class SVM misidentifies 50312 events as GE (which is the minimum among rows in
Table 4). In contrast to these approaches, our approach wrongly identifies
only $7968$ events as GE—a reduction by a factor of 7 in the number of false
positives relative to the best anomaly detection approaches in the literature.
The additional steps utilized in our approach not only improves the detection
of GE but dramatically reduces the number of false positive—a strength of our
methodology.

##  7 Conclusion

We have developed a method for the strip-wise classification of specific
$\alpha$-induced reactions taking place inside the MUSIC detector. The method
consists of four phases where statistical and machine-learning methods are
combined.

We applied the new method to two experimental data sets, where the golden
events (events of interest) correspond to 17F$(\alpha,p)^{20}$Ne and
17O$(\alpha,n)^{20}$Ne reactions. For most strips, the new method yields
results such that the percent differences with respect to results found with
traditional methods are within $\pm 20\%$. In the strips closest to the ends
of the detector, strip 1 and strip 13 or higher, we observe the strongest
disagreements and thus conclude that the new method is not yet valid for those
strips. That said, we found that our newly developed method outperforms
several out-of-the-box methods for outlier detection.

We expect that our method will serve as a foundation for the development of
future AI/ML methods for the classification of MUSIC data with enhanced
automation.

## References

  * [1] M. Kuchera, R. Ramanujan, J. Taylor, R. Strauss, D. Bazin, J. Bradt, R. Chen, Machine learning methods for track classification in the AT-TPC, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 940 (2019) 156–167.  doi:https://doi.org/10.1016/j.nima.2019.05.097.
URL http://www.sciencedirect.com/science/article/pii/S0168900219308046

  * [2] K. Raghavan, P. Balaprakash, A. Lovato, N. Rocco, S. M. Wild, Machine-learning-based inversion of nuclear responses, Phys. Rev. C 103 (2021) 035502.  doi:10.1103/PhysRevC.103.035502.
URL https://link.aps.org/doi/10.1103/PhysRevC.103.035502

  * [3] A. Boehnlein, M. Diefenthaler, C. Fanelli, M. Hjorth-Jensen, T. Horn, M. P. Kuchera, D. Lee, W. Nazarewicz, K. Orginos, P. Ostroumov, L.-G. Pang, A. Poon, N. Sato, M. Schram, A. Scheinker, M. S. Smith, X.-N. Wang, V. Ziegler, Artificial intelligence and machine learning in nuclear physics (2021).  arXiv:2112.02309.
  * [4] P. F. F. Carnelli, S. Almaraz-Calderon, K. E. Rehm, M. Albers, M. Alcorta, P. F. Bertone, B. Digiovine, H. Esbensen, J. O. F. Niello, et al., Multi-Sampling Ionization Chamber (MUSIC) for measurements of fusion reactions with radioactive beams, Nucl. Inst. and Meth. in Phys. Res. Section A 799 (2015) 197–202.  doi:10.1016/j.nima.2015.07.030.
URL http://www.sciencedirect.com/science/article/pii/S0168900215008591

  * [5] M. Avila, K. Rehm, S. Almaraz-Calderon, A. Ayangeakaa, C. Dickerson, C. Hoffman, C. Jiang, B. Kay, J. Lai, O. Nusair, et al., Study of $(\alpha,p)$ and $(\alpha,n)$ reactions with a MUlti-sampling ionization chamber, Nucl. Inst. and Meth. in Phys. Res. Section A 859 (2017) 63–68.  doi:https://doi.org/10.1016/j.nima.2017.03.060.
URL
https://www.sciencedirect.com/science/article/pii/S0168900217304187?via%3Dihub

  * [6] P. F. F. Carnelli, S. Almaraz-Calderon, K. E. Rehm, M. Albers, M. Alcorta, P. F. Bertone, B. Digiovine, H. Esbensen, J. O. F. Niello, et al., Measurements of fusion reactions of low-intensity radioactive carbon beams on ${}^{12}\mathrm{C}$ and their implications for the understanding of X-ray bursts, Phys. Rev. Lett. 112 (2014) 192701.  doi:10.1103/PhysRevLett.112.192701.
URL https://link.aps.org/doi/10.1103/PhysRevLett.112.192701

  * [7] M. Avila, K. Rehm, S. Almaraz-Calderon, A. Ayangeakaa, C. Dickerson, C. Hoffman, C. Jiang, B. Kay, J. Lai, O. Nusair, et al., Study of the astrophysically important 23Na$(\alpha,p)^{26}$Mg and 23Na$(\alpha,n)^{26}$Al reactions, Phys. Rev. C 94 (2016) 065804.  doi:10.1103/PhysRevC.94.065804.
URL https://journals.aps.org/prc/abstract/10.1103/PhysRevC.94.065804

  * [8] R. Talwar, M. J. Bojazi, P. Mohr, K. Auranen, M. L. Avila, A. D. Ayangeakaa, J. Harker, C. R. Hoffman, C. L. Jiang, S. A. Kuvin, et al., Experimental study of ${}^{38}\mathrm{Ar}+\alpha$ reaction cross sections relevant to the ${}^{41}\mathrm{Ca}$ abundance in the solar system, Phys. Rev. C 97 (2018) 055801.  doi:10.1103/PhysRevC.97.055801.
URL https://link.aps.org/doi/10.1103/PhysRevC.97.055801

  * [9] M. L. Avila, S. Almaraz-Calderon, A. D. Ayangeakaa, C. Dickerson, C. R. Hoffman, C. L. Jiang, B. P. Kay, J. Lai, O. Nusair, et al., Study of $\alpha$-particle induced reactions using the MUSIC detector, AIP Conference Proceedings 2038 (2018) 020018.
  * [10] H. Bethe, Zur Theorie des durchgangs schneller Korpuskularstrahlen durch Materie, Annalen der Physik 397 (3) (1930) 325–400.
  * [11] D. Santiago-Gonzalez, musicsim: an open-source program to simulate events in the MUSIC detector, https://gitlab.phy.anl.gov/music/sim (2020).
  * [12] J. F. Ziegler, M. Ziegler, J. Biersack, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 268 (11) (2010) 1818–1823, 19th International Conference on Ion Beam Analysis.  doi:https://doi.org/10.1016/j.nimb.2010.02.091, [link].
URL https://www.sciencedirect.com/science/article/pii/S0168583X10001862

  * [13] H. Weick, H. Geissel, N. Iwasa, C. Scheidenberger, J. L. Rodríguez-Sánchez, A. Prochazka, S. Purushotaman (10 2018).
  * [14] A. Soylemezoglu, S. Jagannathan, C. Saygin, Mahalanobis–Taguchi system as a multi-sensor based decision making prognostics tool for centrifugal pump failures, IEEE Transactions on Reliability 60 (4) (2011) 864–878.
  * [15] R. A. Johnson, D. W. Wichern, et al., Applied multivariate statistical analysis, Vol. 6, Pearson London, UK:, 2014.
  * [16] K. Khan, S. U. Rehman, K. Aziz, S. Fong, S. Sarasvady, Dbscan: Past, present and future, in: The fifth international conference on the applications of digital information and web technologies (ICADIWT 2014), IEEE, 2014, pp. 232–238.
  * [17] M. Hardt, B. Recht, Y. Singer, Train faster, generalize better: Stability of stochastic gradient descent, in: International Conference on Machine Learning, PMLR, 2016, pp. 1225–1234.
  * [18] P. Ram, K. Sinha, Revisiting kd-tree for nearest neighbor search, in: Proceedings of the 25th acm sigkdd international conference on knowledge discovery & data mining, 2019, pp. 1378–1388.
  * [19] C. M. Bishop, Pattern recognition and machine learning, Springer, 2006.
  * [20] B. Schölkopf, R. C. Williamson, A. Smola, J. Shawe-Taylor, J. Platt, Support vector method for novelty detection, Advances in neural information processing systems 12 (1999).
  * [21] J. K. Bair, F. X. Haas, Total neutron yield from the reactions ${}^{13}\mathrm{C}(\alpha,n)^{16}\mathrm{O}$ and ${}^{17,18}\mathrm{O}(\alpha,n)^{20,21}\mathrm{Ne}$, Phys. Rev. C 7 (1973) 1356–1364.  doi:10.1103/PhysRevC.7.1356.
URL https://link.aps.org/doi/10.1103/PhysRevC.7.1356

  * [22] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al., TensorFlow: Large-scale machine learning on heterogeneous distributed systems, arXiv preprint arXiv:1603.04467 (2016).
  * [23] Z. Cheng, C. Zou, J. Dong, Outlier detection using isolation forest and local outlier factor, in: Proceedings of the conference on research in adaptive and convergent systems, 2019, pp. 161–168.
  * [24] M. Antonini, M. Vecchio, F. Antonelli, P. Ducange, C. Perera, Smart audio sensors in the internet of things edge for anomaly detection, IEEE Access 6 (2018) 67594–67610.


11institutetext: Electrical Engineering Faculty, Częstochowa University of
Technology,
Częstochowa, Poland
11email<EMAIL_ADDRESS>
# Ensembles of Randomized Neural Networks for Pattern-based Time Series
Forecasting††thanks: Supported by Grant 2017/27/B/ST6/01804 from the National
Science Centre, Poland.

Grzegorz Dudek  0000-0002-2285-0327 Paweł Pełka  0000-0002-2609-811X

###### Abstract

In this work, we propose an ensemble forecasting approach based on randomized
neural networks. Improved randomized learning streamlines the fitting
abilities of individual learners by generating network parameters in
accordance with the data and target function features. A pattern-based
representation of time series makes the proposed approach suitable for
forecasting time series with multiple seasonality. We propose six strategies
for controlling the diversity of ensemble members. Case studies conducted on
four real-world forecasting problems verified the effectiveness and superior
performance of the proposed ensemble forecasting approach. It outperformed
statistical models as well as state-of-the-art machine learning models in
terms of forecasting accuracy. The proposed approach has several advantages:
fast and easy training, simple architecture, ease of implementation, high
accuracy and the ability to deal with nonstationarity and multiple seasonality
in time series.

###### Keywords:

Ensemble forecasting Pattern representation of time series Randomized neural
networks Short-term load forecasting Time series forecasting.

##  1 Introduction

Time series (TS) forecasting plays an important role in many fields, including
commerce, industry, public administration, politics, health, medicine, etc.
[1]. TS expressing different phenomena and processes may include the nonlinear
trend, multiple seasonal cycles of different lengths and random fluctuations.
This makes the relationship between predictors and output variables very
complex and places high demands on forecasting models. Over the years, many
sophisticated forecasting models have been proposed including statistical,
machine learning (ML) and hybrid solutions.

Among ML models, neural networks (NNs) are the most commonly used. There are a
huge number of forecasting models based on different NN architectures [2]. In
addition to classic NNs such as multilayer perceptron (MLP), radial basis
function NN, generalized regression NN (GRNN), and self-organizing maps [3],
many models based on deep learning (DL) have been developed recently. These
are composed of combinations of basic structures, such as MLPs, convolutional
NNs, and recurrent NNs (RNNs) [4], [5]. Their success can be largely
attributed to increased model complexity and the ability to perform
representation learning and cross-learning.

An effective way to increase the performance of a predictive model is
ensembling. Ensemble methods combine in some fashion multiple learning
algorithms to produce a common response, hopefully improving accuracy and
stability compared to a single learner. The challenge in ensemble learning,
which to a large extent determines its success, is achieving a good tradeoff
between the performance and the diversity of the ensemble members [6]. The
more accurate and the more diverse the individual learners are, the better the
ensemble performance. Depending on the base model type, the diversity of
learners can be achieved using different strategies. For example, the winning
submission to the M4 Makridakis forecasting competition, was a hybrid model
combining exponential smoothing (ETS) and long-short term memory [7], which
uses three sources of diversity. The first is a stochastic training process,
the second is similar to bagging, and the third is training ensemble members
using different initial parameters. In another state-of-the-art DL forecasting
model, N-Beats [8], in order to achieve diversity, each of the ensemble
members is trained using a different random initialization and a different
random sequence of batches.

Randomization-based NNs are especially suitable for ensembling as they are
highly unstable and extremely fast trained [9]. In randomized NNs, the hidden
node parameters are selected randomly and the output weights are tuned using a
simple and fast least-squares method. This avoids the difficulties associated
with gradient-based learning such as slow convergence, local minima problems,
and model uncertainties caused by the initial parameters. Ensemble methods
based on randomized NNs were presented recently in several papers. In [10], an
ensemble based on decorrelated random vector functional link (RVFL) networks
was proposed using negative correlation learning to control the trade-off
among the bias, variance and covariance in the ensemble learning. A selective
ensemble of randomization-based NNs using successive projections algorithm
(SPA) was put forward in [11]. SPA improves diversity by selecting
uncorrelated members. In [12], to enhance the generalization capacities of the
ensemble, a novel framework for building an ensemble model by selecting
appropriate representatives from a number of randomization-base models (RVFL
or stochastic configuration networks) was proposed.

In the forecasting domain, ensembles of randomization-based NNs are widely
used. Examples include: [13], where ensemble members, which are extreme
learning machines, learn on TS decomposed by wavelet transform; [14], where a
hybrid incremental learning approach is presented for ensemble forecasting,
which is composed of discrete wavelet transform, empirical mode decomposition
and RVFL networks as base learners; [15], where a new bagging ensemble
approach based on NN with random weights for online data stream regression is
proposed; [16], where a data-driven evolutionary ensemble forecasting model is
proposed using two optimization objectives (error and diversity) and RVFL as a
base learner.

Motivated by the good performance of the randomized NN based ensemble
forecasting models mentioned above, and new improvements in randomized
learning [17] as well as pattern-based TS representation suitable for TS with
multiple seasonality [3], this work contributes to the development of
forecasting models in the following ways:

  1. 1.

A new ensemble forecasting model based on randomized NNs is presented.
Enhanced randomized learning improves the fitting abilities of individual
learners. To deal with multiple seasonality and nonstationarity, the model
applies pattern representation of TS.

  2. 2.

Six strategies for generating ensemble diversity are proposed. Each strategy
governs diversity in a different way.

  3. 3.

An experimental study using four real-world datasets demonstrates the superior
performance of the proposed ensemble forecasting approach when compared to
statistical and machine learning forecasting models.

The remaining sections of this work are organized as follows. Section 2
presents the base forecasting model, $\operatorname{\textsc{randNN}}$. Section
3 describes strategies for generating ensemble diversity. The performance of
the proposed ensemble forecasting model is evaluated in Section 4. Concluding
remarks are given in Section 5.

##  2 Forecasting Model

Fig. 1 shows $\operatorname{\textsc{randNN}}$, a base forecasting model, which
was designed for forecasting TS with multiple seasonality [18]. It is used as
an ensemble member. $\operatorname{\textsc{randNN}}$ is composed of an
encoder, a decoder and a randomized feedforward NN (FNN).

Figure 1: Block diagram of the proposed forecasting model.

The encoder transforms the original TS into unified input and output patterns
of its seasonal cycles. To create input patterns, the TS expressing multiple
seasonality, $\\{E_{k}\\}_{k=1}^{K}$, is divided into seasonal sequences of
the shortest length. Let these sequences be expressed by vectors
$\mathbf{e}_{i}=[E_{i,1},E_{i,2},…,E_{i,n}]^{T}$, where $n$ is the seasonal
sequence length and $i=1,2,...,K/n$ is the sequence number. These sequences
are encoded in input patterns $\mathbf{x}_{i}=[x_{i,1},x_{i,2},…,x_{i,n}]^{T}$
as follows:

$\mathbf{x}_{i}=\frac{\mathbf{e}_{i}-\overline{e}_{i}}{\widetilde{e}_{i}}$ (1)

where $\overline{{e}}_{i}$ is the mean value of sequence $\mathbf{e}_{i}$, and
$\widetilde{e}_{i}=\sqrt{\sum_{t=1}^{n}(E_{i,t}-\overline{e}_{i})^{2}}$ is a
measure of sequence $\mathbf{e}_{i}$ dispersion.

Note that the x-patterns are normalized versions of centered vectors
$\mathbf{e}_{i}$. All x-patterns, representing successive seasonal sequences,
have a zero mean, the same variance and the same unity length. However, they
differ in shape. Thus, the original seasonal sequences, which have a different
mean value and dispersion, are unified (see Fig. 2 in [18]).

The output patterns $\mathbf{y}_{i}=[y_{i,1},y_{i,2},…,y_{i,n}]^{T}$ represent
the forecasted sequences
$\mathbf{e}_{i+\tau}=[E_{i+\tau,1},E_{i+\tau,2},…,E_{i+\tau,n}]^{T}$, where
$\tau\geq 1$ is a forecast horizon. The y-patterns are determined as follows:

$\mathbf{y}_{i}=\frac{\mathbf{e}_{i+\tau}-\overline{e}_{i}}{\widetilde{e}_{i}}$
(2)

where $\overline{{e}}_{i}$ and $\widetilde{e}_{i}$ are the same as in (1).

Note that in (2), for the $i$-th output pattern, we use the same coding
variables $\overline{{e}}_{i}$ and $\widetilde{e}_{i}$ as for the $i$-th input
pattern. This is because the coding variables for the forecasted sequence,
$\overline{e}_{i+\tau}$ and $\widetilde{e}_{i+\tau}$, are unknown for the
future period.

The decoder transforms a forecasted output pattern into a TS seasonal cycle.
The output pattern predicted by randomized FNN is decoded using the coding
variables of the input query pattern, $\mathbf{x}$, using transformed equation
(2):

$\widehat{\mathbf{e}}=\widehat{\mathbf{y}}\widetilde{{e}}+\overline{{e}}$ (3)

where $\widehat{\mathbf{e}}$ is the forecasted seasonal sequence,
$\widehat{\mathbf{y}}$ is the forecasted output pattern, $\widetilde{e}$ and
$\overline{e}$ are the coding variables determined from the TS sequence
encoded in query pattern $\mathbf{x}$.

The randomized FNN is composed of $n$ inputs, one hidden layer with $m$
nonlinear nodes, and $n$ outputs. Logistic sigmoid activation functions are
employed for hidden nodes. The training set is composed of the corresponding
input and output patterns:
$\Phi=\left\\{(\mathbf{x}_{i},\mathbf{y}_{i})|\mathbf{x}_{i},\mathbf{y}_{i}\in\mathbb{R}^{n},i=1,2,\ldots,N\right\\}$.
The randomized learning algorithm consists of three steps (we use an improved
version of the randomized learning algorithm proposed in [17]):

  1. 1.

Generate hidden node parameters. The weights are selected randomly:
$a_{t,j}=\sim U(-u,u)$ and biases are calculated from:

$b_{j}=-\mathbf{a}_{j}^{T}\mathbf{x}^{*}_{j}$ (4)

where $j=1,2,...,m$; $t=1,2,...,n$; $\mathbf{x}^{*}_{j}$ is one of the
training x-patterns selected for the $j$-th hidden node at random.

  2. 2.

Calculate hidden layer output matrix $\mathbf{H}$.

  3. 3.

Calculate the output weights:

$\boldsymbol{\beta}=\mathbf{H}^{+}\mathbf{Y}$ (5)

where $\boldsymbol{\beta}\in\mathbb{R}^{m\times n}$ is a matrix of output
weights, $\mathbf{Y}\in\mathbb{R}^{N\times n}$ is a matrix of target output
patterns, and $\mathbf{H}^{+}\in\mathbb{R}^{m\times N}$ is the Moore-Penrose
generalized inverse of matrix $\mathbf{H}$.

In the first step, the weights are selected randomly from symmetrical interval
$[-u,u]$. The interval bounds, $u$, decide about the steepness of the
sigmoids. To make the bounds interpretable, let us express them by the sigmoid
slope angle [17]: $u=4\tan\alpha_{max}$, where $\alpha_{max}$ is the upper
bound for slope angles. Hyperparameter $\alpha_{max}$ as well as number of
hidden nodes $m$ control the bias-variance tradeoff of the model. Both these
hyperparameters should be adjusted to the target function complexity.

The biases calculated according to (4) ensure the introduction of the steepest
fragments of the sigmoid into the input hypercube [17]. These fragments are
most useful for modeling target function fluctuations. Thus, there are no
saturated sigmoids and wasted nodes in the network.

##  3 Ensembling

An ensemble is composed of $M$ individual learners
($\operatorname{\textsc{randNN}s}$). Each ensemble member learns from the
training set $\Phi$ using one of the strategies for generating diversity,
$\operatorname{\textsc{ens1}}-\operatorname{\textsc{ens6}}$ described below.
The ensemble prediction is an average of individual member predictions
$\mathbf{\widehat{e}}_{k}$:

$\mathbf{\widehat{e}}_{\text{ens}}=\frac{1}{M}\sum_{k=1}^{M}\mathbf{\widehat{e}}_{k}$
(6)

As an ensemble diversity measure, we define the average standard deviation of
forecasts produced by the individual learners:

$Diversity=\frac{1}{n|\Psi|}\sum_{i\in\Psi}\sum_{t=1}^{n}\sqrt{\frac{1}{M}\sum_{k=1}^{M}(\widehat{E}^{k}_{i,t}-\overline{\widehat{E}}_{i,t})^{2}}$
(7)

where $\Psi$ is a test set, $\widehat{E}^{k}_{i,t}$ is a forecast of the
$t$-th element of the $i$-th seasonal sequence produced by the $k$-th learner,
and $\overline{\widehat{E}}_{i,t}$ is an average of forecasts produced by $M$
learners.

An ensemble diversity is generated using one of the following six strategies:

  * $\operatorname{\textsc{ens1}}$

generates diversity by using different parameters of hidden nodes. For each
learner, new weights are randomly selected taking $\alpha_{\max}$ as the upper
bound for the sigmoid slope angles. Then, biases are calculated from (4) and
output weights form (5). The diversity level is controlled by $\alpha_{\max}$.
For larger $\alpha_{\max}$ we get steeper sigmoids and higher diversity.

  * $\operatorname{\textsc{ens2}}$

controls diversity by training individual learners on different subsets of the
training set. For each ensemble member, a random sample from the training set
is selected without replacement. The sample size is $N^{\prime}=\eta N$, where
$\eta\in(0,1)$ is a diversity parameter. Each learner has the same hidden node
parameters. Its output weights are tuned to the training subset.

  * $\operatorname{\textsc{ens3}}$

controls diversity by training individual learners on different subsets of
features. For each ensemble member the features are randomly sampled without
replacement. The sample size is $n^{\prime}=\kappa n$, where $\kappa\in(0,1)$
is a diversity parameter. The ensemble members share the hidden node
parameters. Their output weights are tuned to the training set.

  * $\operatorname{\textsc{ens4}}$

is based on hidden node pruning. The learners are created from the initial
$\operatorname{\textsc{randNN}}$ architecture including $m$ hidden nodes. For
each learner, $m^{\prime}=\rho m$ nodes are randomly selected and the
remaining are pruned. Output weights $\beta$ are determined anew for each
learner. Parameter $\rho\in(0,1)$ controls the diversity level.

  * $\operatorname{\textsc{ens5}}$

is based on hidden weight pruning. The learners are created from the initial
$\operatorname{\textsc{randNN}}$ architecture including $m$ hidden nodes. For
each learner, $p=\lambda mn$ hidden node weights are randomly selected and set
to zero. Output weights $\beta$ are determined anew for each learner.
Parameter $\lambda\in(0,1)$ controls the diversity level.

  * $\operatorname{\textsc{ens6}}$

generates diversity by noising training data. For each learner the training
patterns are perturbed by Gaussian noise as follows:
$x_{i,t}=x_{i,t}(1+\zeta_{i,t}),y_{i,t}=y_{i,t}(1+\xi_{i,t})$, where
$\zeta_{i,t},\xi_{i,t}\sim N(0,\sigma)$. Standard deviation of the noise,
$\sigma$, is a diversity parameter. Each learner has the same hidden node
parameters. Its output weights are tuned to the noised training data.

##  4 Experiments and Results

We evaluate the performance of the proposed ensembles of
$\operatorname{\textsc{randNN}s}$ on four real-world forecasting problems.
These concern forecasting electricity demand for four European countries:
Poland (PL), Great Britain (GB), France (FR) and Germany (DE) (data was
collected from www.entsoe.eu). The hourly electrical load TS express three
seasonalities: yearly, weekly and daily (see Fig. 2 in [18]). The data period
covers the 4 years from 2012 to 2015. Atypical days such as public holidays
were excluded from the data (between 10 and 20 days a year). The forecast
horizon $\tau$ is one day, i.e. 24 hours. We forecast the daily load profile
for each day of 2015. For each forecasted day, a new training set is created
which includes historical pairs of corresponding input and output patterns
representing the same days of the week as the pair in the query pattern and
forecasted pattern. The number of ensemble members was $M=100$.

In the first experiment, to assess the impact of the
$\operatorname{\textsc{randNN}s}$ hyperparameters on the forecasting accuracy
of $\operatorname{\textsc{ens1}}$, we train the
$\operatorname{\textsc{randNN}}$ members with $m=10,20,...,70$ and
$\alpha_{\max}=0^{\circ},10^{\circ},...,90^{\circ}$. Fig. 2 shows the mean
absolute percentage errors (MAPE) depending on the hyperparameters. The
optimal values of $m/\alpha_{\max}$ were: $50/70^{\circ}$ for PL,
$30/60^{\circ}$ for GB, $40/70^{\circ}$ for FR, and $40/80^{\circ}$ for DE.
Based on these results, we select $m=40$ and $\alpha_{\max}=70^{\circ}$ for
each ensemble variant and dataset, except $\operatorname{\textsc{ens1}}$,
where we control diversity by changing $\alpha_{\max}$, and
$\operatorname{\textsc{ens4}}$, where we control diversity by changing the
number of pruning nodes. In this latter case we use $m=80$ nodes.

From Fig. 3, we can assess the impact of the hyperparameters on the diversity
(7) of $\operatorname{\textsc{ens1}}$. As expected, greater numbers of hidden
nodes and steeper sigmoids (higher $\alpha_{\max}$) cause an increase in
diversity.

Figure 2: MAPE depending on $\operatorname{\textsc{randNN}}$ hyperparameters
for $\operatorname{\textsc{ens1}}$. Figure 3: Diversity of
$\operatorname{\textsc{ens1}}$ depending on $\operatorname{\textsc{randNN}}$
hyperparameters.

Fig. 4 shows MAPE and ensemble diversity depending on the diversity parameters
for all ensemble variants and datasets. From this figure, we can conclude
that:

  * •

for $\operatorname{\textsc{ens1}}$ diversity increases with $\alpha_{\max}$.
The optimal $\alpha_{\max}$ is around $70^{\circ}-80^{\circ}$ for all
datasets.

  * •

for $\operatorname{\textsc{ens2}}$ diversity decreases with $\eta$. For lower
$\eta$ (small training sets), the diversity as well as MAPE are very high. The
MAPE curves have an irregular character, so it is difficult to select the
optimal $\eta$ value. The error levels are much higher than for
$\operatorname{\textsc{ens1}}$.

  * •

for $\operatorname{\textsc{ens3}}$ diversity decreases with $\kappa$. MAPE
reaches its minima for: $\kappa\in(8/24,16/24)$ for PL, $\kappa=20/24$ for GB,
$\kappa\in(14/24,18/24)$ for FR, and $\kappa\in(10/24,18/24)$ for DE. Thus,
the best ensemble solutions use around $40^{\circ}-84^{\circ}$ of features.

  * •

for $\operatorname{\textsc{ens4}}$ diversity changes slightly with $\rho$.
MAPE reaches its minima for: $\rho\in(40/80,48/80)$ for PL, $\rho=32/80$ for
GB, $\rho\in(40/80,48/80)$ for FR, and $\rho=40/80)$ for DE. Thus, the best
choice for the number of selected nodes is around $m^{\prime}=40$.

  * •

for $\operatorname{\textsc{ens5}}$ diversity has its maximum for $\lambda=0.1$
(10% of hidden weights are set to 0). MAPE increases with $\lambda$, having
its lowest values for $\lambda=0$ (no weight pruning). Thus, any weight
pruning makes the error higher. Increasing the hidden node number to 80 did
not improve the results.

  * •

for $\operatorname{\textsc{ens6}}$ we observe low sensitivity of diversity as
well as MAPE to diversity parameter $\sigma$. For smaller $\sigma$, the error
only slightly decreases as $\sigma$ increases. Then, for higher $\sigma$, the
error starts to increase.

Figure 4: MAPE (solid lines) and ensemble diversity (dashed lines) depending
on the diversity parameters.

Table 1 shows the quality metrics of the forecasts for the optimal values of
the diversity parameters: MAPE, median of APE, root mean square error (RMSE),
mean percentage error (MPE), and standard deviation of percentage error
(Std(PE)) as a measure of the forecast dispersion. For comparison, the results
for a single $\operatorname{\textsc{randNN}}$ are shown. In this case, for
each forecasting task (each day of 2015) the hyperparameters of
$\operatorname{\textsc{randNN}}$, $m$ and $u$, were optimized using grid
search and 5-fold cross-validation [18]. The results for single
$\operatorname{\textsc{randNN}}$ shown in Table 1 are averaged over 100
independent training sessions. $\operatorname{\textsc{ens7}}$ shown in Table 1
is an ensemble of these 100 runs (calculated from (6)).
$\operatorname{\textsc{ens7}}$ is similar to $\operatorname{\textsc{ens1}}$.
The only difference is that for $\operatorname{\textsc{ens7}}$ the
$\operatorname{\textsc{randNN}}$ hyperparameters were optimized for each
forecasting task, and for $\operatorname{\textsc{ens1}}$ we set $m=40$ and
$\alpha_{\max}=70^{\circ}$ for all forecasting tasks. Distributions of APE are
shown in Fig. 5.

Table 1: Forecasting results. |  | $\operatorname{\textsc{ens1}}$ | $\operatorname{\textsc{ens2}}$ | $\operatorname{\textsc{ens3}}$ | $\operatorname{\textsc{ens4}}$ | $\operatorname{\textsc{ens5}}$ | $\operatorname{\textsc{ens6}}$ | $\operatorname{\textsc{randNN}}$ | $\operatorname{\textsc{ens7}}$
---|---|---|---|---|---|---|---|---|---
PL | MAPE | 1.14 | 1.32 | 1.13 | 1.16 | 1.39 | 1.34 | 1.32 | 1.24
| Median(APE) | 0.79 | 0.94 | 0.78 | 0.81 | 0.98 | 0.94 | 0.93 | 0.89
| RMSE | 304 | 342 | 299 | 313 | 380 | 351 | 358 | 333
| MPE | 0.31 | 0.43 | 0.31 | 0.25 | 0.36 | 0.45 | 0.40 | 0.40
| Std(PE) | 1.67 | 1.88 | 1.64 | 1.72 | 2.07 | 1.94 | 1.94 | 1.80
GB | MAPE | 2.51 | 2.74 | 2.55 | 2.49 | 2.92 | 2.65 | 2.61 | 2.52
| Median(APE) | 1.76 | 2.01 | 1.82 | 1.77 | 2.09 | 1.95 | 1.88 | 1.80
| RMSE | 1151 | 1205 | 1160 | 1131 | 1325 | 1167 | 1187 | 1147
| MPE | -0.53 | -0.83 | -0.53 | -0.53 | -0.55 | -0.65 | -0.61 | -0.61
| Std(PE) | 3.48 | 3.64 | 3.51 | 3.42 | 4.02 | 3.56 | 3.57 | 3.44
FR | MAPE | 1.57 | 1.80 | 1.59 | 1.61 | 1.90 | 1.74 | 1.67 | 1.60
| Median(APE) | 1.04 | 1.24 | 1.07 | 1.06 | 1.33 | 1.24 | 1.15 | 1.07
| RMSE | 1378 | 1505 | 1398 | 1402 | 1597 | 1459 | 1422 | 1385
| MPE | -0.28 | -0.26 | -0.25 | -0.31 | -0.33 | -0.27 | -0.42 | -0.42
| Std(PE) | 2.53 | 2.80 | 2.58 | 2.57 | 2.92 | 2.72 | 2.60 | 2.53
DE | MAPE | 1.18 | 1.31 | 1.19 | 1.20 | 1.45 | 1.33 | 1.38 | 1.29
| Median(APE) | 0.81 | 0.91 | 0.80 | 0.82 | 0.98 | 0.96 | 0.96 | 0.91
| RMSE | 1077 | 1155 | 1111 | 1097 | 1423 | 1159 | 1281 | 1206
| MPE | 0.11 | 0.03 | 0.10 | 0.09 | 0.02 | 0.17 | 0.14 | 0.14
| Std(PE) | 1.89 | 2.01 | 1.96 | 1.92 | 2.44 | 2.04 | 2.22 | 2.11

Figure 5: Distribution of APE.

To confirm that the differences in accuracy between ensemble variants are
statistically significant, we performed a one-sided Giacomini-White test for
conditional predictive ability [19]. This is a pairwise test that compares the
forecasts produced by different models. We used a Python implementation of the
Giacomini-White test in multivariate variant from [20, 21]. Fig. 6 shows
results of the Giacomini-White test. The resulting plots are heat maps
representing the obtained $p$-values. The closer they are to zero the
significantly more accurate the forecasts produced by the model on the
$X$-axis are than the forecasts produced by the model on the $Y$-axis. The
black color indicates $p$-values larger than 0.10.

Figure 6: Results of the Giacomini-White test for the proposed
$\operatorname{\textsc{randNN}}$ ensembles.

From the results presented in Table 1 and Figs. 5 and 6 we can conclude that
the most accurate ensembles are $\operatorname{\textsc{ens1}}$,
$\operatorname{\textsc{ens3}}$, and $\operatorname{\textsc{ens4}}$. For PL and
DE they significantly outperformed all other models including single
$\operatorname{\textsc{randNN}}$ and $\operatorname{\textsc{ens7}}$, although,
for GB and FR, $\operatorname{\textsc{ens7}}$ can compete with them in terms
of accuracy. Note that the APE distributions for
$\operatorname{\textsc{ens1}}$, $\operatorname{\textsc{ens3}}$, and
$\operatorname{\textsc{ens4}}$ shown in a simplified form in Fig. 5 are very
similar. The worst ensemble solution was $\operatorname{\textsc{ens5}}$ based
on weight pruning. It was beaten by all other models for all datasets.
$\operatorname{\textsc{ens2}}$ and $\operatorname{\textsc{ens6}}$ were only
slightly better in terms of accuracy than $\operatorname{\textsc{ens5}}$.

MPE shown in Table 1 allows us to assess the bias of the forecasts. A positive
value of MPE indicates underprediction, while a negative value indicates
overprediction. As can be seen from Table 1, for PL and DE all the models
underpredicted, whilst for GB and FR they overpredicted. The forecasts
produced by the three best ensemble models were less biased than the forecast
produced by a single $\operatorname{\textsc{randNN}}$ and
$\operatorname{\textsc{ens7}}$.

In the next experiment, we compare $\operatorname{\textsc{ens1}}$ performance
(one of the best ensemble solutions) with that of other models based on
classical statistical methods and ML methods. The baseline models that we use
in our comparative studies are outlined below (see [22] for details). Their
hyperparameters were selected on the training set in grid search procedures.

  * •

naive – naive model: $\widehat{\mathbf{e}}_{i+\tau}=\mathbf{e}_{i+\tau-7}$,

  * •

arima – autoregressive integrated moving average model,

  * •

ets – exponential smoothing model,

  * •

prophet – a modular additive regression model with nonlinear trend and
seasonal components [23],

  * •

mlp – perceptron with a single hidden layer and sigmoid nonlinearities,

  * •

svm – linear epsilon insensitive support vector machine ($\epsilon$-SVM) [24],

  * •

anfis – adaptive neuro-fuzzy inference system,

  * •

lstm – long short-term memory,

  * •

fnm – fuzzy neighborhood model,

  * •

n-we – Nadaraya–Watson estimator,

  * •

grnn – general regression NN.

Table 2 shows MAPE for $\operatorname{\textsc{ens1}}$ and the baseline models.
Note that for each dataset, $\operatorname{\textsc{ens1}}$ returned the lowest
MAPE. To confirm the statistical significance of these findings, we performed
a Giacomini-White test. The test results depicted in Fig. 7, clearly show that
$\operatorname{\textsc{ens1}}$ outperforms all the other models in terms of
accuracy. Only in two cases out of 44 were the baseline models close to
$\operatorname{\textsc{ens1}}$ in terms of accuracy, i.e. N-WE for PL data and
SVM for GB data.

Table 2: MAPE for $\operatorname{\textsc{ens1}}$ and baseline models. | $\operatorname{\textsc{ens1}}$ | naive | arima | ets | prophet | mlp | svm | anfis | lstm | fnm | n-we | grnn
---|---|---|---|---|---|---|---|---|---|---|---|---
PL | 1.14 | 2.96 | 2.31 | 2.14 | 2.63 | 1.39 | 1.32 | 1.64 | 1.57 | 1.21 | 1.19 | 1.22
GB | 2.51 | 4.80 | 3.50 | 3.19 | 4.00 | 2.84 | 2.54 | 2.80 | 2.92 | 3.02 | 3.12 | 3.01
FR | 1.57 | 5.53 | 3.00 | 2.79 | 4.71 | 1.93 | 1.63 | 2.12 | 1.81 | 1.84 | 1.86 | 1.81
DE | 1.18 | 3.13 | 2.31 | 2.10 | 3.23 | 1.58 | 1.38 | 2.48 | 1.57 | 1.30 | 1.29 | 1.30

Figure 7: Results of the Giacomini-White test for the proposed
$\operatorname{\textsc{ens1}}$ and baseline models.

##  5 Conclusion

Challenging forecasting problems, such as forecasting nonstationary TS with
multiple seasonality, need sophisticated models. In this work, to deal with
multiple seasonal periods, we employ a pattern representation of the TS, in
order to simplify the relationship between input and output data and make the
problem easier to solve using simple regression models such as randomized NNs.
$\operatorname{\textsc{randNN}s}$ have three distinct advantages: a simple
single-hidden layer architecture, extremely fast learning (no gradient-based
learning, no problems with vanishing and exploding gradients), and ease of
implementation (no complex optimization algorithms, no additional mechanisms
such as dilation, attention, and residual connections).
$\operatorname{\textsc{randNN}}$ produces a vector output, which means the
model is able to forecast the time series sequence at once.

We propose an ensemble of $\operatorname{\textsc{randNN}s}$ with six different
methods of managing diversity. Among these, the most promising ones turned out
to be: learning the ensemble members using different parameters of hidden
nodes ($\operatorname{\textsc{ens1}}$), training individual learners on
different subsets of features ($\operatorname{\textsc{ens3}}$), and hidden
node pruning ($\operatorname{\textsc{ens4}}$). These three methods brought
comparable results in an experimental study involving four real-world
forecasting problems (short-term load forecasting).
$\operatorname{\textsc{ens1}}$ has the additional advantage of having only two
hyperparameters to tune, i.e. the number of hidden nodes and the interval
bounds for random weights. In the comparative study,
$\operatorname{\textsc{ens1}}$ significantly outperformed both classical
statistical methods and machine learning methods.

In our further work, we plan to introduce an attention mechanism into our
randomization-based forecasting models to select training data and develop
probabilistic forecasting models based on $\operatorname{\textsc{randNN}s}$.

## References

  * [1] Palit, A. K., Popovic, D.: Computational Intelligence in Time Series Forecasting: Theory and Engineering Applications. Springer-Verlag, London 2005.
  * [2] Benidis, K., Rangapuram, S.S., Flunkert, V., Wang, B., Maddix, D., Turkmen, C., Gasthaus, J., Bohlke-Schneider, M., Salinas, D., Stella, L., Callot L. and Januschowski, T.: Neural forecasting: Introduction and literature overview. arXiv:2004.10240 (2020)
  * [3] Dudek, G.: Neural networks for pattern-based short-term load forecasting: A comparative study. Neurocomputing 205, 64–74 (2016)
  * [4] Torres, J.F., Hadjout, D. Sebaa, A., Martínez-Álvarez, F., Troncoso, A.: Deep Learning for Time Series Forecasting: A Survey. Big Data 9(1), 3–21 (2021)
  * [5] Hewamalage, H., Bergmeir, C., Bandara, K.: Recurrent Neural Networks for Time Series Forecasting: Current Status and Future Directions. International Journal of Forecasting 37(1), 388–427 (2021)
  * [6] Reeve H.W.J., Brown, G.: Diversity and degrees of freedom in regression ensembles. Neurocomputing 298, 55–68 (2018)
  * [7] Smyl, S.: A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting. International Journal of Forecasting 36(1), 75–85 (2020)
  * [8] Oreshkin, B.N., Carpov, D., Chapados, N. and Bengio, Y.: N-BEATS: Neural basis expansion analysis for interpretable time series forecasting. 8th International Conference on Learning Representations, ICLR, (2020).
  * [9] Ren, Y., Zhang, L., Suganthan, P.N.: Ensemble Classification and Regression – Recent Developments, Applications and Future Directions. IEEE Computational Intelligence Magazine, 11(1), 41–53 (2016)
  * [10] Alhamdoosh, M., Wang, D.: Fast decorrelated neural network ensembles with random weights. Information Sciences 264, 104–117 (2014)
  * [11] Mesquita, D.P.P., Gomes, J.P.P, Rodrigues, L.R., Oliveira, S.A.F., Galvão, R.K.H.: Building selective ensembles of Randomization Based Neural Networks with the successive projections algorithm. Applied Soft Computing 70, 1135–1145 (2018)
  * [12] Huang, C., Li, M., Wang, D.: Stochastic configuration network ensembles with selective base models. Neural Networks, 264, 106–118 (2021)
  * [13] Li, S., Goel, L., Wang, P.: An ensemble approach for short-term load forecasting by extreme learning machine. Applied Energy 170, 22–29 (2016)
  * [14] Qiu, X., Suganthan, P.N., Amaratunga, G.A.J., Ensemble incremental learning Random Vector Functional Link network for short-term electric load forecasting. Knowledge-Based Systems 145, 182–196 (2018)
  * [15] de Almeida, R., Goh, Y.M., Monfared, R. et al. An ensemble based on neural networks with random weights for online data stream regression. Soft Computing 24, 9835–9855 (2020)
  * [16] Hu, Y., Qu, B., Wang, J., Liang, J., Wang, Y., Yu, K., Li, Y., Qiao, K.: Short-term load forecasting using multimodal evolutionary algorithm and random vector functional link network based ensemble learning. Applied Energy 285, 116415 (2021)
  * [17] Dudek, G.: Generating random parameters in feedforward neural networks with random hidden nodes: Drawbacks of the standard method and how to improve it. Neural Information Processing, ICONIP 2020, Springer CCIS 1333, 598–606 (2020)
  * [18] Dudek G.: Randomized Neural Networks for Forecasting Time Series with Multiple Seasonality. 16th International Work-Conference on Artificial Neural Networks IWANN 2021 (in print).
  * [19] Giacomini, R., White, H.: Tests of conditional predictive ability. Econometrica 74(6), 1545–1578 (2006)
  * [20] Epftoolbox library - https://github.com/jeslago/epftoolbox, epftoolbox documentation - https://epftoolbox.readthedocs.io
  * [21] Lago, J. Marcjasz, G., De Schutter, B., Weron, R.: Forecasting day-ahead electricity prices: A review of state-of-the-art algorithms, best practices and an open-access benchmark. Applied Energy 293, 116983 (2021)
  * [22] Dudek, G., Pełka, P.: Pattern Similarity-based Machine Learning Methods for Mid-term Load Forecasting: A Comparative Study. Applied Soft Computing 104, 107223 (2021)
  * [23] Taylor, S.J., Letham, B. Forecasting at Scale. The American Statistician 72(1), 37–45 (2018)
  * [24] Pełka P.: Pattern-based Forecasting of Monthly Electricity Demand using Support Vector Machine. International Joint Conference on Neural Networks IJCNN 2021 (in print).