# Predicting the redshift of gamma-ray loud AGNs using Supervised Machine
Learning: Part 2

Aditya Narendra  Jagiellonian University, Poland  Spencer James Gibson
Carnegie Mellon University, USA Maria Giovanna Dainotti
<EMAIL_ADDRESS>Maria Giovanna Dainotti  National Astronomical
Observatory of Japan, Mitaka  Space Science Institute, 4750 Walnut St, Suite
205, Boulder,CO,80301,USA  Malgorzata Bogdan  Department of Mathematics,
University of Wroclaw, Poland  Department of Statistics, Lund University,
Sweden  Agnieszka Pollo  Astronomical Observatory of Jagiellonian University,
Krakow  National Centre for Nuclear Research, Warsaw  Ioannis Liodakis
Finnish Centre for Astronomy with ESO (FINCA), University of Turku, Finland
Artem Poliszczuk  National Centre for Nuclear Research, Warsaw

(Received June 1, 2019; Revised January 10, 2019)

###### Abstract

Measuring the redshift of active galactic nuclei (AGNs) requires the use of
time-consuming and expensive spectroscopic analysis. However, obtaining
redshift measurements of AGNs is crucial as it can enable AGN population
studies, provide insight into the star formation rate, the luminosity
function, and the density rate evolution. Hence, there is a requirement for
alternative redshift measurement techniques.  In this project, we aim to use
the Fermi gamma-ray space telescope’s 4LAC Data Release 2 (DR2) catalog to
train a machine learning model capable of predicting the redshift reliably. In
addition, this project aims at improving and extending with the new 4LAC
Catalog the predictive capabilities of the machine learning (ML) methodology
published in Dainotti et al. (2021). Furthermore, we implement feature
engineering to expand the parameter space and a bias correction technique to
our final results. This study, uses additional machine learning techniques
inside the ensemble method, the SuperLearner, previously used in Dainotti et
al. (2021). Additionally, we also test a novel ML model called Sorted L-One
Penalized Estimation (SLOPE). Using these methods we provide a catalog of
estimated redshift values for those AGNs that do not have a spectroscopic
redshift measurement. These estimates can serve as a redshift reference for
the community to verify as updated Fermi catalogs are released with more
redshift measurements.

AGNs, Machine learning, redshift

††journal: APJ

##  1 Introduction

The Fermi $\gamma$-ray space telescope observes the sky in the energy range of
50 MeV to 1 TeV. With its $\gamma$-ray observational capabilities, it has
observed 3511 $\gamma$-ray loud active galactic nuclei (AGNs) to this day.
These AGNs are gathered in the latest Fourth Fermi-LAT Data Release 2 (4LAC)
catalog (Ajello et al., 2020; Abdollahi et al., 2020), along with their
measured $\gamma$-ray properties. The majority ($>98\%$) of AGNs in the 4LAC
catalog are blazars, which according to the Unified Model of AGN, are
identified with galaxies that have their relativistic emission jets pointed
directly along our line of sight (Blandford et al., 2019). Two major
subcategories of blazars are present in the 4LAC catalog, BL Lacertae Objects
(BLL) and Flat Spectrum Radio Quasars (FSRQ). BLLs are classified as sources
with weak or no emission lines in their spectra, whereas FSRQs are classified
as sources that have broad emission lines (Madejski & Sikora, 2016; Ghisellini
et al., 2011).

Measuring the redshift of blazars proves to be problematic as it often
requires multiwavelength observations and coordination across multiple
facilities. Further complications arise because the most dominant category of
$\gamma$-ray loud blazars, BLLs, often do not show detectable emission lines.
Not surprisingly, only 50% of $\gamma$-ray AGNs have redshift estimates,
highlighting the difficulty in measuring it, as well as the need for a
reliable photometric redshift estimation method. In light of these hurdles, a
predictive approach for redshift estimation using the more readily available
photometric $\gamma$-ray features is an attractive proposition.

Photometric redshifts can play a significant role in helping us understand the
origin of Extragalactic Background Light (EBL), crucial for blazars emitting
GeV-TeV $\gamma$-rays, which are typically BLLs (Wakely & Horan,
2008)111http://tevcat.uchicago.edu/. This can, in turn, enhance our knowledge
about the cosmic evolution of blazars (Singal et al., 2012, 2014; Singal,
2015; Singal et al., 2013a; Chiang et al., 1995; Ackermann et al., 2015;
Singal et al., 2013b; Ackermann et al., 2012), the structure of the magnetic
field in the intergalactic medium, (Marcotulli et al., 2020; Venters &
Pavlidou, 2013; Fermi-LAT Collaboration et al., 2018) as well as constraining
cosmological parameters (Domínguez et al., 2019; Petrosian, 1976; Singal et
al., 2013b).  In recent years the number of studies focusing on photometric
redshift estimation of high redshift AGN, using a machine learning
(MLapproach, has increased significantly (Jones & Singal, 2017; Cavuoti et
al., 2014; Fotopoulou & Paltani, 2018; Logan & Fotopoulou, 2020; Yang et al.,
2017; Zhang et al., 2019; Curran, 2020; Nakoneczny et al., 2020; Pasquet-Itam
& Pasquet, 2018; Jones & Singal, 2017). This is primarily due to the
availability of large data sets from all-sky surveys like the Sloan Digital
Sky Survey (SDSS) (Aihara et al., 2011) and Wide-field Infrared Survey
Explorer (WISE) (Brescia et al., 2019; Ilbert et al., 2008; Hildebrandt et
al., 2010; Brescia et al., 2013; Wright et al., 2010; D’Isanto & Polsterer,
2018). However, when it comes to $\gamma$-ray loud AGNs, ML methods have been
primarily used for AGN classification (Doert & Errando, 2014; Hassan et al.,
2013; Einecke, 2016; Lefaucheur & Pita, 2017; Kang et al., 2019; Chiaro et
al., 2016; Liodakis & Blinov, 2019).
Dainotti et al. (2021) remains the only work in the field, as of this writing,
that investigates the capabilities of Fermi in estimating photometric
redshifts. However, in that study, a significant hindrance was the small
dataset, which also contained a large number of missing values. Since AGNs
need to have both observed spectroscopic redshift and complete observational
data to be included in the ML models, this led to only 50% of the data set
being eligible. As a result, the small training sample contained 730 AGNs.
Hence, in this research one of our primary focuses was on increasing the
training sample size, which was achieved by including AGNs from the 4LAC
catalog that are present below the galactic plane. A detailed description of
the differences between the data set of the current work and Dainotti et al.
(2021) is provided in Sec. 6.
We also increase the number of predictors using feature engineering techniques
(see Sec. 3.6). We also use additional ML models in our ensemble (see Sec.
3.3), chosen specifically for this updated data set, and an independent method
called Sorted L-One Penalized Estimation (SLOPE) (see Sec. 3.5). Furthermore,
we investigate the effects of bias correction methods (see Sec. 3.7) on our
results. Finally, we present a list of photometric redshift predictions for
the AGNs in the 4LAC that do not have a spectroscopic redshift measurement
(see Sec. 7).
With these updates, we were able to maintain low numbers of catastrophic
outlier percentage (5%), small values of the root mean square error (RMSE),
and the normalized median absolute deviation ($\sigma_{NMAD}$) as in Dainotti
et al. (2021). We also managed to increase the Pearson correlation value, r,
between our predicted and observed redshifts. The results on the validation
set are also improved by 10%. We here stress that this has been achieved on
training and validation data samples 52% larger than Dainotti et al. (2021)
(see Fig. 1).

##  2 The sample

Since 2008 the instruments onboard Fermi (LAT, GBM) have been continuously
detecting $\gamma$-ray sources (e.g., Gamma-Ray Bursts, AGNs, among others).
Out of the 3511 AGNs observed in the latest 4LAC catalog, only 1767 have
measured redshifts, 50.3% of the total data set. In the 4LAC catalog, two
categories of $\gamma$-ray loud AGNs dominate, BLLs and FSRQs. The catalog
contains a total of 1190 BLLs and 733 FSRQs. A significant percentage (43%) of
the AGNs are unclassified and are denoted as Blazar Candidate of Unknown Type
(BCU). We also have 9 Narrow Line Seyferts (NLSY1), 44 radio galaxies, and 19
AGNs of other categories. Keeping this analysis consistent with the work
previously done in Dainotti et al. (2021), and to improve the readability and
interpretability of the method, we choose to omit all the non-BLL and non-FSRQ
type AGNs from our data.

Figure 1: The histogram showing the overlapped distributions of z in Dainotti
et al. (2021) training data set (white) and the current training data set
(red).  Figure 2: The overlapped histogram between Dainotti et al. (2021) data
set and the current one, in the $1/(z+1)$ scale.

Furthermore, Fermi measures 13 properties of the 4LAC AGNs. These observations
are used as predictors for the redshift in our ML models. The corresponding
Gaia g-band magnitudes (Jordi et al., 2010) of the AGNs are also included. Out
of these 13 features or predictors, we consider 11, excluding
$Fractional\\_Variability$ and $Highest\\_Energy$. Both of them were excluded
due to a lack of observational data in 43% of the AGNs. The definitions of the
11 predictors are provided below.

  * •

Flux \- The integral photon flux, in photons/cm2/s, from 1 to 100 GeV.

  * •

Energy_Flux \- The energy flux, the units are in erg cm-2 s-1, in the 100 MeV
- 100 GeV range obtained by the spectral fitting in this range.

  * •

Significance \- The source detection significance in Gaussian sigma units, on
the range from 50 MeV to 1 TeV.

  * •

Variability_Index \- Summation of the log(likelihood) differences between the
fitted flux in every interval of time and the mean flux over the 50 MeV to 1
TeV range.

  * •

$\nu$ \- The synchrotron peak frequency, measured in Hz, in the observer
frame.

  * •

$\nu$f$\nu$ \- Spectral energy distribution at the synchrotron peak frequency
in the $\nu$fν spectrum.

  * •

PL_Index \- It is the photon index when fitting the spectrum with a power law,
in the energy range from 50 MeV to 1 TeV.

  * •

Pivot_Energy \- The energy, in MeV, at which the error in the differential
photon flux is minimal, is derived from the likelihood analysis in the range
from 100 MeV - 1 TeV.

  * •

LP_Index \- Photon index at pivot energy ($\alpha$) when fitting the spectrum
(100 MeV to 1 TeV) with Log Parabola.

  * •

LP_$\beta$ \- the spectral parameter ($\beta$) when fitting with Log Parabola
spectrum from 50 MeV to 1 TeV.

  * •

Gaia_G_Magnitude \- Gaia Magnitude at the g-band provided by the 4LAC, taken
from the Gaia Survey.

Figure 3: The full scatter matrix plot of all the variables defined above.
Here InvRedshift denotes $1/(z+1)$ scaled data. The generalization set refers
to samples which do not have a redshift observation.

###  2.1 Training, generalization and validation sets

Before we can train the ML models, we need to transform our predictors and
create the training and generalization sets. For the predictors’
$Flux,Energy\\_Flux,Significance,Variability\\_Index,\nu,\nu$f$\nu$, and
$Pivot\\_Energy$, a transformation to their logarithmic scale is necessary as
they span over several orders of magnitude. Also, note that we do not train
our ML models to predict the redshift directly, instead, it is trained to
predict $1/(z+1)$, where z is the redshift. With such a transformation, we
predict the scale factor $a$, which is equal to $1/(z+1)$; furthermore, we
found this transformation to be better for the overall performance of our ML
models. The 1/(z+1) is indeed a much more physical quantity to estimate, as it
is directly used for cosmological and redshift evolution studies in both AGNs
and GRBs. (Dainotti et al., 2013, 2017, 2021; Singal, 2015; Dainotti et al.,
2013a, b, 2021; Dainotti & Del Vecchio, 2017; Dainotti & Amati, 2018).
Furthermore, we also make cuts in the sample to remove data points with an
$LP\\_\beta$ $<$ 0.7, $LP\\_Index$ $>$ 1, and $LogFlux$ $>$ -10.5, as they are
outliers of the distribution. These cuts leave us with 1897 AGNs, out of which
1444 have observed redshift and the remaining 453 do not. The scatter matrix
plot in Fig. 3 shows the distribution of these 1897 samples.

The training set or the training sample is the data on which our ML models are
trained to learn the underlying relations between the response variable and
the predictors. Thus, it consists of AGNs with observed redshift (the response
variable) and predictors with complete data. Due to missing values in 3
predictors, namely $Log\nu$, $Log\nu f\nu$, and $Gaia\\_G\\_Magnitude$, we had
to remove an additional 332 AGNs. This leads to a training sample of 1112 data
points, which is 52% larger than Dainotti et al. (2021). The overlapped
histogram of the redshifts of the entire data and the training data is
presented in Fig. 4.

Figure 4: Overlapped histogram showing the training set and the total 4LAC
redshift distributions.

On the other hand, the generalization set contains AGNs without observed
redshift but with complete observational data. These are the AGNs for which we
want to estimate the photometric redshift. Similar to the training set, if we
remove the data points with missing values then we are left with 320 AGNs in
the generalization sample. These 320 AGNs are made up of 319 BLLs and 1 FSRQ.
The validation set consists of AGNs, which contain measured redshift and
complete predictor data similar to the training set. Hence, we take a cut of
10% from the training set to create the validation set. This set of 111 AGNs
will not be used in the model’s training and will serve as an independent test
to verify the cross-validated results.

##  3 Methodology

For this project, we expand upon the methodology published in Dainotti et al.
(2021), which includes a more accurate version of feature selection where we
select the best predictors for estimating $1/(z+1)$, additional ML models,
feature engineering, and bias corrections. All the techniques and methods are
described in the following subsections.

The metrics we use for measuring the goodness of our results are the
following:

  * •

Bias: Mean of the difference between the observed and predicted values.

  * •

$\sigma_{NMAD}$: Normalized median absolute deviation between the predicted
and observed measurements.

  * •

Pearson correlation coefficient, $r$, between the predicted and observed
measurements.

  * •

Root Mean Square Error (RMSE).

  * •

Standard Deviation $\sigma$.

We quote the measured values of these parameters for $\Delta z_{norm}$ and
$\Delta z$. Where:

$\Delta z=z_{observed}-z_{predicted}$ $\Delta z_{norm}=\frac{\Delta
z}{1+z_{observed}}$

Along with these we also quote the catastrophic outlier percentage, which is
defined as the percentage of predictions that lies beyond the 2$\sigma$ error.
The quoted metrics used here are the same as in Dainotti et al. (2021) and
other works in the field, thus allowing for easy comparison.

###  3.1 Cross-Validation

Cross-validation (CV) is a frequently used method (Browne, 2000; Refaeilzadeh
et al., 2009) which estimates the expected prediction error of ML models,
which in our case is the RMSE. Commonly, a $10$-fold cross-validation (10fCV)
is used. In 10fCV, ten equal partitions or folds are created at random. At
each step, nine folds are used as training data to fit a model and the 10th
fold is used as the test set, on which the prediction error is evaluated. This
process is iterated such that each fold is used exactly once as a test set.
The averaged errors of each test set provide a metric to judge the performance
of the ML model on novel data samples. Since the folds are generated randomly,
there is a risk of over-fitting to a specific data partition. Therefore, we
perform 10fCV 100 times, and for each iteration, the folds are generated
randomly. The results are an average of these 100 iterations. Such a step is
necessary to mitigate the uncertainty of 10fCV and derandomize our results.
However, it should also be noted that this step is computationally expensive,
as we need to train multiple ML algorithms 10$\times$100 = 1000 times.

###  3.2 Feature selection : LASSO inside 10fCV

One of the most crucial steps in any ML project is feature selection. It is
the process where we choose a subset of all the predictors present in the data
set that are better at predicting the response variable. Feature selection is
crucial because it makes the model easier to interpret and reduces the size of
the high-dimensional data, which enables algorithms to work faster. As a
consequence, this also reduces the risk of over-fitting. Here we use the Least
Absolute Shrinkage and Selection Operator (LASSO) method (Tibshirani, 1996).
The LASSO algorithm works as a shrinkage method by constraining the
$\ell^{1}$-norm of the solution vector to be less than or equal to a positive
number known as the tuning parameter ($\lambda$). A vector’s $\ell^{1}$-norm
is defined as the sum of the magnitudes of the components of that vector. Such
a constraint on the solution vector penalizes the predictor coefficients,
resulting in some of the coefficients being set to 0, effectively eliminating
those predictors. The $\lambda$ parameter controls the shrinkage that is
applied to the coefficients of the predictors.  Our analysis, uses the R
implementation of the GLMNET function and the LASSO feature selection mode by
setting the parameter $alpha$ to 1 (Hastie et al., 2017; Tibshirani et al.,
2012). The function varies $\lambda$ while fitting a linear model to the data
and assigns coefficients to the predictors to reduce the error. As mentioned
previously, the actual features that are chosen depending on the value of
$\lambda$. Thus, we select those predictors with a non-zero coefficient for
the highest $\lambda$ value such that the prediction errors are within 1
standard deviation (Friedman et al., 2010a; Birnbaum, 1962; Hastie &
Tibshirani, 1987, 1990; Friedman et al., 2010b).

In this project, we perform LASSO feature selection during the cross-
validation step. This means that instead of choosing a feature set for the
entire training data, we pick optimal features from the data of the nine folds
that will be used to train the model at that stage. Then the model predicts on
the test set (10th fold) using these features only, such that each of the ten
folds acts as a test set exactly once. Note that the LASSO feature selection
technique uses the response variable to choose the optimal subset of
predictors. Hence, this method assures us that we are not using information
from the test set of the 10fCV, and allows us to substantiate the claim that
the CV results accurately estimate the expected prediction error for this
method.

###  3.3 The ML algorithms used in our analysis

The choice of the best ML algorithms depends significantly on the problem
being addressed and the data set being used. That holds in our study as well.
Previously in Dainotti et al. (2021), we tested multiple ML algorithms in an
ensemble and selected the ones with the highest coefficients. We perform the
same process once again, but with an expanded data set and different
predictors. Unsurprisingly, different ML algorithms were chosen to be best
suited for this analysis compared to Dainotti et al. (2021). In Fig. 5, we
present the coefficients assigned to the various ML models we tested. Those
with a higher coefficient are better at predicting the redshift. Therefore, we
define a threshold value of 0.5 and choose all the models with a coefficient
greater than this. Consequently, we end up with six algorithms in our final
ensemble. These are Enhanced Adaptive Regression Through Hinge (EARTH), KSVM,
Cforest, Ranger, Random Forest, and Linear Model. We describe these methods
below.

Figure 5:  Left Panel: The coefficients assigned by SuperLearner to the
algorithms tested. We select the algorithms that have a coefficient above 0.5
to be incorporated into our ensemble. Right Panel: The RMSE error (risk) of
each of the algorithms, scaled to show the minimum risk algorithm at 0, which
is cforest.

####  3.3.1 EARTH

Enhanced Adaptive Regression Through Hinges (EARTH) is an implementation of
the multivariate adaptive regression splines method (MARS), based on the work
by Friedman & Roosen (1995). MARS is a non-parametric regression technique
that can better model the interactions of predictors and non-linearities than
a simple linear model. MARS starts with a constant term then fits a weighted
sum of functions where each function is either another constant term, a hinge
function, or a product of hinge functions. Each term added to the model
results in the maximal reduction of sum-of-square residual error.

####  3.3.2 KSVM

In this study, we use the support vector regression method in the ksvm
algorithm from the SuperLearner package. Support vector regression (SVR) is
the continuous predictive counterpart to support vector machine classification
(Cortes & Vapnik, 1995). SVR uses a kernel function that maps an input vector
onto a space of higher-dimension. Then it fits a hyperplane to the data, such
that the predictions are within a certain error margin. In this study, the
default parameters of the ksvm model in the SuperLearner package are used.
Specifically, the kernel used is the radial basis kernel, i.e. the Gaussian
kernel.

####  3.3.3 Cforest, Ranger, and Random Forest

The Cforest, Ranger, and Random Forest algorithms belong to the tree class of
supervised machine learning models. A single decision tree algorithm may
perform very well on the training data, but it is prone to overfitting.
Besides the true relationships, such models learn the noise present in the
training set and are referred to as models with high variance. To obtain an
optimal bias-variance trade-off, one can often use ensembles made of multiple
decision trees. In such models, several stages of randomness are introduced to
reduce the variance of the final prediction.
Random Forest (Breiman, 2001) works by generating multiple decision trees.
While building the tree the data is partitioned based on the values of a
randomly chosen subset of features in the training data (Ho, 1995). This
partitioning repeats until the pre-specified depth of the tree is reached. At
the nodes, the values of the response variables are averaged over and that
average is assigned as the prediction for an AGN that falls into that node.
The algorithm finally outputs the prediction as an average over all the
independent trees.

The Ranger algorithm is nearly equivalent to Random Forest with the difference
being a slightly quicker implementation and the ability to use extremely
randomized trees (ERT). ERTs are built in the same top-down procedure as
decision trees in Random Forest with the difference being that the cuts are
randomly generated, assuming a uniform distribution for each of the features
in the training set. The cut that reduces prediction error the most is chosen
(Geurts et al., 2006). ERTs have the advantage of reducing variance further
than Random Forest due to the increased randomness in the node-splitting
procedure.
The Cforest algorithm (Hothorn et al., 2006) is similar to Random Forest and
Ranger except that it uses conditional inference trees (Hothorn et al., 2015)
which splits based on significance tests instead of information gain. These
tree-based regression methods can fit very complex, sophisticated models.

####  3.3.4 Linear Model

The linear model is the simplest regression model. Given a set of features
$X_{1},..,X_{n}$, the linear model fits a function of the form
$f(X_{1},..,X_{n})=\beta_{0}+\sum_{i=1}^{n}\beta_{i}*X_{i}$, where the
coefficients $\beta_{i}$ are real numbers. This project uses the default
ordinary least squares (OLS) linear model from the SuperLearner package. Here,
the optimal $\beta$ parameter vector is picked such that it minimizes the
average squared error.

###  3.4 SuperLearner

SuperLearner (Van der Laan et al., 2007) is an ensemble learning algorithm
that utilizes 10fCV. Ensemble learning is an ML technique that combines
multiple individual algorithms into a single one. There are multiple methods
to do this (Opitz & Maclin, 1999). The advantage of such a method is that the
shortcomings of one constituent algorithm can be mitigated by leveraging other
algorithms. This overall has a positive impact, thereby boosting performance.
In the case of SuperLearner, this ensemble is created by assigning
coefficients to its constituent models. The coefficients are determined based
on the RMSE of each algorithm in 10fCV. Algorithms with lower RMSE are
assigned higher coefficients. The sum of all the coefficients always equals
one, and they are greater than or equal to 0. This ensemble created by
SuperLearner intrinsically has the lowest RMSE among all its constituent
models (Polley & Van der Laan, 2010).

###  3.5 SLOPE

Sorted L-One Penalized Estimation, or SLOPE (Bogdan et al., 2015), is a
supervised learning method for both prediction and feature selection. Given a
vector $b$ $\in$ $\mathbb{R}^{p}$, and a sequence of non-negative scalars
$\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}\geq...\geq\lambda_{p}\geq 0$ where
$\lambda\neq 0$, the Sorted L1-Norm is defined as
$\lambda_{1}b_{1}+\lambda_{2}b_{2}+...+\lambda_{p}b_{p}$. SLOPE is a natural
extension of the LASSO method. It aims to minimize a cost function with the
constraint that the Sorted L-1 Norm is less than or equal to a positive
constant (Opitz & Maclin, 1999). It does this by assigning larger weights for
larger parameter values. We test this independent method that is not included
in the SuperLearner ensemble.

###  3.6 Feature Engineering

It is common in machine learning to create new predictors from the pre-
existing features in the data sample (Ng, 2013). These techniques, called
feature engineering, can help boost the performance of the final model. The
importance of feature engineering is clear when considering that the response
variable may depend on combinations of the initial features that certain
machine learning algorithms cannot learn. These include, but are not limited
to, features such as ratios for neural networks, and count variables for tree
algorithms (Heaton, 2016). Continuing the work from Dainotti et al. (2021)
where we considered models trained with the features present in 4LAC, here we
extended the set of predictors.  Our engineered features are simple cross-
multiplications of the existing features, as described in Sec. 2 to
distinguish the feature sets, we denote the original features obtained from
the 4LAC catalog as Order 1 (O1). The cross multiplied features are called
Order 2 (O2). Additionally, the feature selection with LASSO and the
algorithms used is done in the same way for O1 and O2 cases. The results are
presented for both.

###  3.7 Bias Correction

We define bias as the difference between the predicted and the observed values
of a response variable. Training models on an unbalanced sample often leads to
predictions that have a high bias, although it should be noted that an
unbalanced sample is not the sole reason for predictions to be biased. This
outcome can occur for other reasons such as the ML methods used, and the
transformations applied to the predictors and the response variable. In the
case of an unbalanced sample, the model is better trained on the redshift
range that contains more observations. As a result, the model predicts more
frequently in those particular redshift ranges and may provide erroneous
results.
We apply the technique of Optimal Transport for our bias correction. In this
technique, the predicted and observed 1/(z+1) values are sorted in ascending
order, and then a linear fit is obtained between them. The bias is corrected
using the slope and intercept of this linear fit. The analysis is performed
for results which are obtained in the case of O1 predictors with SuperLearner
and is applied for BLLs and FSRQs independently. We correct the respective
predictions using the following formula:

$C_{preds}=B+A*U_{phot}$

where $C_{preds}$ are the corrected 1/(z+1) predictions, $U_{phot}$ is the
Superlearner 1/(z+1) predictions, $A$ and $B$ are the slope and intercept of
the linear fits, respectively. $C_{pred}$ values are then converted to the
z-scale to obtain the bias corrected redshift estimates.
The effect of these corrections is presented in Sec. 5.4. Following this
procedure, we also correct the predictions on the generalization set and we
present these results in Table 3, along with the predictions without the bias
correction.

##  4 Comparison

As mentioned, this study aims at extending the method developed by Dainotti et
al. (2021). Here we provide a comparison between the two studies.
Previously, in Dainotti et al. (2021), some of us worked with $2863$ sources
from the 4LAC, comprising of $658$ FSRQs, $1067$ BLLs, $1074$ BCUs, and $64$
sources categorized as radio galaxies, Narrow line Seyferts, or other non-
blazar AGNs. Out of those $2863$ sources, only BLLs and FSRQs, with no missing
data, were used, resulting in $657$ AGNs for the training set and $73$ AGNs
for the validation set. In the current study, the 4LAC catalog has been
significantly expanded from $2863$ total AGNs to $3511$, an increment of
$52\%$. This increment has been achieved by including the lower galactic
latitude $\gamma$-ray loud AGNs. Out of the total number of AGNs under
consideration, $50.3\%$ have a measured redshift. Thus, the first difference
between the current and Dainotti et al. (2021) is the enlarged sample. Indeed,
in Sec. 2.1, we describe how we obtained a training set of $1112$ AGNs, a
$52\%$ increase compared to Dainotti et al. (2021).
In addition, part of our extension of the work in Dainotti et al. (2021) is
the consideration of new ML algorithms. They had obtained the highest weights
for the algorithms Extreme Gradient Boosting (XGBoost), Random Forest, Big
LASSO, and Bayes Generalized Linear Model (Bayes GLM) using SuperLearner. For
this study, we obtain the highest weights for KSVM, Random Forest, Linear
Model, Ranger, Cforest, and EARTH (see Sec. 3.3) within the same SuperLearner
framework.
Previously, Dainotti et al. (2021) used the property $Highest\\_Energy$ as one
of the predictors and did not take into account for $Log\nu f\nu$, as this was
already considered a cross-product term. In this study, we are including
$Log\nu f\nu$ and excluding $Highest\\_Energy$ (see Sec. 2 for additional
details), along with introducing the use of O2 predictors (see Sec. 3.6).
In Dainotti et al. (2021), LASSO feature selection was performed on the entire
training set. Here, we perform LASSO feature selection on a fold-by-fold basis
inside the 10fCV leading to more accurate results (see Sec. 3.2). This
procedure was applied in all the trails that use SuperLearner, regardless of
the predictor set.

##  5 Results

In this study, we tested multiple models with different parameter spaces to
obtain the best possible model for predicting the redshift of $\gamma$-ray
loud AGNs. For each of the methods, we perform 10fCV 100 times. The method
which leads to the most promising results is the SuperLearner ensemble using
O1 predictors. We then contrast this with two other experiments: SuperLearner
ensemble with O2 predictors, and the SLOPE method with O1 predictors.
The experiment with SuperLearner using O2 predictors is intended to test how
feature engineering affects the predictive capabilities of the model. For this
we have 78 O2 predictors compared to 12 O1 predictors.
With SLOPE, we want to test how this single algorithm performs with respect to
the SuperLearner ensemble. We are using the O1 predictor set with SLOPE and
the same feature selection method using LASSO.

###  5.1 SuperLearner ensemble with O1 predictors

Using the O1 predictor set with the SuperLearner ensemble presents us the best
results amongst the experiments we tested.

Figure 6: 10fCV results. The orange filled circles represent the BLLs and the
blue filled triangles represent the FSRQs. The empty circles and triangles
represent the BLLs and FSRQs that are catastrophic outliers, respectively. Top
left panel: cross-validation correlation plot in the 1/(z+1) scale. Top right
panel: cross-validation correlation plot in the linear scale. Bottom left
panel: Validation set correlation plot in the 1/(z+1) scale. Bottom right
panel: Validation set correlation plot in the linear scale.

The results obtained are shown in the Fig. 6. The top left and right panels
show the cross-validated results, in the 1/(z+1) and linear scales,
respectively. In the 1/(z+1) scale we obtain a Pearson correlation of $78.1\%$
between the predicted and observed values, an RMSE of $0.119$, $\sigma$ of
$0.119$, and $\sigma_{NMAD}$ of $0.113$. In the linear scale, the correlation
between predicted and observed redshift is $74\%$, RMSE is $0.467$, bias is
$0.094$, and $\sigma_{NMAD}$ is $0.321$. We obtain a low catastrophic outlier
percentage of $5\%$: $1052$ AGNs are predicted within the $2\sigma$ range.

The lines in blue depict the 2$\sigma$ curves for each plot, where the
$\sigma$ is calculated in the 1/(z+1) scale as follows:

$\frac{1}{z_{p}+1}=\frac{1}{z_{s}+1}\pm 2\sigma,$

where $z_{s}$ is the spectroscopic redshift and $z_{p}$ is the photometric
redshift. The choice of our scaling causes the 2$\sigma$ curve to not be a
straight line on the linear scale. Thus, we use the formula mentioned below to
draw them in the linear-scale.

$z_{p}=z_{s}\left[\frac{1\pm 2\sigma(z_{p}+1)}{1\mp
2\sigma}\right]\pm\frac{2\sigma}{1\mp 2\sigma.}$

Figure 7:  Top Left Panel: Distribution of $\Delta z$ with average bias (red)
and sigma lines (blue). Top Right Panel: Distribution of the $\Delta z_{norm}$
with the average bias (red) and sigma values (blue). Bottom Left Panel:
Relative influence of the O1 predictors. Bottom right panel: Distribution of
the correlations in the linear scale from the 100 iterations of 10fCV.

In the Fig. 7 the top left and right panels show the histogram distribution of
the $\Delta$z and $\Delta$znorm. The bias is denoted by the red line and the
blue lines show the $\sigma$ values. The bottom left panel shows the relative
influence of our predictors with $LP\\_\beta$ with the highest relative
influence, followed by the $Log\nu$ and $Log\nu f\nu$. The categorical
predictor $LabelNo$, which distinguishes between BLL and FSRQ, is the 7th most
influential predictor. These relative influence values presented here are an
average over 100 iterations. Such a step is necessary as it serves to
derandomize and stabilize the relative influence results, and thus they can be
compared with Dainotti et al. (2021). Finally, the bottom right panel shows
the distribution of the correlation value between observed and predicted
redshift across the 100 iterations of SuperLearner. As shown, these
correlations lie completely within 0.73 to 0.745, with the majority falling
between 0.735 and 0.74.

Figure 8: These distribution plots are derived from the 100 iterations of
10fCV. They represent how dispersed these metrics are for the data and
algorithms at hand. Top left panel: Distribution of $\sigma_{NMAD}$ in the
1/(z+1) scale. The $\sigma_{NMAD}$ is written as NMAD on the plots. Top right
panel: Distribution of the $\sigma_{NMAD}$ in the linear scale. Bottom left
panel: Distribution of the RMSE in the 1/(z+1) scale. Bottom right panel:
Distribution of the RMSE in the linear scale.

In Fig. 8 we are presenting the distributions of two parameters, namely RMSE
and $\sigma_{NMAD}$. The top left and right panels show the distribution of
$\sigma_{NMAD}$ in the 1/(z+1) and the linear scales, respectively. Similarly,
the bottom left and right panels show the distribution of RMSE in 1/(z+1) and
the linear scales, respectively.

###  5.2 Comparison with the other methods

In this section, we present a comparative analysis between the O1 predictors
with SuperLearner and the two experiments we mentioned previously.

When using the O2 predictor set with the SuperLearner ensemble we obtain
comparable results to the O1 case. In the 1/(z+1) scale, the correlation is
0.782, which is within 1%. Similarly, we obtain a RMSE of 0.119, bias of
4.5$\times 10^{-4}$, and $\sigma_{NMAD}$ of 0.114.
In the linear scale, the correlation is 0.74, which is equal to the value
obtained with the O1 predictors. The RMSE is 0.463, which is a 1% improvement.
The bias obtained with O2 is 0.088, which is similarly a 1% improvement. The
$\sigma_{NMAD}$ obtained is 0.328.

As mentioned previously, we also used the SLOPE method. Compared to the
SuperLearner ensemble with O1 predictors, it performs nearly as well despite
being a singular algorithm. Specifically, we obtain a Pearson correlation of
76% in the 1/(z+1) scale, a RMSE of 0.124, bias of 3$\times 10^{-5}$, and
$\sigma_{NMAD}$ of 0.12.
In the linear scale, the correlation is 72%, the RMSE is 0.479, bias is 0.1
and $\sigma_{NMAD}$ is 0.348. Compared to the O1 case, in the linear scale,
the R is down by 2.7%, RMSE is up by 2.5%, and the $\sigma_{NMAD}$ is up by
8.4%. In both the linear scale and the 1/(z+1) scale the O1 predictors with
SuperLearner outperforms SLOPE, although by only a small margin.

All the results are presented in Table. 1 and Table. 2, for the 1/(z+1) and
the linear scales, respectively.

Metric | SL with O1 | SL with O2 | SLOPE with O1
---|---|---|---
R | 0.781 | 0.782 | 0.762
RMSE | 0.119 | 0.119 | 0.124
Bias | $4\times 10^{-4}$ | $4.5\times 10^{-4}$ | $3\times 10^{-5}$
$\sigma_{NMAD}$ | 0.113 | 0.114 | 0.12
$\sigma$ | 0.119 | 0.119 | 0.124
Table 1: Metrics to compare the different experiments performed in the 1/(z+1) scale. Metric | SL with O1 | SL with O2 | SLOPE with O1
---|---|---|---
R | 0.74 | 0.74 | 0.72
RMSE ($\Delta z$) | 0.467 | 0.463 | 0.479
Bias ($\Delta z$) | 0.095 | 0.089 | 0.1
$\sigma_{NMAD}$ ($\Delta z$) | 0.321 | 0.328 | 0.348
$\sigma$ ($\Delta z$) | 0.458 | 0.455 | 0.47
Bias ($\Delta z_{norm}$) | $9.6\times 10^{-4}$ | $-5.1\times 10^{-4}$ | $6.9\times 10^{-4}$
$\sigma_{NMAD}$ ($\Delta z_{norm}$) | 0.195 | 0.195 | 0.205
$\sigma$ ($\Delta z_{norm}$) | 0.208 | 0.209 | 0.217
Table 2: Metrics to compare the different experiments performed in the linear
scale.

###  5.3 Prediction on the generalization set

As described in Sec. 2.1 the generalization set consists of AGNs that do not
have a measured redshift.

Figure 9: The scatter matrix plot for the BLLs in the training (red) and
generalization (blue) samples.

When predicting on this set it is crucial to ensure that the data points fall
within the parameter space of the trained model. Otherwise, the model will
extrapolate those predictions and our confidence in them will be drastically
reduced. We ensure this does not happen by only keeping those AGNs in the
generalization set whose maximum and minimum values for each predictor fall
within the range of the same predictor in the training set. This is
demonstrated in Fig. 9, where the blue points (Generaliztion set BLLs) lie
within the range of the red points (Training set BLLs). This leaves us with
305 AGNs, all of which are BLLs.
In Fig. 10, we show the distribution of the predicted values of 305 BLLs,
overlapped with the redshift distribution of the BLLs in the training sample.
These predictions are made by using the O1 predictor set with the SuperLearner
ensemble.

Figure 10: The distribution of the predicted redshifts for the generalization
set from the model trained on O1 predictors, overlapped with the redshift
distribution of the BLLs in the training set.

####  5.3.1 Comparative results

Here we consider the intersection of the generalization sets of both the
current study and Dainotti et al. (2021), which consists of 240 AGNs. By
plotting the estimated redshifts of these AGNs from the two studies, we see
that there is a Pearson correlation of 87%. This correlation plot is presented
in Fig. 11. We also observe a low RMSE and $\sigma$ of 0.078 and 0.079,
respectively. The red line denotes the equality line, while the blue lines
denote the 2$\sigma$ error curves. Based on the position of the points in the
plot, we note here that the redshift estimate of only 8 AGNs (3%) lie beyond
the 2$\sigma$ error limit. It should be noted that both these sets of
estimates have not been bias corrected.

Figure 11: Comparision between the redshift estimates of 240 generalization
set AGNs, common between this study and Dainotti et al. (2021). We obtain a
correlation of 0.87, RMSE of 0.078 and $\sigma$ of 0.079.

###  5.4 Bias Correction for the Generalization set

Following the bias-correction described in Sec. 3.7, the modified 10fCV
results are presented in Fig. 13. The left and right panels of Fig. 12 present
the optimal transport linear fit of BLLs and FSRQs, respectively.

Figure 12: Optimal transport linear fits of BLLs, left panel, and FSRQs, right
panel.

We note that there is a 8.3% reduction in the correlation in the linear scale,
15.6% increase in the RMSE, and a 3.4% decrease in the $\sigma_{NMAD}$. We
also see an increase in catastropic outliers from 5% to 6%. However, in actual
numbers there are only two additional AGNs that are outside the 2$\sigma$
cones, as compared to the O1 scenario. These results are quite expected, since
with bias correction we are reducing the accuracy in the bulk of the redshift
distribution, in order to get better results for the tails of the
distribution. The application of this bias correction to the predictions on
the generalization set lead to the results shown in Fig. 14. We present the
full set of bias-corrected predictions in Table. 3.

Figure 13: The bias-corrected predictions of the 10fCV. The symbols here
retain the same meaning as those in Fig. 6. Left Panel: Correlation plot
between the bias-corrected prediction and observed redshift in the 1/(z+1)
scale. Right Panel: Correlation plot between the bias-corrected prediction and
observed redshift in the linear scale Figure 14: Histogram distributions of
bias-corrected predictions of the generalization set AGNs, overlapped with the
training set AGNs.

##  6 Conclusions & Discussion

In Dainotti et al. (2021), some of us had obtained a Pearson correlation of
0.71 and an RMSE of 0.434 in the linear scale. In this study, using the O1
predictors with the SuperLearner we achieve a Pearson correlation of 0.74 and
an RMSE of 0.467. Specifically, an improvement of 4.2% in the correlation and
an increase of 7.6% in RMSE. At the same time, the percentage of catastrophic
outliers remains the same at 5% between this study and Dainotti et al. (2021).
We found that using O2 predictors leads to performance equivalent to O1
predictors. Thus, we can deduce that O1 predictors contain all the information
needed to predict the redshift. This is evident in Fig. 15 as the top ten O2
predictors with the highest relative influence are either O1 predictors or
cross products of other high influence O1 predictors (see bottom left panel of
Fig. 7). Such an outcome could possibly be exclusive to the 4LAC properties,
and indeed cross-product terms in other data sets might lead to better
performance.
Additionally we tested SLOPE independently of SuperLearner. On average we
observe that the results from SLOPE are comparable to those from the ensemble
via the metrics presented in Table 1. However, the advantage of SLOPE is that
it works much faster than the SuperLearner ensemble, along with the added
benefit of high interpretability of the algorithm. Hence, when the requirement
is to obtain estimates quickly, SLOPE can be used as a reliable alternative.

Here we summarise about the predictors used in our study.

  1. 1.

$LP\\_\beta$, $Log\nu f\nu$, $Log\nu$, $LogPivot\\_Energy$,
$LogEnergy\\_Flux$, $PL\\_Index$, $LabelNo$, and $LogSignificance$ have
highest influence in this study.

  2. 2.

$LP\\_\beta$, $LogPivot\\_Energy$, $LogSignificance$, $LogEnergy\\_Flux$,
$Gaia\\_G\\_Magnitude$, $Log\nu$, and $Log\\_Highest\\_Energy$ had the highest
influence in Dainotti et al. (2021).

  3. 3.

In this study and Dainotti et al. (2021), $LP\\_\beta$ has the highest
influence at 25%.

  4. 4.

With O2 predictors, $LP\\_\beta$ has a relative influence of 3% at the 10th
position. Three O2 predictors more influential than $LP\\_\beta$ are cross-
products of it. These are $LP\\_\beta^{2}$, $LogSignificance*LP\\_\beta$ and
$LogVariability\\_Index*LP\\_\beta$, with $LogVariability\\_Index*LP\\_\beta$
being the second most influential at 10%.

  5. 5.

Hence, we can conclude that for redshifts estimation of 4LAC AGNs,
$LP\\_\beta$ is crucial.

  6. 6.

$Gaia\\_G\\_Magnitude$ has a relative influence of 1% at 12th position. In
Dainotti et al. (2021) however it was the 5th most influential at 12% relative
influence. This shift in the relative influence of $Gaia\\_G\\_Magnitude$ can
be attributed to the enlarged data set. However, this is not conclusive and
merits further investigation.

  7. 7.

$PL\\_Index$ has the highest influence when paired with O2 predictors at 15%.
With O1 predictors, it was the 6th most influential at 11% influence.

Finally, to make the 4LAC more complete, we present two lists containing the
predictions on the generalization set. One contains the predictions from the
SuperLearner ensemble with O1 predictors (shown in Fig. 10) and the second
contains the bias-corrected results (shown in Fig. 14). Additionally, we
provide a comparison between the redshift estimates of 240 AGNs that are
common between the generalization sets of this study and Dainotti et al.
(2021) (see Fig. 11). These two results have a Pearson correlation of 0.87,
RMSE and $\sigma$ of 0.078 and 0.079, respectively.
This is the most complete investigation into the predictive capabilities of
the 4LAC dealing with redshift estimation. Compared to the Dainotti et al.
(2021), we have been able to maintain the predictive capabilities of our
methodology with an even larger data sample. These estimates will likely serve
as a reference point for the community to be used as tools for cosmological
investigations. The final aim is to use the estimation of the redshift to
investigate AGN with interesting properties.

##  7 Appendix

Here we are presenting Table. 3 containing the predictions of fifty
generalization set AGNs. It contains the names of AGNs, as assigned in the
4LAC catalog, followed by the predicted redshift from the SuperLearner
ensemble with O1, and finally the bias-corrected redshifts. The predictions
for the complete generalization set can be found in the data file provided.

Name | Predicted $z_{pred}$ | Bias-corrected $z_{pred}$
---|---|---
4FGL J0009.3+5030 | 0.57 | 0.62
4FGL J0001.2-0747 | 0.39 | 0.39
4FGL J0013.1-3955 | 0.7 | 0.90
4FGL J0019.3-8152 | 0.24 | 0.17
4FGL J0019.6+2022 | 0.79 | 1.09
4FGL J0021.5-2552 | 0.48 | 0.51
4FGL J0022.1-1854 | 0.57 | 0.65
4FGL J0022.5+0608 | 0.76 | 1.03
4FGL J0023.9+1603 | 0.63 | 0.78
4FGL J0029.0-7044 | 0.59 | 0.72
4FGL J0031.3+0726 | 0.32 | 0.27
4FGL J0045.3+2128 | 0.34 | 0.28
4FGL J0056.4-2118 | 0.51 | 0.54
4FGL J0058.0-3233 | 0.54 | 0.62
4FGL J0058.3+1723 | 0.37 | 0.36
4FGL J0107.4+0334 | 0.4 | 0.48
4FGL J0112.8-7506 | 0.39 | 0.38
4FGL J0115.6+0356 | 0.48 | 0.52
4FGL J0116.0-2745 | 0.31 | 0.26
4FGL J0116.2-6153 | 0.39 | 0.37
4FGL J0120.4-2701 | 0.34 | 0.28
4FGL J0125.3-2548 | 0.71 | 0.95
4FGL J0127.2+0324 | 0.34 | 0.29
4FGL J0132.8-4413 | 0.37 | 0.35
4FGL J0136.5+3906 | 0.34 | 0.28
4FGL J0138.0+2247 | 0.26 | 0.19
4FGL J0139.0+2601 | 0.33 | 0.29
4FGL J0142.7-0543 | 0.36 | 0.32
4FGL J0143.7-5846 | 0.21 | 0.11
4FGL J0144.6+2705 | 0.73 | 0.95
4FGL J0153.0+7517 | 0.31 | 0.26
4FGL J0155.0+4433 | 0.54 | 0.62
4FGL J0158.5-3932 | 0.4 | 0.39
4FGL J0159.7-2740 | 0.39 | 0.37
4FGL J0202.7+4204 | 0.51 | 0.57
4FGL J0203.6+7233 | 0.47 | 0.51
4FGL J0208.3-6838 | 0.44 | 0.46
4FGL J0209.3+4449 | 0.29 | 0.23
4FGL J0219.5+0724 | 0.42 | 0.44
4FGL J0226.5-4441 | 0.32 | 0.28
4FGL J0232.9+2608 | 0.35 | 0.30
4FGL J0241.0-0505 | 0.58 | 0.69
4FGL J0243.4+7119 | 0.66 | 0.84
4FGL J0245.1-0257 | 0.45 | 0.48
4FGL J0258.1+2030 | 0.6 | 0.73
4FGL J0304.9-0606 | 0.48 | 0.50
4FGL J0309.7-0745 | 0.32 | 0.26
4FGL J0310.6-5017 | 0.28 | 0.21
4FGL J0312.5-2221 | 0.37 | 0.34
4FGL J0314.3-5103 | 0.55 | 0.64
4FGL J0316.2-6437 | 0.41 | 0.39
Table 3: List of predictions for the generalization set. The first column has
the names of the AGNs as they appear in the 4LAC catalog. Second column
contains the redshift predictions from the SuperLearner ensemble using O1
predictors. Third column contains the same predictions with bias correction.
Here $z_{pred}$ denotes the predicted redshift.

The relative influence of the top 10 O2 predictors are shown in Fig. 15

Figure 15: This plot presents the relative influence of top 10 the O2
predictors.

Due to BLLs being the major category of AGNs present in the generalization set
of the 4LAC catalog, we here present the distribution of predicted redshift
made from a Superlearner model trained exclusively on the BLLs of the training
set (501 BLL). We perform the Kolmogorov Smirnov (KS) test to compare the
distribution of the predicted redshifts of generalization samples of BLLs,
obtained from the model trained only on BLLs and the O1 SuperLearner model. We
conclude based on the KS test (p-value=0.66) that we cannot reject the null
hypothesis that the two distributions are drawn from the same population.
Thus, we have stable results, which are resistant to the construction of the
training sample, either with only BLLs or with the total 4LAC sample.

Figure 16: Distribution of the redshifts of the training data (red) and the
predicted redshifts of the generalization set, from a Superlearner trained
only on BLLs (turquoise).

This work presents results from the European Space Agency (ESA) space mission,
Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis
Consortium (DPAC). Funding for the DPAC is provided by national institutions,
in particular the institutions participating in the Gaia MultiLateral
Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia.
The Gaia archive website is https://archives.esac.esa.int/gaia. M.G.D. thanks
Trevor Hastie for the interesting discussion on overfitting problems. We also
thank Raymond Wayne for the initial computation and discussions about balanced
sampling techniques which will be implemented in subsequent papers. This
research was supported by the Polish National Science Centre grant
UMO-2018/30/M/ST9/00757 and by Polish Ministry of Science and Higher Education
grant DIR/WK/2018/12. We would also like to acknowledge Enrico Rinaldi and
Biagio De Simone for their invaluable discussions regarding bias correction
and SLOPE implementation.

## References

  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Ackermann et al. (2012) Ackermann, M., Ajello, M., Allafort, A., et al. 2012, The Astrophysical Journal, 755, 164
  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Albert, A., et al. 2015, The Astrophysical Journal Letters, 813, L41
  * Aihara et al. (2011) Aihara, H., Prieto, C. A., An, D., et al. 2011, The Astrophysical Journal Supplement Series, 193, 29
  * Ajello et al. (2020) Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105, doi: 10.3847/1538-4357/ab791e
  * Birnbaum (1962) Birnbaum, A. 1962, Journal of the American Statistical Association, 57, 269
  * Blandford et al. (2019) Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467, doi: 10.1146/annurev-astro-081817-051948
  * Bogdan et al. (2015) Bogdan, M., Van Den Berg, E., Sabatti, C., Su, W., & Candès, E. J. 2015, The annals of applied statistics, 9, 1103
  * Breiman (2001) Breiman, L. 2001, Machine learning, 45, 5
  * Brescia et al. (2013) Brescia, M., Cavuoti, S., D’Abrusco, R., Longo, G., & Mercurio, A. 2013, The Astrophysical Journal, 772, 140
  * Brescia et al. (2019) Brescia, M., Salvato, M., Cavuoti, S., et al. 2019, MNRAS, 489, 663, doi: 10.1093/mnras/stz2159
  * Browne (2000) Browne, M. W. 2000, Journal of mathematical psychology, 44, 108
  * Cavuoti et al. (2014) Cavuoti, S., Brescia, M., D’Abrusco, R., Longo, G., & Paolillo, M. 2014, Monthly Notices of the Royal Astronomical Society, 437, 968
  * Chiang et al. (1995) Chiang, J., Fichtel, C., Von Montigny, C., Nolan, P., & Petrosian, V. 1995, The Astrophysical Journal, 452, 156
  * Chiaro et al. (2016) Chiaro, G., Salvetti, D., La Mura, G., et al. 2016, MNRAS, 462, 3180, doi: 10.1093/mnras/stw1830
  * Cortes & Vapnik (1995) Cortes, C., & Vapnik, V. 1995, Machine learning, 20, 273
  * Curran (2020) Curran, S. 2020, Monthly Notices of the Royal Astronomical Society: Letters, 493, L70
  * Dainotti & Amati (2018) Dainotti, M. G., & Amati, L. 2018, PASP, 130, 051001, doi: 10.1088/1538-3873/aaa8d7
  * Dainotti et al. (2013a) Dainotti, M. G., Cardone, V. F., Piedipalumbo, E., & Capozziello, S. 2013a, MNRAS, 436, 82, doi: 10.1093/mnras/stt1516
  * Dainotti et al. (2021) Dainotti, M. G., De Simone, B., Schiavone, T., et al. 2021, The Astrophysical Journal, 912, 150
  * Dainotti et al. (2021) Dainotti, M. G., De Simone, B., Schiavone, T., et al. 2021, ApJ, 912, 150, doi: 10.3847/1538-4357/abeb73
  * Dainotti & Del Vecchio (2017) Dainotti, M. G., & Del Vecchio, R. 2017, New A Rev., 77, 23, doi: 10.1016/j.newar.2017.04.001
  * Dainotti et al. (2017) Dainotti, M. G., Hernandez, X., Postnikov, S., et al. 2017, The Astrophysical Journal, 848, 88
  * Dainotti et al. (2013b) Dainotti, M. G., Petrosian, V., Singal, J., & Ostrowski, M. 2013b, ApJ, 774, 157, doi: 10.1088/0004-637X/774/2/157
  * Dainotti et al. (2013) Dainotti, M. G., Singal, J., Ostrowski, M., et al. 2013, The Astrophysical Journal, 774, 157
  * Dainotti et al. (2021) Dainotti, M. G., Bogdan, M., Narendra, A., et al. 2021, Predicting the redshift of gamma-ray loud AGNs using supervised machine learning.  https://arxiv.org/abs/2107.10952
  * D’Isanto & Polsterer (2018) D’Isanto, A., & Polsterer, K. L. 2018, A&A, 609, A111, doi: 10.1051/0004-6361/201731326
  * Doert & Errando (2014) Doert, M., & Errando, M. 2014, The Astrophysical Journal, 782, 41
  * Domínguez et al. (2019) Domínguez, A., Wojtak, R., Finke, J., et al. 2019, ApJ, 885, 137, doi: 10.3847/1538-4357/ab4a0e
  * Einecke (2016) Einecke, S. 2016, Galaxies, 4, 14
  * Fermi-LAT Collaboration et al. (2018) Fermi-LAT Collaboration, Abdollahi, S., Ackermann, M., et al. 2018, Science, 362, 1031, doi: 10.1126/science.aat8123
  * Fotopoulou & Paltani (2018) Fotopoulou, S., & Paltani, S. 2018, Astronomy & Astrophysics, 619, A14
  * Friedman et al. (2010a) Friedman, J., Hastie, T., & Tibshirani, R. 2010a, Journal of Statistical Software, 33, 1.  https://www.jstatsoft.org/v33/i01/
  * Friedman et al. (2010b) —. 2010b, Journal of statistical software, 33, 1
  * Friedman & Roosen (1995) Friedman, J. H., & Roosen, C. B. 1995, Statistical Methods in Medical Research, 4, 197, doi: 10.1177/096228029500400303
  * Geurts et al. (2006) Geurts, P., Ernst, D., & L., W. 2006, Machine Learning, 63, 42, doi: 10.1007/s10994-006-6226-1
  * Ghisellini et al. (2011) Ghisellini, G., Tavecchio, F., Foschini, L., & Ghirlanda, G. 2011, Monthly Notices of the Royal Astronomical Society, 414, 2674
  * Hassan et al. (2013) Hassan, T., Mirabal, N., Contreras, J., & Oya, I. 2013, Monthly Notices of the Royal Astronomical Society, 428, 220
  * Hastie & Tibshirani (1987) Hastie, T., & Tibshirani, R. 1987, Journal of the American Statistical Association, 82, 371
  * Hastie et al. (2017) Hastie, T., Tibshirani, R., & Tibshirani, R. J. 2017, arXiv preprint arXiv:1707.08692
  * Hastie & Tibshirani (1990) Hastie, T. J., & Tibshirani, R. J. 1990, Generalized additive models, Vol. 43 (CRC press)
  * Heaton (2016) Heaton, J. 2016, SoutheastCon 2016, doi: 10.1109/secon.2016.7506650
  * Hildebrandt et al. (2010) Hildebrandt, H., Arnouts, S., Capak, P., et al. 2010, Astronomy & Astrophysics, 523, A31
  * Ho (1995) Ho, T. K. 1995, in Proceedings of the Third International Conference on Document Analysis and Recognition (Volume 1) - Volume 1, ICDAR ’95 (USA: IEEE Computer Society), 278
  * Hothorn et al. (2006) Hothorn, T., Hornik, K., & Zeileis, A. 2006, Journal of Computational and Graphical statistics, 15, 651
  * Hothorn et al. (2015) —. 2015, The comprehensive R archive network, 8
  * Ilbert et al. (2008) Ilbert, O., Capak, P., Salvato, M., et al. 2008, The Astrophysical Journal, 690, 1236
  * Jones & Singal (2017) Jones, E., & Singal, J. 2017, Astronomy & Astrophysics, 600, A113
  * Jordi et al. (2010) Jordi, C., Gebran, M., Carrasco, J. M., et al. 2010, A&A, 523, A48, doi: 10.1051/0004-6361/201015441
  * Kang et al. (2019) Kang, S.-J., Fan, J.-H., Mao, W., et al. 2019, The Astrophysical Journal, 872, 189
  * Lefaucheur & Pita (2017) Lefaucheur, J., & Pita, S. 2017, Astronomy & Astrophysics, 602, A86
  * Liodakis & Blinov (2019) Liodakis, I., & Blinov, D. 2019, MNRAS, 486, 3415, doi: 10.1093/mnras/stz1008
  * Logan & Fotopoulou (2020) Logan, C., & Fotopoulou, S. 2020, Astronomy & Astrophysics, 633, A154
  * Madejski & Sikora (2016) Madejski, G. G., & Sikora, M. 2016, ARA&A, 54, 725, doi: 10.1146/annurev-astro-081913-040044
  * Marcotulli et al. (2020) Marcotulli, L., Ajello, M., & Di Mauro, M. 2020, in American Astronomical Society Meeting Abstracts, Vol. 235, American Astronomical Society Meeting Abstracts #235, 405.06
  * Nakoneczny et al. (2020) Nakoneczny, S. J., Bilicki, M., Pollo, A., et al. 2020, Photometric selection and redshifts for quasars in the Kilo-Degree Survey Data Release 4.  https://arxiv.org/abs/2010.13857
  * Ng (2013) Ng, A. 2013, Accessed: May, 3, 2018
  * Opitz & Maclin (1999) Opitz, D., & Maclin, R. 1999, Journal of artificial intelligence research, 11, 169
  * Pasquet-Itam & Pasquet (2018) Pasquet-Itam, J., & Pasquet, J. 2018, Astronomy & Astrophysics, 611, A97
  * Petrosian (1976) Petrosian, V. 1976, The Astrophysical Journal, 209, L1
  * Polley & Van der Laan (2010) Polley, E. C., & Van der Laan, M. J. 2010
  * Refaeilzadeh et al. (2009) Refaeilzadeh, P., Tang, L., & Liu, H. 2009, Encyclopedia of database systems, 5, 532
  * Singal (2015) Singal, J. 2015, Monthly Notices of the Royal Astronomical Society, 454, 115
  * Singal et al. (2013a) Singal, J., Ko, A., & Petrosian, V. 2013a, Proceedings of the International Astronomical Union, 9, 149
  * Singal et al. (2014) —. 2014, The Astrophysical Journal, 786, 109
  * Singal et al. (2012) Singal, J., Petrosian, V., & Ajello, M. 2012, The Astrophysical Journal, 753, 45
  * Singal et al. (2013b) Singal, J., Petrosian, V., & Ko, A. 2013b, AAS/High Energy Astrophysics Division# 13, 300
  * Tibshirani (1996) Tibshirani, R. 1996, Journal of the Royal Statistical Society Series B, 58, 267
  * Tibshirani et al. (2012) Tibshirani, R., Bien, J., Friedman, J., et al. 2012, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74, 245
  * Van der Laan et al. (2007) Van der Laan, M. J., Polley, E. C., & Hubbard, A. E. 2007, Statistical applications in genetics and molecular biology, 6
  * Venters & Pavlidou (2013) Venters, T. M., & Pavlidou, V. 2013, MNRAS, 432, 3485, doi: 10.1093/mnras/stt697
  * Wakely & Horan (2008) Wakely, S. P., & Horan, D. 2008, in International Cosmic Ray Conference, Vol. 3, International Cosmic Ray Conference, 1341–1344
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, The Astronomical Journal, 140, 1868–1881, doi: 10.1088/0004-6256/140/6/1868
  * Yang et al. (2017) Yang, Q., Wu, X.-B., Fan, X., et al. 2017, The Astronomical Journal, 154, 269
  * Zhang et al. (2019) Zhang, K., Schlegel, D. J., Andrews, B. H., et al. 2019, The Astrophysical Journal, 883, 63