# 3D scattering microphantom sample to assess quantitative accuracy in
tomographic phase microscopy techniques

Wojciech Krauze,1,∗ Arkadiusz Kuś,1 Michał Ziemczonok,1 Max Haimowitz,2
Shwetadwip Chowdhury2 Małgorzata Kujawińska1

(1Warsaw University of Technology, Institute of Micromechanics and Photonics,
Boboli 8 street, Warsaw, Poland, 02-525
2University of Texas at Austin, Department of Electrical Engineering, 2501
Speedway, Austin, TX 78712, USA)

###### Abstract

In this paper we present a method to robustly evaluate the quantitative
accuracy of various tomographic phase microscopy (TPM) methods with a multiple
scattering 3D-printed microphantom with known geometry and refractive index
distribution. We demonstrate this method by fabricating a multiple-scattering
phantom and comparing 3D refractive index results that are output from three
TPM reconstruction methods operating with visible and near-infrared
wavelengths. One of these methods assumes the sample to be weak-scattering,
while the other two take multiple scattering into account. This study can be
readily extended to more complex microphantoms fabricated to more closely
capture scattering characteristics of real-world scattering objects, such as
tissue.

##  1 Introduction

Tomographic phase microscopy (TPM) is a quantitative, label-free imaging
method that utilizes optical projections through a semi-transparent sample
along various illumination angles to reconstruct the sample’s 3D refractive
index (RI). This method has found several applications in biological imaging,
where RI is directly related to the dry mass distribution at the cellular and
subcellular levels. Refractive index and dry mass are known to be crucial
factors in analyzing the current stage of cell cycle [1], cell structure [2,
3], photobiochemical effects on cells [4], influence of external factors on
cellular parameters [5, 6] and many others. However, to reconstruct 3D RI,
traditional TPM methods utilize critical assumptions in their computational
reconstruction methodologies that rely on the sample being weakly scattering
[7]. These assumptions limit samples to having thicknesses on the order of
only tens of microns. Given the wealth of information provided by dry mass at
the single-cell level, there is significant demand to analyze dry mass in
large multicellular clusters, thick tissue slices, or whole microorganisms.
For these types of samples, reconstruction frameworks that accommodate for
multiple scattering must be utilized. To this end, numerous TPM approaches
have been proposed in recent years that introduce new frameworks to
accommodate multiple scattering [8, 9, 10, 11, 12, 13]. Notably, these
approaches utilize nonlinear and nonconvex solvers to iteratively solve for a
sample’s 3D RI. Though these methods have demonstrated impressive results in
reconstructing RI in multiple-scattering samples, their quantitative accuracy
has not yet been robustly characterized experimentally, and the presented
results usually do not allow comparison of different methods in order to
select proper approach for a given scattering level in a sample. The general
strategy to experimentally evaluate quantitative accuracy in TPM methods is to
reconstruct 3D RI in samples with known RI distributions [14, 15, 16].
Unfortunately, existing microphantoms are typically either weak-scattering
(e.g. index-matched microspheres) or overly simplistic (e.g. index-mismatched
microspheres) compared to the types of heterogenously scattering multicellular
samples that multiple-scattering methods are intended for. This is a critical
limitation when characterizing TPM reconstruction methods that utilize
nonconvex solvers, where iterative convergence depends on the complexity of
the energy landscape and directly associates with a sample’s 3D complexity
[17]. A multiple-scattering TPM method that outputs accurate 3D RI
reconstructions of a weakly scattering or overly simplistic microphantom
cannot be expected to output similarly accurate RI reconstructions for more
complex multiple-scattering samples, where the probability of converging to
local minima are drastically higher. To robustly characterize the quantitative
accuracy of multiple-scattering TPM methods, it is imperative to use
microphantoms with known 3D RI that mimic the structural complexity of the
types of samples that the TPM methods are intended to image. To the best of
our knowledge, these types of gold-standard multiple-scattering microphantoms
do not exist.

In this work we present a 3D-printed microphantom with internal multiple-
scattering RI distribution. To do so, we leverage recent developments in 3D
printing via direct laser writing [18, 19, 20, 21]. We specifically utilize a
two-photon polymerization technique that enables 3D printing of microphantom
samples with known geometry and calibrated RI. Additionally, we use these
microphantom samples to quantitatively assess the RI reconstruction quality in
three TPM reconstruction algorithms that have demonstrated impressive
biological imaging results in prior works.

##  2 Methods

We present below the 1) methodology and design with which we 3D-print multiple
scattering microphantoms; 2) the optical design of the TPM imaging systems
that we use to experimentally collect scattering electric-field measurements
of the microphantom; and 3) short theoretical descriptions of three
tomographic algorithms that were used to reconstruct 3D RI from the measured
data.

###  2.1 Design and manufacturing of the 3D scattering microphantom

The phantom is fabricated using two-photon laser lithography, in which a
focused laser beam is scanned within liquid resin. The resin within the
laser’s focal volume is locally polymerized. Adjusting the scanning trajectory
and the exposure time of the laser beam enables simultaneous control over the
3D printed geometry (accuracy at the order of 100 nm) and RI (accuracy at the
order of $5\cdot 10^{-4}$, maximal $\Delta$RI = 0.03 within the structure) in
three dimensions. We used Photonic Professional GT (Nanoscribe GmbH) equipped
with a 1.3 NA 100x microscope objective and piezo scanning stage. The phantom
is fabricated in the IP-Dip resin (Nanoscribe GmbH) on top of a #1.5H
coverslip (dip-in configuration [22]). After fabrication the structure was
developed in PGMEA (Propylene glycol monomethyl ether acetate; 12 min),
followed by isopropyl alcohol (10 min) and then blow-dried.

For the purposes of demonstration, the particular microphantom we 3D printed
here consists of a cell-like target with internal test structures [23]
embedded within a pseudo-random distribution of rods that switch their
orientation across various layers (Fig. 1a). The width and height of each rod
is equal to $0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$1.8\text{\,}\mathrm{\SIUnitSymbolMicro m}$ respectively. The lateral distance
between the rods in each layer is randomized between
$0.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$ to
$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and the layers are stacked
vertically every $1.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The resulting
structure is transparent (over 99% transmittance for the extinction
coefficient of $0.1\text{\,}\mathrm{m}\mathrm{m}^{-1}$ [24]) and multiple
scattering [25]. The final scattering region is a
$$60\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$60\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ cube with a fill-factor
of $\sim$25% (fraction of volume occupied by the polymer).

Figure 1: a) Half-section view of the scattering microphantom design. b) SEM
image of a scattering layer made out of quasi-randomly distributed rods. c)
Individual rod that comprise a scattering layer. d) Visualization of the 3D RI
distribution of the imaging target – cell phantom. Subcellular features, such
as resolution targets (shown in the inset) and cell nucleoli, are enclosed in
truncated ellipsoid with the external dimensions of
$$30\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$25\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$12\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, which is then embedded
in the center of the $$60\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$60\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ scattering cube.

The cell-like target embedded within the layers of randomly distributed rods
mimics a single biological cell encased within a scattering cube. The cell-
target comprises of substructures that enable assessment of quantitative
accuracy in TPM imaging systems [26]. Methodology for fabrication and
validation of these features can be found in our previous work [23]. The main
features within the cell target include resolution test targets, nucleoli
suspended in a nucleus, and a region of slow RI variation (see Fig. 1d).
Notably, the resolution test target comprises of lines with increasing spatial
frequency [27] up to 1667 lp/mm. By assessing the maximum spatial frequency of
lines that can be distinguished within the cell test target, the imaging
resolution of a TPM method of choice can be characterized. The pseudo-random
distribution of rods that compose the scattering portion of the whole
microphantom are suppressed within $0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$
of the cell target and do not intersect with any of the test structures. To
conduct our TPM imaging experiments, the microphantom was immersed in Zeiss
Immersol 518F oil (RI${}_{$632\text{\,}\mathrm{nm}$}$ = 1.5123), which
provides similar RI contrast as in the case of cells immersed in culture
medium. By using immersion oils with varying RI, it is possible to adjust the
scattering properties of the microphantom post-fabrication.

###  2.2 Experimental setup

In this work, an optical system, as shown in Fig. 2a) was used in order to
study the scattering phantom. The system is a Mach-Zehnder-based TPM
microscope[28], working in a limited-angle configuration with stationary
sample and illumination rotated with a galvo mirror (Thorlabs GVS212/M) [29].
The research was performed with two wavelengths and thus there were two
modified versions of the presented tomographic microscope. First version,
$\text{TPM}_{633}$ works with wavelength $\lambda=$ $633\text{\,}\mathrm{nm}$
and the second, $\text{TPM}_{835}$ with $\lambda=$ $835\text{\,}\mathrm{nm}$.
The input beam (S in Fig.2a) is delivered with an optical fiber, collimated
and then split into object and reference arm. In the $\text{TPM}_{633}$ system
the light source was a volume Bragg grating laser (Necsel NovaTru Chroma 633
SLM), $\text{S}_{633}$, providing a single longitudinal mode and offering long
coherence length. The $\text{TPM}_{835}$ system utilized a swept source
(Superlum Broadsweeper BS-840-2-HP, $\Delta\lambda$= 800-870 nm),
$\text{S}_{835}$ set at $\lambda=$ $835\text{\,}\mathrm{nm}$. Due to
difference in coherence length, an additional delay module was placed in the
reference beam for $\text{TPM}_{835}$ measurements. The beam-splitting cubes
in this work were either coated for 400-700nm or 700-1100nm depending on the
wavelength used. The focal length of the tube lens TL1 was
$\text{EFL}_{633}\text{=}$ $150\text{\,}\mathrm{mm}$ and
$\text{EFL}_{835}\text{=}$ $200\text{\,}\mathrm{mm}$ respectively. Both
microscope objectives (MO1 and MO2) in Fig. 2 were 100x NA 1.3 Semi plan-
apochromatic, infinity-corrected objectives. The second tube lens TL2 used was
either $\text{EFL}_{633}\text{=}$ $200\text{\,}\mathrm{mm}$ or
$\text{EFL}_{835}\text{=}$ $300\text{\,}\mathrm{mm}$. This provided
magnifications $M_{633}=-48.5$ and $M_{835}=-72.7$. The camera used in the
system was a CMOS sensor in both cases, with 3.45$\mu m$ pixel size (JAI
BM500GE) in case of $\text{CAM}_{633}$ and 5.5$\mu m$ pixel size (Basler
acA2040-180km) in case of $\text{CAM}_{835}$. The minimum magnification, which
is imposed by pixel size and wavelength in order to assure correct hologram
recording for each projection [30] is $M_{633min}=-44.2$ and
$M_{835min}=-53.5$, which is satisfied in both cases. A sample hologram is
presented in Fig.2b). Both systems were set to illuminate the sample with a
circular scanning scenario (see Fig. 2d ) at zenith angle $\theta=47$°and
provided 180 projections spaced at $\varphi=2$°. A sample of the phase and
amplitude provided by the system at 835nm is presented in Fig.2 c) and e).

Figure 2: a) Mach-Zehnder-based TPM measurement system. S - light source, GM -
galvo system, TL1 and TL2 - tube lenses, MO1 and MO2 - microscope objective,
SPL - sample plane, CAM - camera; b) hologram acquired at axial illumination
of the sample; c) amplitude of a projection at $\varphi=304$°; d) circular
scanning scenario used in the measurement. Rotation angle: $\varphi$, zenith
angle: $\theta=47$; e) phase of the projection.

###  2.3 Reconstruction algorithms

We compare three different TPM reconstruction methods for computing 3D RI of
the microphantom from the multi-angle scattering measurements captured as
described above. We provide a short description of these methods below.
Complete descriptions of these methods is given in respective references.

####  2.3.1 Gerchberg-Papoulis with object support

To provide a baseline standard to compare to TPM reconstruction algorithms
utilizing multiple-scattering models, we first reconstruct the microphantom’s
3D RI using a weak-scattering TPM method. We specifically use the Gerchberg-
Papoulis algorithm enhanced with additional finite object support constraint
(GPSC) [31]. Complex-valued electric-fields measured by our TPM systems are
used as inputs. The procedure is performed in two steps. First, an initial
tomographic 3D RI distribution is reconstructed from the electric-field
measurements with strong total-variation regularization [32]. This is
performed through the Chambolle-Pock [33] optimization method and implemented
with ASTRA tomography toolbox [34]. The result undergoes binarization and a
finite object support is generated. Secondly, a classic Gerchberg-Papoulis
algorithm is used which is an iterative version of Direct Inversion method
(also known as the Wolf transform) [35]. This iterative procedure is based on
the Fourier Diffraction Theorem [36] and utilizes first-order scattering
approximation. Here, the reconstruction and its Fourier transform are
calculated alternately and constraints are applied: nonnegativity and finite
object support in the signal domain, and replenishment of original projections
in the frequency domain.

####  2.3.2 Multi-slice beam-propagation with electric-field measurements

We utilize the multi-slice beam-propagation method (MSBP) as our main method
to model multiple-scattering [37], which has recently shown promising results
for biological imaging [38, 39]. In our first implementation of MSBP, we use
the same exact electric-field dataset utilized by GPSC from above. An initial
guess of the sample’s 3D RI is selected to start off the iterative procedure.
Afterwards, the MSBP method is used to simulate scattering measurements
resulting from plane waves propagating through the sample’s initial estimated
3D RI. The scattered fields resulting from this simulation are compared with
those obtained experimentally with our TPM systems. The error computed between
the simulated and experimental measurements is back-propagated through each
layer of the 3D sample estimate to incrementally modify the RI value of each
voxel. Continued iterations repeating these steps eventually result in the 3D
sample estimate converging to a stable steady-state solution. We implemented
the MSBP with electric-field measuemerements (MSBP-E) through the Learnining
Tomography algorithm (LT) [38]. The LT procedure is an iterative optimization
algorithm with additional weak total-variation (TV) regularization applied in
each iteration to ensure convergence. We found that the method works best when
an initial-guess is chosen as a starting point for the iterative process. In
this paper, we use the Direct Inversion method to provide the initial guess.

####  2.3.3 Multi-slice beam-propagation with intensity-only measurements

Recent works have demonstrated that the gradient-update step within the MSBP
method can be reformulated to reconstruct 3D RI from only non-interferometric
intensity measurements [13]. The key advantages of this method include the use
of a non-interferometric imaging system, which are resistant to mechanical
instabilities that often limit long-term use of dual-arm interferometers
without realignment. Furthermore, the light source can be partially coherent,
to avoid coherent speckle artifacts in the measurements, while retaining
sufficient coherence necessary for RI reconstruction. For the purposes of
demonstrating 3D RI reconstruction using this intensity-only variant of MSBP
(which we refer to here as MSBP-I), we simply use the amplitude component of
the electric-field measurements used for the GPSC and MSBP-E reconstructions,
described above. Similarly to MSBP-E, total-variation regularization is
applied at every iteration. The starting point for MSBP-I is a matrix of
zeros.

####  2.3.4 Stopping criterion

All described TPM reconstruction methods are iterative procedures that use the
same stopping criterion to automatically terminate the computations. This
criterion is a modification of a method presented earlier [40], and is
described in Alg. 1 below. The general intuition behind this procedure is to
terminate the iterative computation process when the dynamics of the change
between 3D sample estimates outputted from consecutive iterations drops below
a certain saturation level $\epsilon$. In order to be less dependent on
outliers, the median value of the dynamics from the last 10 iterations is
calculated. The value for $\epsilon$ is chosen empirically for each algorithm.
For GPSC $\epsilon=0.02$, for LT and MS $\epsilon=0.01$.

Algorithm 1 Stopping criterion

for $i=1:\textit{NIter}$ do $\triangleright$ NIter \- max number of iterations

Calculate reconstruction $x_{i}$ with a given method

if $i>2$ then

$d_{i}=||x_{i}-x_{i-1}||_{2}/||x_{i}||_{2}$ $\triangleright$ $d_{i}$ \-
"distance" between i-th and (i-1)-th reconstruction

end if

if $i\geq 10$ then

$p_{i}=|d_{i}-d_{i-1}|/|\widehat{d(3:10)}|$ $\triangleright$ $p_{i}$ \-
dynamics of "distance"

end if$\triangleright$ $\widehat{(...)}$ \- median value

if $\widehat{p(\text{end}-11:\text{end})}<\epsilon$ then $\triangleright$
$\epsilon$ \- saturation level

$i=\textit{NIter}$

end if

end for

##  3 Results

Two sinograms were composed from the complex-valued scattered field
measurements of the microphantom being illuminated at different angles, using
both $633\text{\,}\mathrm{nm}$ and $835\text{\,}\mathrm{nm}$ wavelengths. For
GPSC and MSBP-E, these sinograms were used directly. For MSBP-I, only the
amplitude components of these sinograms were used. The sinograms and their
corresponding 3D RI reconstructions have been made open-source and are
available in Dataset 1 [41]. Figure 3 shows the reconstruction results for
both wavelengths. We characterize lateral (x and y) resolution of the TPM
reconstructions by visualizing the resolution test lines inside the
microphantom. To quantitatively compare the three TPM reconstruction methods,
horizontal and vertical cross-sectional plots across these resolution tests
were generated by computing the average and standard-deviation of the pixel-
values across rows or columns adjacent to the dashed white a-a and b-b lines,
by $\pm 4$ pixels. These cross-sectional plots are shown below.

Figure 3: Comparison of tomographic reconstructions of the microphantom
measured with $633\text{\,}\mathrm{nm}$ and $835\text{\,}\mathrm{nm}$
wavelength and calculated with 3 algorithms. The shaded colored regions
surrounding each of the 1D plots at the bottom represent the standard
deviation.

##  4 Discussion

The results show that changing the wavelength from visible into near-infrared
impacts the applicability of TPM reconstruction methods [42, 43]. The
magnified views of the resolution test region after 3D reconstructing the
microphantom using the GPSC algorithm with $633\text{\,}\mathrm{nm}$
wavelength light reveals significant grainy artifacts which occlude the line
features within the microphantom’s test region. We note that these artifacts
are heavily decreased when reconstructing with $835\text{\,}\mathrm{nm}$
wavelength light. This suggests that, at $633\text{\,}\mathrm{nm}$ wavelength,
the microphantom is too highly-scattering for GPSC to be applicable. The
decreased GPSC reconstruction artifacts for $835\text{\,}\mathrm{nm}$ matches
conventional knowledge that longer wavelengths of light are more resistant to
scattering than shorter wavelengths. In order to confirm that this observation
is due to reduced scattering and not different noise characteristics between
the two light sources, we analyzed the standard deviation of phase-noise in an
object-free region within the input sinograms. For $633\text{\,}\mathrm{nm}$
wavelength light, a phase-noise standard deviation of $\sigma=0.10$ radians
was observed, while for $835\text{\,}\mathrm{nm}$ wavelength light, a phase-
noise standard deviation of $\sigma=0.08$ radians was observed. Given such a
small variation between these noise characteristics, we conclude that the
major factor in the GPSC reconstruction quality that affects its ability to
visualize the microphantom’s test lines is the scattering strengths of the
microphantoms at the two different wavelengths.

Notably, both MSBP-E and MSBP-I use total-variation (TV) regularization in
order to stabilize the convergence of the nonconvex iterative solver in the
presence of noise [32]. Especially in the case of using
$633\text{\,}\mathrm{nm}$ wavelength light, TV regularization results in 3D RI
reconstructions with less noise compared to the 3D reconstructions computed
via GPSC, which does not use TV regularization. This can be directly
visualized in the 2D cross-sections in Fig. 3, and is confirmed by the bounds
of standard deviation shown in the 1D cross-sectional plots. However, the
drawback of regularization is that it has a blurring effect on high-resolution
features. Because the microphantom is less scattering in
$835\text{\,}\mathrm{nm}$ wavelength light, only GPSC managed to reconstruct
the high spatial-frequency test lines within the microphantom. Typically, the
strength of regularization is manually tuned to fit experimental factors and
balance between the tradeoff between achieving iterative stability versus high
imaging resolution.

In terms of the average RI, the knowledge about the ground truth RI
distribution of the phantom makes it possible to quantitatively compare
reconstruction results with GPSC, MSBP-E, and MSBP-I. We see that all methods
successfully capture the bulk characteristics of the microphantom. As
described above, the 3D reconstruction via GPSC exhibits grainy artifacts when
using $633\text{\,}\mathrm{nm}$ wavelength light, likely due to the
microphantom being multiple-scattering at that wavelength. Furthermore, MSBP-I
outputs slightly overestimated RI values, and also suffers from low-frequency
spatial artifacts (which have been observed in other intensity-only phase-
imaging techniques [44, 45]). Other works have shown that MSBP-I demonstrates
higher accuracy when using partially-coherent illumination, which drastically
reduces coherent noise [13]. Future work may include repeating this analysis
across a larger range of TPM reconstruction techniques with more complex
scattering microphantoms.

##  5 Conclusions

Here we presented a 3D-printed microphantom with known geometry and RI
distribution, for the purpose of quantitatively assessing reconstruction
quality and accuracy between various TPM reconstruction methodologies. For the
microphantom presented here, the cell-like portion within the center is weakly
scattering. However, this cell-like portion is surrounded with a dense pseudo-
random grid of heterogenous RI structures, which scatters light. This enables
the total microphantom to be an appropriate 3D calibration object to robustly
evaluate the capability of various TPM methodologies in reconstructing 3D
multiple-scattering samples. It should be noted that it is possible to change
the scattering properties of the phantom by changing the design of the grid or
by changing the refractive index of the immersion liquid.

We have shown experimentally that changing the illumination wavelength affects
the scattering nature of the microphantom. Specifically, though the
microphantom is multiple scattering with $633\text{\,}\mathrm{nm}$ wavelength
light, it is weak scattering with $835\text{\,}\mathrm{nm}$ wavelength light.
This naturally indicates that the optimum choice for illumination wavelength
must balance between resolution ($\frac{\lambda}{NA}$ for single projection)
and scattering strength. As has been shown, there are cases when instead of
applying multiple-scattering methods, it is advantageous to increase the
illumination wavelength (thus decreasing the scattering strength of the
sample) and apply a method based on the Fourier Diffraction Theorem that does
not utilize the total-variation constraint. More fundamentally, however, we
showed that 3D printed microphantoms enable quantitative assessment of 3D RI
reconstruction accuracy across various TPM methodologies. This capability is
important when choosing a TPM method optimized for specific classes of samples
and imaging conditions. An important and natural extension of this work is to
develop a methodology to 3D-print complex and realistic microphantoms that
mimic a wide range of biological specimens, ranging from multicellular
clusters to bulk tissues and small organisms.

##  6 Backmatter

###  6.1 Funding

The research leading to the described results was carried out within the
program TEAM TECH/2016-1/4 of Foundation for Polish Science, co-financed by
the European Union under the European Regional Development Fund. Calibration
object development was funded by FOTECH-1 project granted by Warsaw University
of Technology under the program Excellence Initiative: Research University
(ID-UB). We also gratefully acknowledge support from the University of Texas
at Austin, Cockrell School of Engineering, and the Chan Zuckerberg Initiative.

###  6.2 Acknowledgments

The authors want to thank Demetri Psaltis and Joowon Lim from École
Polytechnique Fédérale de Lausanne for providing the code for Learning
Tomography algorithm.

###  6.3 Disclosures

The authors declare no conflicts of interest.

###  6.4 Data availability

Data underlying the results presented in this paper are available in Dataset 1
[41].

###  6.5 CRediT

Wojciech Krauze: Conceptualization, Formal analysis, Methodology, Software,
Visualization, Writing – original draft. Arkadiusz Kuś: Conceptualization,
Funding acquisition, Investigation, Methodology, Project administration,
Visualization, Writing – original draft. Michał Ziemczonok: Data curation,
Investigation, Methodology, Resources, Visualization, Writing – original
draft. Max Haimowitz: Software, Writing – review & editing. Shwetadwip
Chowdhury: Funding acquisition, Software, Writing – review & editing.
Małgorzata Kujawińska: Funding acquisition, Project administration, Writing –
review & editing.

## References

  * [1] Pinhas Girshovitz and Natan T. Shaked.  Generalized cell morphological parameters based on interferometric phase microscopy and their application to cell life cycle characterization.  Biomedical Optics Express, 3(8):1757–1773, August 2012.
  * [2] Youngchan Kim, Hyoeun Shim, Kyoohyun Kim, HyunJoo Park, Ji Han Heo, Jonghee Yoon, Chulhee Choi, Seongsoo Jang, and YongKeun Park.  Common-path diffraction optical tomography for investigation of three-dimensional structures and dynamics of biological cells.  Optics Express, 22(9):10398–10407, 2014.
  * [3] Tess A Stanly, Rakesh Suman, Gulab Fatima Rani, Peter J O’Toole, Paul M Kaye, and Ian S Hitchcock.  Quantitative optical diffraction tomography imaging of mouse platelets.  Frontiers in physiology, 11:568087, 2020.
  * [4] Maria Baczewska, Piotr Stępień, Martyna Mazur, Wojciech Krauze, Natalia Nowak, Jędrzej Szymański, and Małgorzata Kujawińska.  Method to analyze effects of low-level laser therapy on biological cells with a digital holographic microscope.  Applied Optics, 61(5):B297–B306, Feb 2022.
  * [5] Maria Baczewska, Kai Eder, Steffi Ketelhut, Björn Kemper, and Małgorzata Kujawińska.  Refractive index changes of cells and cellular compartments upon paraformaldehyde fixation acquired by tomographic phase microscopy.  Cytometry Part A, 99(4):388–398, 2021.
  * [6] Kai Eder, Anne Marzi, Álvaro Barroso, Björn Kemper, and Jürgen Schnekenburger.  Medical nanoparticle impact on macrophage temporal dry mass development quantified in vitro by digital holographic microscopy.  In Quantitative Phase Imaging VII, volume 11653, pages 132–138. SPIE, 2021.
  * [7] Bingquan Chen and Jakob J Stamnes.  Validity of diffraction tomography based on the first born and the first rytov approximations.  Applied Optics, 37(14):2996–3006, 1998.
  * [8] Moosung Lee, Herve Hugonnet, and YongKeun Park.  Inverse problem solver for multiple light scattering using modified born series.  Optica, 9(2):177–182, 2022.
  * [9] Hsiou-Yuan Liu, Dehong Liu, Hassan Mansour, Petros T Boufounos, Laura Waller, and Ulugbek S Kamilov.  Seagle: Sparsity-driven image reconstruction under multiple scattering.  IEEE Transactions on Computational Imaging, 4(1):73–86, 2017.
  * [10] Yu Sun, Zhihao Xia, and Ulugbek S Kamilov.  Efficient and accurate inversion of multiple scattering with deep learning.  Optics Express, 26(11):14678–14688, 2018.
  * [11] Emmanuel Soubies, Thanh-An Pham, and Michael Unser.  Efficient inversion of multiple-scattering model for optical diffraction tomography.  Optics Express, 25(18):21786–21800, 2017.
  * [12] Michael Chen, David Ren, Hsiou-Yuan Liu, Shwetadwip Chowdhury, and Laura Waller.  Multi-layer born multiple-scattering model for 3d phase microscopy.  Optica, 7(5):394–403, 2020.
  * [13] Shwetadwip Chowdhury, Michael Chen, Regina Eckert, David Ren, Fan Wu, Nicole Repina, and Laura Waller.  High-resolution 3D refractive index microscopy of multiple-scattering samples from intensity images.  Optica, 6(9):1211, September 2019.
  * [14] JooWon Lim, KyeoReh Lee, Kyong Hwan Jin, Seungwoo Shin, SeoEun Lee, YongKeun Park, and Jong Chul Ye.  Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography.  Optics Express, 23(13):16933–16948, 2015.
  * [15] Yongjin Sung and Ramachandra R Dasari.  Deterministic regularization of three-dimensional optical diffraction tomography.  JOSA A, 28(8):1554–1561, 2011.
  * [16] JooWon Lim, Abdul Wahab, GwangSik Park, KyeoReh Lee, YongKeun Park, and Jong Chul Ye.  Beyond born-rytov limit for super-resolution optical diffraction tomography.  Optics Express, 25(24):30445–30458, 2017.
  * [17] Haitham Hindi.  A tutorial on convex optimization.  In Proceedings of the 2004 American Control Conference, volume 4, pages 3252–3265. IEEE, 2004.
  * [18] Andrew C Lamont, Michael A Restaino, Abdullah T Alsharhan, Zhuolin Liu, Daniel X Hammer, Ryan D Sochol, and Anant Agrawal.  Direct laser writing of a titanium dioxide-laden retinal cone phantom for adaptive optics-optical coherence tomography.  Optical Materials Express, 10(11):2757–2767, 2020.
  * [19] Christopher N LaFratta and Tommaso Baldacchini.  Two-photon polymerization metrology: Characterization methods of mechanisms and microstructures.  Micromachines, 8(4):101, 2017.
  * [20] Hannah Horng, Kerry O’Brien, Andrew Lamont, Ryan D Sochol, T Joshua Pfefer, and Yu Chen.  3d printed vascular phantoms for high-resolution biophotonic image quality assessment via direct laser writing.  Optics Letters, 46(8):1987–1990, 2021.
  * [21] Matthias Eifler, Julian Hering, Georg Von Freymann, and Jörg Seewig.  Calibration sample for arbitrary metrological characteristics of optical topography measuring instruments.  Optics Express, 26(13):16609–16623, 2018.
  * [22] Ada-Ioana Bunea, Nuria del Castillo Iniesta, Ariadni Droumpali, Alexandre Emmanuel Wetzel, Einstom Engay, and Rafael Taboryski.  Micro 3D Printing by Two-Photon Polymerization: Configurations and Parameters for the Nanoscribe System.  Micro, 1(2):164–180, sep 2021.
  * [23] Michał Ziemczonok, Arkadiusz Kuś, Piotr Wasylczyk, and Małgorzata Kujawińska.  3D-printed biological cell phantom for testing 3D quantitative phase imaging systems.  Scientific Reports, 9(1):18872, dec 2019.
  * [24] Michael Schmid, Dominik Ludescher, and Harald Giessen.  Optical properties of photoresists for femtosecond 3D printing: refractive index, extinction, luminescence-dose dependence, aging, heat treatment and comparison between 1-photon and 2-photon exposure.  Optical Materials Express, 9(12):4564–4577, dec 2019.
  * [25] Evangelos Marakis, Ravitej Uppu, Maryna L. Meretska, Klaas-Jan Gorter, Willem L. Vos, and Pepijn W. H. Pinkse.  Deterministic and Controllable Photonic Scattering Media via Direct Laser Writing.  Advanced Optical Materials, 8(24):2001438, dec 2020.
  * [26] Michał Ziemczonok, Arkadiusz Kuś, and Małgorzata Kujawińska.  Optical diffraction tomography meets metrology — Measurement accuracy on cellular and subcellular level.  Measurement, 195(January):111106, may 2022.
  * [27] Roarke Horstmeyer, Rainer Heintzmann, Gabriel Popescu, Laura Waller, and Changhuei Yang.  Standardizing the resolution claims for coherent microscopy.  Nature Photonics, 10(2):68–71, 2016.
  * [28] Arkadiusz Kuś.  Illumination-related errors in limited-angle optical diffraction tomography.  Applied Optics, 56(33):9247–9256, 2017.
  * [29] Arkadiusz Kuś, Wojciech Krauze, Piotr L. Makowski, and Małgorzata Kujawińska.  Holographic tomography: hardware and software solutions for 3D quantitative biomedical imaging (Invited paper).  ETRI Journal, 41(1):61–72, February 2019.
  * [30] Emilio Sánchez-Ortiga, Ana Doblas, Genaro Saavedra, Manuel Martínez-Corral, and Jorge Garcia-Sucerquia.  Off-axis digital holographic microscopy : practical design parameters for operating at diffraction limit.  Applied Optics, 53(10):2058–2066, 2014.
  * [31] Wojciech Krauze.  Optical diffraction tomography with finite object support for the minimization of missing cone artifacts.  Biomedical Optics Express, 11(4):1919–1926, 2020.
  * [32] David Strong and Tony Chan.  Edge-preserving and scale-dependent properties of total variation regularization.  Inverse problems, 19(6):S165, 2003.
  * [33] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of mathematical imaging and vision, 40(1):120–145, 2011\.
  * [34] Wim Van Aarle, Willem Jan Palenstijn, Jan De Beenhouwer, Thomas Altantzis, Sara Bals, K Joost Batenburg, and Jan Sijbers.  The astra toolbox: A platform for advanced algorithm development in electron tomography.  Ultramicroscopy, 157:35–47, 2015.
  * [35] Emil Wolf.  Three-dimensional structure determination of semi-transparent objects from holographic data.  Optics communications, 1(4):153–156, 1969.
  * [36] O. Haeberlé, K. Belkebir, H. Giovaninni, and A. Sentenac.  Tomographic diffractive microscopy: basics, techniques and perspectives.  Journal of Modern Optics, 57(9):686–699, 2010.
  * [37] J Van Roey, J Van der Donk, and PE Lagasse.  Beam-propagation method: analysis and assessment.  Josa, 71(7):803–810, 1981.
  * [38] Ulugbek S. Kamilov, Ioannis N. Papadopoulos, Morteza H. Shoreh, Alexandre Goy, Cedric Vonesch, Michael Unser, and Demetri Psaltis.  Learning approach to optical tomography.  Optica, 2(6):517, 2015.
  * [39] Ulugbek S Kamilov, Ioannis N Papadopoulos, Morteza H Shoreh, Alexandre Goy, Cedric Vonesch, Michael Unser, and Demetri Psaltis.  Optical tomographic image reconstruction based on beam propagation and sparse regularization.  IEEE Transactions on Computational Imaging, 2(1):59–70, 2016.
  * [40] Piotr L. Makowski and Michał Ziemczonok.  Projection extrapolation routine for tight-frame limited-angle optical diffraction tomography.  Optics Letters, 44(14):3442–3445, Jul 2019.
  * [41] Wojciech Krauze, Arkadiusz Kuś, Michał Ziemczonok, Max Haimowitz, Schwetadwip Chowdhury, and Małgorzata Kujawińska.  Optical tomography measurements and reconstructions of a multiple-scattering 3d-printed microphantom.  Zenodo, 2022.  https://doi.org/10.5281/zenodo.6652333.
  * [42] Fritjof Helmchen and Winfried Denk.  Deep tissue two-photon microscopy.  Nature methods, 2(12):932–940, 2005.
  * [43] Michael Rubart.  Two-photon microscopy of cells and tissue.  Circulation research, 95(12):1154–1166, 2004.
  * [44] Zhuo Wang, Larry Millet, Mustafa Mir, Huafeng Ding, Sakulsuk Unarunotai, John Rogers, Martha U Gillette, and Gabriel Popescu.  Spatial light interference microscopy (slim).  Optics Express, 19(2):1016–1026, 2011.
  * [45] Lei Tian and Laura Waller.  Quantitative differential phase contrast imaging in an led array microscope.  Optics Express, 23(9):11394–11403, 2015.


# Quantile Extreme Gradient Boosting for Uncertainty Quantification ††thanks:
Citation: Authors. Title. Pages…. DOI:000000/11111.

Xiaozhe Yin, Masoud Fallah-Shorshani, Rob McConnell, Scott Fruin
Department of Population and Public Health Sciences, Keck School of Medicine
University of Southern California
Los Angeles, CA 90032 USA
{xiaozhey, shorshan, rmcconne<EMAIL_ADDRESS>
&Yao-Yi Chiang
Department of Computer Science & Engineering
University of Minnesota
Minneapolis, MN 55455 USA
<EMAIL_ADDRESS>
Meredith Franklin
Department of Statistical Sciences and School of Environment
University of Toronto
Toronto Ontario, Canada
<EMAIL_ADDRESS>

###### Abstract

As the availability, size and complexity of data have increased in recent
years, machine learning (ML) techniques have become popular for modeling.
Predictions resulting from applying ML models are often used for inference,
decision-making, and downstream applications. A crucial yet often overlooked
aspect of ML is uncertainty quantification, which can significantly impact how
predictions from models are used and interpreted.

Extreme Gradient Boosting (XGBoost) is one of the most popular ML methods
given its simple implementation, fast computation, and sequential learning,
which make its predictions highly accurate compared to other methods. However,
techniques for uncertainty determination in ML models such as XGBoost have not
yet been universally agreed among its varying applications. We propose
enhancements to XGBoost whereby a modified quantile regression is used as the
objective function to estimate uncertainty (QXGBoost). Specifically, we
included the Huber norm in the quantile regression model to construct a
differentiable approximation to the quantile regression error function. This
key step allows XGBoost, which uses a gradient-based optimization algorithm,
to make probabilistic predictions efficiently.

QXGBoost was applied to create 90% prediction intervals for one simulated
dataset and one real-world environmental dataset of measured traffic noise.
Our proposed method had comparable or better performance than the uncertainty
estimates generated for regular and quantile light gradient boosting. For both
the simulated and traffic noise datasets, the overall performance of the
prediction intervals from QXGBoost were better than other models based on
coverage width-based criterion.

_K_ eywords uncertainty quantification $\cdot$ xgboost $\cdot$ quantile
regression $\cdot$ huber function

##  1 Introduction

Modeling complex data such as those arising from environmental processes is
challenging given the mix of spatial scales, interactions between components,
and underlying spatiotemporal patterns. With the explosion of geospatial data
from a variety of sources, including satellites, monitoring networks, and
field campaigns, environmental scientists have been adopting data-driven
approaches rather than relying solely on more traditional deterministic
techniques such as earth system (climate), dispersion, and acoustic modeling.
To handle large data sets, artificial intelligence (AI) and machine learning
(ML) have become increasingly popular, as simple classification or linear
regression methods do not provide the flexibility that allow for non-
linearities and complex interactions [1]. Nevertheless, though AI and ML
methods such as neural networks and gradient boosting algorithms have achieved
high performance in point predictions, they are often still operated as black
boxes. More concerningly, predictions from these models have been shown to be
miscalibrated, and can yield arbitrarily high confidence [2]; they also
usually do not come with any estimation of uncertainty, an important gap. A
recent report by the National Academy of Sciences indicated that to make
advances at the intersection of AI/ML and Earth system science, understanding
and communicating uncertainty is essential for decision-making [3].

Studies have been conducted that quantify uncertainty in deep learning, among
which Bayesian methods [4, 5], the delta method [6], bootstrap [7, 8], and
lower upper bound estimation [5] are the most commonly used approaches. These
uncertainty quantification methods typically give a probability density
function as the output, which can be unsatisfactory. For example, some studies
argued that the PIs of the bootstrap method were large while coverage was poor
([9, 10, 11]). Unfortunately these methods are difficult to generalize to
tree-based methods, which is why other approaches need to be developed.

Quantile regression for uncertainty quantification has recently gained
popularity in gradient-based optimization algorithms, where it is introduced
to construct point-wise prediction intervals (PIs) of specific quantiles. For
example, Taylor (2000) [12] used a quantile regression neural network (QRNN)
to estimate the multiple quantiles of multi-period financial returns. QRNN has
also been used in predicting wind power [13], electricity consumption [14],
rainfall extremes [15], and to evaluate online teaching [16]. Xu et al. (2017)
[17] developed a composite quantile regression neural network and applied it
to three standard datasets, showing that it had better performance than QRNN.

Quantile regression has a solid statistical foundation and can be easily
customized and incorporated into the cost function of many ML models. For
example, quantile regression random forest (QRRF) has been developed, and some
studies have suggested that this method is competitive in terms of its
predictive power [18]. There have been applications of QRRF in statistical
load forecasting [19], uncertainty estimation of digital soil mapping products
[20] and suspended sediment concentration [21]. Quantile regression has also
been implemented in gradient-boosting-tree-based methods, such as gradient
boosting machine (GBM) [22] and light gradient boosting machine (LightGBM).
Quantile GBM has been used in environmental probabilistic forecasting
problems, including to forecast wind power [23] and solar irradiation [24].
Both studies found that including quantiles significantly improved the models
by enabling prediction of quantiles for uncertainty quantification. As an
emerging machine learning method, LightGBM has been applied to predict
particulate matter air pollution from satellite observations [25]. However, we
have not seen an application of the quantile LightGBM.

Despite quantile regression gaining popularity in neural networks and some
tree-based machine learning methods, it has never been used in extreme
gradient boosting (XGBoost) for two reasons. First, the quantile regression
function is not differentiable at 0, meaning that the gradient-based XGBoost
method might not converge properly and lead to sub-optimal model parameters
[26]. Second, different from GBM, a second-order approximation is required to
optimize the regularized objective function in XGBoost, and the second
derivative of the quantile regression function is 0 everywhere. To solve these
limitations, we incorporated the use of the Huber norm in combining the
quantile regression method with XGBoost and we tested the proposed quantile
XGBoost method in both a simulated and a real-world traffic noise dataset.
Even though we used environmental data as an example in this study given the
complexity in geospatial data, QXGBoost can be universally applied to any type
of data for which uncertainty analysis is desired.

##  2 Extreme Gradient boosting

Built upon the widely used tree-based boosting method, XGBoost was proposed by
Chen and Guestrin (2016) [27] to use additive linear functions to predict the
output as shown in equation 1.

$\hat{y_{i}}=\phi(x_{i})=\sum_{k=1}^{k}f_{k}(x_{i}),~{}f_{k}\in\mathcal{F}$
(1)

where $\mathcal{F}=\\{f(x)=w_{q(x)}\\}(q:\mathbb{R}^{m}\to
T,w\in\mathbb{R}^{T})$ is the space of regression trees and $q$ is the
distinctive structure of each tree that is represented by $f(k)$ to map an
input to the corresponding leaf index; $K$ is the number of trees and $T$ is
the number of leaves in each tree; $w$ represents the leaf weights (scores).
The final prediction of a specific input is calculated by summing up the
scores ($w$) in the corresponding leaves.

As with all supervised machine learning models, XGBoost uses an objective
function to train the data and assess the model fit. Instead of using just the
loss function, XGBoost adds an additional regularization term to prevent
overfitting. For regression problems, mean squared error is a common choice
for loss function and the corresponding objective function can be expressed
as:

$\displaystyle\mathcal{L}(\theta)^{(t)}$
$\displaystyle=L(\theta)+\Omega(\theta)$ (2)
$\displaystyle=\sum_{i}l(\hat{y_{i}},y_{i})+\sum_{k}\Omega(f_{k})$
$\displaystyle=\sum_{i=1}^{n}(y_{i}-(\hat{y_{i}}^{(t-1)}+f_{t}(x_{i})))^{2}+\Omega(f_{t})$

Here $L$ represents the differentiable convex loss function that measures the
difference between the prediction $\hat{y_{i}}$ and the target $y_{i}$. The
second term $\Omega$ is the regularization term that penalizes the complexity
of the model, and $\Omega(f)=\gamma T+\frac{1}{2}\lambda w^{2}$, where $T$
refers to the number of leaves in each tree; $\gamma$ and $\lambda$ are the
$L1$ and $L2$ regularization terms that control weights at leaves.

The loss function for regression problems can easily be calculated using
first-order and quadratic terms (equation 2). However, the loss functions for
other losses of interest, such as the logistic loss, are not friendly. To
overcome this issue, XGBoost employs the Taylor expansion, originally proposed
by Friedman et al. (2000) [28] to quickly optimize the objective function in
equation 2 to the second order. The Taylor expansion of function $f(x)$ that
is infinitely differentiable at a real number $a$ can be written as:

$\displaystyle f(x)$
$\displaystyle=f(a)+\frac{f^{{}^{\prime}}(a)}{1!}(x-a)+\frac{f^{{}^{\prime\prime}}(a)}{2!}(x-a)^{2}+\frac{f^{{}^{\prime\prime\prime}}(a)}{3!}(x-a)^{3}+...$
(3) $\displaystyle=\sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^{n}$

In practice, the higher orders of the Taylor expansion function are usually
omitted to get the approximation of the function. XGBoost keeps the constant,
first, and second derivatives of $f(x)$ evaluated at $a$. By replacing $f(x)$
with $(f_{t}(x_{i})+\hat{y_{i}}^{(t-1)})-y_{i}$ and expanding
$(f_{t}(x_{i})+\hat{y_{i}}^{(t-1)})-y_{i}$ at $(\hat{y_{i}}^{(t-1)})-y_{i})$
in equation 2, the objective function can be rewritten as:

$\mathcal{L}^{(t)}\simeq\sum_{i=1}^{n}[l(y_{i},\hat{y_{i}}^{(t-1)})+g_{i}{f_{t}(x_{i})}+\frac{1}{2}h_{i}f_{t}^{2}(x_{i})]+\Omega(f_{t})$
(4)

where the $g_{i}$ and $h_{i}$ are the first and second order gradient
statistics of the loss function and they are defined as

$\displaystyle
g_{i}=\partial_{\hat{y_{i}}^{(t-1)}}l(y_{i},\hat{y_{i}}^{(t-1)})$ (5)
$\displaystyle
h_{i}=\partial_{\hat{y_{i}}^{(t-1)}}^{2}l(y_{i},\hat{y_{i}}^{(t-1)})$

In XGBoost, only the functions $f_{t}(x_{i})$ are treated as parameters. Thus,
the $l(y_{i},\hat{y_{i}}^{(t-1)})$ can be considered as constant terms. Given
that all the instances on the same leaf get the same score, if we define
$I_{j}=\\{i|q(x_{i})=j\\}$ as the set of indices of data points assigned to
the $j$-th leaf, $G_{i}=\sum_{i\in I_{j}}g_{i}$ and $H_{i}=\sum_{i\in
I_{j}}h_{i}$, and substitute the regularization term $\Omega(f_{t})$, the
objective function of equation LABEL:eq:g_i_and_h_i can be further compressed
as

$\displaystyle\mathcal{\tilde{L}}^{(t)}$
$\displaystyle\simeq\sum_{i=1}^{n}[g_{i}{f_{t}(x_{i})}+\frac{1}{2}h_{i}f_{t}^{2}(x_{i})]+\gamma
T+\frac{1}{2}\lambda w^{2}$ (6) $\displaystyle=\sum_{j=1}^{T}[(\sum_{i\in
I_{j}}g_{i})w_{j}+\frac{1}{2}(\sum_{i\in I_{j}}h_{i}+\lambda)w_{j}^{2}]+\gamma
T$
$\displaystyle=\sum_{j=1}^{T}[G_{j}w_{j}+\frac{1}{2}(H_{j}+\lambda)w_{j}^{2}]+\lambda
T$

For a fixed structure, the optimal weight $w_{j}^{*}$ of leaf $j$ can be
obtained by

$w_{j}^{*}=-\frac{G_{j}}{H_{j}+\lambda}$ (7)

The corresponding objective function when plugging $w_{j}^{*}$ back to
equation 6 is

$\mathcal{\tilde{L}}^{(t)}=-\frac{1}{2}\sum_{j=1}^{T}\frac{G_{j}^{2}}{H_{j}+\lambda}+\lambda
T$ (8)

In practice, function 8 was used to evaluate the tree structure. The smaller
the score, the better the tree structure. From equation 8, we can see that the
value of the objective function only depends on $g_{i}$ and $h_{i}$, meaning
that any function with $g_{i}$ and $h_{i}$ can be used as a potential
objective function. Thus, XGBoost supports the customized objective function.
XGBoost uses a greedy algorithm to search for the best split of a tree.
Specifically, the following score gain function is used to split a leaf into
two leaves

$\mathcal{L}_{splitgain}=\frac{1}{2}[\frac{G_{L}^{2}}{H_{L}+\lambda}+\frac{G_{R}^{2}}{H_{R}+\lambda}-\frac{(G_{L}+G_{R})^{2}}{H_{L}+H_{R}+\lambda}]-\gamma$
(9)

where $I_{L}$ and $I_{R}$ are the instance sets of the new left and right leaf
after the split and $I=I_{L}+I_{R}$ is the score on the original leaf. The
larger the split score, the better the split. The best split of a leaf happens
when the split gain score is the largest among all possible split candidates.
For equation 9, if the split gain scores are all smaller than $\gamma$, which
means that the split score is negative, then there will be no split at that
leaf. Equation 9 is usually used in practice to evaluate all split candidates.

##  3 Quantile regression

Quantile regression was first introduced by Koenker and Bassett Jr (1978) [29]
to estimate and conduct inference about conditional quantile functions.
Instead of minimizing the least square, quantile regression minimizes the mean
absolute error, which applies asymmetric weights to positive and negative
errors to get the entire conditional distribution [29]. Thus, quantile
regression is the process of changing the mean square error loss function to
one that predicts conditional quantiles rather than conditional means [30].
The absolute value function, which is also called the check, tick, or pinball
loss function, is given by

$\rho_{\tau}(x)=\begin{cases}(\tau-1)x&x<0\\\ \tau x&otherwise\end{cases}$
(10)

where $\tau$ is a specific quantile that we want to calculate ($0<\tau<1$) and
$x$ represents the prediction error $y_{i}-\hat{y_{i}}$. Given that the mean
absolute error is $\frac{1}{n}\sum_{i=1}^{n}|y_{i}-\hat{y_{i}}|$, the
corresponding mean quantile loss can be expressed as
$\frac{1}{n}\sum_{i=1}^{n}\rho_{q}(y_{i}-\hat{y_{i}})$. To prove that the
optimal estimate of this loss function is actually the quantiles, we calculate
the expectation of the absolute value function as

$\displaystyle E[\rho_{\tau}(y-\hat{y})]$
$\displaystyle=(\tau-1)\int_{-\infty}^{\hat{y}}f(t)(t-\hat{y})dt+\tau\int_{\hat{y}}^{\infty}f(t)(t-\hat{y})dt$
(11) $\displaystyle=\tau\int
f(t)(t-\hat{y})dt-\int_{-\infty}^{\hat{y}}f(t)(t-\hat{y})dt$
$\displaystyle=\tau\int f(t)tdt-\hat{y}\tau\int
f(t)dt-\int_{-\infty}^{\hat{y}}f(t)tdt+\hat{y}\int_{-\infty}^{\hat{y_{i}}}f(t)dt$
$\displaystyle=\tau\int
f(t)tdt-\hat{y}\tau-\int_{-\infty}^{\hat{y}}f(t)tdt+\hat{y}F(\hat{y})$

Here $f$ and $F$ represent the probability density and the cumulative
distribution functions of $Y$, respectively. If we differentiate the equation
11 in terms of $\hat{y}$ and set the expression equal to 0, we get

$\displaystyle\frac{\partial E[\rho_{\tau}(y-\hat{y})]}{\partial\hat{y}}=0$
$\displaystyle=\frac{\partial(\tau\int
f(t)tdt-\hat{y}\tau-\int_{-\infty}^{\hat{y}}f(t)tdt+\hat{y}F(\hat{y}))}{\partial\hat{y}}$
(12)
$\displaystyle=\frac{\partial(\hat{y}F(\hat{y}))}{\partial\hat{y}}+\frac{\partial(\tau\int
f(t)tdt)}{\partial\hat{y}}-\frac{\partial(\hat{y}\tau)}{\partial\hat{y}}-\frac{\partial\int_{-\infty}^{\hat{y}}f(t)tdt}{\partial\hat{y}}$
$\displaystyle=\frac{F(\hat{y})\partial\hat{y}}{\partial\hat{y}}+\frac{\hat{y}\partial
F(\hat{y})}{\partial\hat{y}}+0-\tau\frac{\partial\hat{y}}{\partial\hat{y}}-\hat{y}f(\hat{y})$
$\displaystyle=F(\hat{y})+\hat{y}f(\hat{y})-\tau-\hat{y}f(\hat{y})$
$\displaystyle=F(\hat{y})-\tau$

Solving the above equation with respect to $\hat{y}$, we get
$\hat{y}=F^{-1}(\tau)$, which is the specific quantile that we want to
calculate. Equation 12 proves that the prediction $\hat{y}$ is the $\tau^{th}$
quantile. The absolute value function is a non-negative function as shown in
Figure 1

Figure 1: Quantile plot.

##  4 Quantile extreme gradient boosting

The proposed quantile extreme gradient boosting (QXGBoost) method combines
quantile regression and XGBoost to construct prediction intervals (PIs).
Specifically, instead of using the mean square error as the loss function in
XGBoost for regression problems, QXGBoost employs the absolute value function
as the loss function to calculate the PIs, which quantify uncertainties. By
doing so, QXGBoost can create specific prediction quantiles that can be
customized depending on the application area.

In XGBoost, the optimal weights (equation 7), structure score (equation 8),
and the split gain function (equation 9) all depend on the gradient statistics
of each instance. Thus, the objective function in XGBoost must meet two
requirements. First, it must be continuously differentiable in the defined
sample space. Second, the first and second derivatives of the loss function
should exist and not be 0 everywhere. However, the absolute value function
(equation 10) is not differentiable at the origin since it has a “break” point
and the second derivative is 0 everywhere (figure 1). Thus, the original
quantile regression function cannot be used directly as the customized
objective function in XGBoost.

To solve the two problems imposed for implementing the quantile regression in
XGBoost, we need a differentiable approximation to the quantile regression
error function so that the gradient-based optimization algorithm can work. We
followed the work by Chen (2007) [31], who introduced the Huber function [32]
to construct the smooth approximation of the absolute value function and thus
has a smooth transition at 0. The Huber norm function has been used as a loss
function in neural networks [33], but this is the first use of it in another
algorithm, including XGBoost. The Huber function is written as:

$h_{\upsilon}(t)=\begin{cases}\frac{t^{2}}{2\upsilon}&|t|\leq\upsilon\\\
|t|-\frac{\upsilon}{2}&|t|>\upsilon\end{cases}$ (13)

where the threshold $\upsilon$ is a positive real number. By replacing the
quantile error function $x$ in the absolute value function 10 with the Huber
function 13, we get the new objective function and its first and second
derivatives as follows:

$\begin{aligned}
\rho_{\tau}=\begin{cases}(\tau-1)(|t|-\frac{\upsilon}{2})&t<-\upsilon\\\
(\tau-1)(\frac{t^{2}}{2\upsilon})&-\upsilon\leq t<0\\\
\tau(\frac{t^{2}}{2\upsilon})&0\leq t<\upsilon\\\
\tau(|t|-\frac{\upsilon}{2})&t>\upsilon\end{cases}\end{aligned}\quad\text{,}\quad\begin{aligned}
\nabla_{\rho}=\begin{cases}(1-\tau)&t<-\upsilon\\\
(\tau-1)(\frac{t}{\upsilon})&-\upsilon\leq t<0\\\
\tau(\frac{t}{\upsilon})&0\leq t<\upsilon\\\
\tau&t>\upsilon\end{cases}\end{aligned}\quad\text{, and}\quad\begin{aligned}
\nabla_{\rho}^{2}=\begin{cases}0&t<-\upsilon\\\
\frac{\tau-1}{\upsilon}&-\upsilon\leq t<0\\\ \frac{\tau}{\upsilon}&0\leq
t<\upsilon\\\ 0&t>\upsilon\end{cases}\end{aligned}$ (14)

Our proposed QXGBoost method uses equation 14 as the objective function for
XGBoost to construct PIs. Equation 14 proves that the proposed objective
function is differentiable in the defined space, and the second derivative
exists and is not 0 everywhere. Figure 2 compares the original absolute value
function and the Huber smoothed absolute value functions with different
thresholds when the quantile $\tau=0.95$. As we can see, the Huber smoothed
absolute value function $h_{\upsilon}(t)$ is continuously differentiable and
provides a smooth transition between absolute and square errors around the
origin [33]. When the threshold ($\upsilon$) is smaller than 1, the Huber
smoothed absolute value functions approach the absolute value function
infinitely. As the threshold $\tau$ gets larger, the Huber smoothed absolute
value functions depart further away from the original absolute value function.
Thus, in practice, we prefer to use smaller threshold $\tau$ values.

Figure 2: The difference between the absolute value function and Huber
approximation function with different thresholds when the quantile
$\tau=0.95$.

In summary, inspired by the idea of quantile neural networks, GBM, LightGBM,
and random forest, we came up with the idea of combining quantile regression
with XGBoost for the latter to do uncertainty quantification. For the quantile
regression function to be a qualified objective function (have first and
second derivatives everywhere in the defined sample space and not all 0), we
introduced the Huber smoothed absolute value function, which is a
differentiable approximation to the absolute value function that provides a
smooth transition between absolute and squared errors around the origin.
Accordingly, two sets of parameters need to be tuned for QXGBoost. One set of
parameters is inhabited in the XGBoost model, the most important parameters of
which include the number of estimators, the maximum tree depth, and the
learning rate. The other set comes with the threshold parameter in the Huber
function ($\upsilon$) inside the smoothed quantile regression function. Both
sets of parameters need to be tuned to achieve the desired model performance.

##  5 Evaluation of prediction intervals

The evaluation of traditional point predictions usually adopts statistical
methods, such as root mean square error, mean absolute error, and mean
absolute percentage error. However, the evaluation of PIs is not as
straightforward and no established approach exists. In practice, the most
commonly used criteria for quantifying PIs are the prediction interval
coverage probability (PICP) and prediction interval average width (PIAW) [4,
5, 34, 35].

PICP measures the percentage of the target data points that fall in the bounds
of the lower and upper limits of the PIs inclusive. The mathematical
expression of PICP is represented in equation 15

$PICP=\frac{1}{n}\sum_{i=1}^{n}c_{i}$ (15)

where

$c_{i}=\begin{cases}1&t_{i}\in[\underline{y}_{i},\bar{y}_{i}]\\\
0&otherwise\end{cases}$ (16)

Here, $n$ represents the number of samples. $t_{i}$ is the true $i^{th}$
target value. $\underline{y}_{i}$ and $\bar{y}_{i}$ are the lower and upper
bound of the $i^{th}$ sample, individually. Based on equation 15, PICP depends
on the width of PIs. Usually wider width of PIs means a larger PICP. However,
a very large (for example, 100%) is meaningless as it conveys no information
about the target values.

PIAW calculates the average width of PIs. Given that different datasets
usually have different sample ranges, PIAW needs to be normalized so that it
can be used to compare across different datasets. The normalized PIAW (PINAW)
is given as

$PINAW=\frac{1}{R\times n}\sum_{j=1}^{n}(\bar{y}_{i}-\underline{y}_{i})$ (17)

where $R$ is the range of the observed data. PICP and PINAW are positively
correlated and typically, a large PICP with small PINAW is preferred. Thus,
PICP and PINAW should be used together to co-evaluate the quality of PIs. In
practice, we could control one criterion and compare the other criteria among
different models. For example, to compare PIs from different models, we could
try to set the PIAPs all the same and compare PINAWs. The smaller PINAW is,
the better the model is. However, most of the time, Both PICP and PINAW vary
among models. To solve the dilemma of PICP and PINAW going in opposite
directions, a new criterion called the coverage width-based criterion (CWC),
which combines both PICP and PINAW is proposed [6]

$CWC=\begin{cases}PINAW(1+e^{\eta(\tau-PICP)})&PICP<\tau\\\
PINAW&PICP\geq\tau\end{cases}$ (18)

where $\eta=50$ is a hyperparameter and $\tau$ is the desired quantile [4].
CWC is used in practice to assess the quality of PIs across different models.
All three methods were adopted to quantify the PIs in the following
applications.

##  6 Application to simulated and real-world data

In this section, we applied QXGBoost on simulated data and real-world traffic
noise data. To evaluate the performance of QXGBoost, we compared the quality
of PIs from QXGBoost with two other uncertainty quantification methods,
quantile GBM, and quantile LightGBM. For each method, 90% PIs were obtained
using two models, one predicting the upper bound of PIs with $\tau=0.95$ and
one for the lower bound with $\tau=0.05$. Data were randomly split and all
models were trained on 75% of the data and tested on the rest 25%.

###  6.1 Simulated data

The simulation data that we chose was adapted from the Scikit-learn library
example which shows how quantile regression can be used to create PIs for GBM
and it was generated using the equation 19

$y=1.5\times x\times sin(x)+\varepsilon$ (19)

where $\varepsilon$ has a normal distribution with mean 0 and standard
deviation that was randomly drawn between 1.5 and 2.5 so that we have
heteroscedastic noise for the simulated data. The relationship between $x$ and
$y$ is not linear.

The hyperparameters tuning is shown in table 1. Quantile GBM and LightGBM
shared the same tree-based parameters tuning. QXGBoost used a larger number of
trees and maximum tree depth. All three methods used the same learning rate
(0.05) for the simulated data. The extra threshold parameter $\tau$ in
QXGBoost was set to 2, meaning that the smoothed loss function approximates
well to the original absolute value function (figure 2).

Table 1: Hyperparameters tuning for the simulated data. Method | learning rate | Number of estimators | Max depth | Huber threshold
---|---|---|---|---
Quantile GBM | 0.05 | 200 | 2 | –
Quantile LightGBM | 0.05 | 200 | 2 | –
QXGBoost | 0.05 | 300 | 3 | 2

Table 2 shows that even though QXGBoost has the smallest training PICP
(0.891), it has the largest testing PICP (0.892), meaning that QXGBoost does
not overfit and can be generalized to other datasets better than quantile
LightGBM and quantile GBM. The higher PICP of QXGBoost also comes with a
slightly larger PINAW (0.777), which is unsurprising. Overall, QXGBoost
outperformed quantile GBM and LightGBM with the smallest CWC (1.937 VS 3.704
and 3.682). Figure 3 shows that all three models have roughly the same PIs
when $x<7$. QXGBoost has smaller values for the lower quantile and higher
values for the upper quantile when $7<x<9$ and $9<x<10$, respectively.

However, this is a very simple example with only one independent variable,
$x$. For XGBoost, this means that there is only one variable used to split and
build all the trees. The strong suit of ML methods is their capability to
untangle the complex non-linear relationship among the independent variables.
Thus, a more complicated example is illustrated using a real-life
environmental dataset.

Table 2: Prediction intervals results for the traffic noise data. Method | Training PICP | Testing PICP | testing PINAW | Testing CWC
---|---|---|---|---
Quantile GBM | 0.904 | 0.872 | 0.733 | 3.704
Quantile LightGBM | 0.901 | 0.872 | 0.728 | 3.682
QXGBoost | 0.891 | 0.892 | 0.777 | 1.937
Figure 3: Comparison of prediction intervals for the simulated data using
three methods. The blue dots is the observations. The red line is the point
prediction of the traffic noise. The two black lines represent the upper and
lower bounds of PIs and the blue areas in between are the PIs.

###  6.2 Traffic noise data

Traffic noise data came from Yin et al. (2020) [36]. In brief, traffic noise
was measured while walking along 16 predefined routes in Long Beach, CA with
two hand-held Brüel & Kjær 2250 sound analyzers, which are high-end sound
meters that allow the user to measure A-weighted, equivalent noise (LAeq) at
intervals as short as one second. The sound meters were coupled with Strava to
get real-time GPS location data. After data processing, a total of 6,647 geo-
referenced observations were obtained. Based on their latitude and longitude,
the 44 predictor variables, including meteorology, traffic, land use, and road
features, were extracted to build and train the models.

Before performing the uncertainty quantification analysis, we first run the
point prediction analysis of traffic noise using GBM, LightGBM, and XGBoost.
The parameters of each algorithm were tuned by a 5-fold CV random grid search
process that included measured average LAeq from all of the observations
together. After tuning, each model was trained on 70% of the data and
predicted on the remaining 30%; the prediction results are shown in Figure 4.
GBM and LightGBM have very similar prediction performance with the same R2 =
0.96 and RMSE = 2.04. XGBoost has a better model performance with a slightly
higher R2 = 0.97 and a lower RMSE = 1.89.

Figure 4: Point prediction for traffic noise using GBM, LightGBM, and XGBoost.
Blue points represent each observed and predicted sample; solid red lines
represent the regression line; and dashed black lines represent the 1:1
correspondence between measured and predicted values.

The tuned hyperparameters for Quantile GBM, Quantile Light GBM, and QXGBoost
are shown in table 3. Quantile GBM and LightGBM have the same tree-based
tuning parameters. By comparison, QXGBoost used more trees (300 vs 200) to
build the models but also a higher learning rate (0.05 vs 0.03). The Huber
threshold for the traffic noise data was tuned at 0.07, which gave us a very
close approximation to the real quantile regression function according to
figure 2.

Table 3: Hyperparameters tuning for the traffic noise data. Method | learning rate | Number of estimators | Max depth | Huber threshold
---|---|---|---|---
Quantile GBM | 0.03 | 200 | 3 | –
Quantile LightGBM | 0.03 | 200 | 3 | –
QXGBoost | 0.05 | 300 | 3 | 0.07

Table 4 summarizes the training PICP and the testing PICP, PINAW, and CWC from
all three models. Among them, quantile LightGBM has the largest training PICP
(0.907) and the smallest testing PICP (0.885), indicating that the LightGBM
might have an overfitting problem. Even though quantile GBM has the largest
testing PICP (0.890), it also has the biggest PI width (0.331). By comparison,
QXGBoost has a reasonable testing PICP (0.889) with the smallest PINAW
(0.319). Overall, QXGBoost has the smallest testing CWC (4.580), followed by
quantile LightGBM and quantile GBM (4.610 and 4.751, respectively).

Table 4: Prediction intervals results for the traffic noise data. Method | Training PICP | Testing PICP | testing PINAW | Testing CWC
---|---|---|---|---
Quantile GBM | 0.905 | 0.890 | 0.331 | 4.751
Quantile LightGBM | 0.907 | 0.885 | 0.320 | 4.610
QXGBoost | 0.902 | 0.889 | 0.319 | 4.580

The PIs generated for all three models are shown in Figure 5. The top row of
Figure 5 shows that all three models have wider PIs for lower traffic noise
levels, and LightGBM tends to over-predict traffic noise compared to GBM and
XGBoost. In terms of the mid-level values of traffic noise, the PIs of
QXGBoost cover most of the extremely high and low point predictions whereas
quantile GBM and LightGBM do not. This is why the QXGBoost PIs look wider than
the Quantile GBM and LightGBM in the middle noise level. One thing to notice
for QXGBoost is that it has several very wide PIs at the high end of traffic
noise for some observations, whereas the Quantile GBM and LightGBM have
relatively consistent lengths of PIs, indicating that quantile GBM and
LightGBM might be more stable than the QXGBoost in predicting the high traffic
noise levels. Both the first and second rows of the figure 5 show that PIs
vary in length, with some being much shorter than others. The relatively wider
PIs at both ends of the traffic noise levels are due to a smaller number of
observations. Figure 4 shows that most of the traffic noise values were
distributed between 50 and 75 dB and the ones that are below 50 dB or above 75
dB consisted only 2.9% and 2.6% of all the observations, respectively. Thus,
it is not unexpected that with less training data, models will generate PIs of
high uncertainty with a wider length.

Figure 5: Prediction intervals (turquoise) for the traffic noise data from
quantile gradient boosting (left), quantile light gradient boosting (center),
and quantile extreme gradient boosting (right). The red dots are traffic noise
point predictions. The figures on the top row show the prediction intervals
for the test set; the figures on the bottom row show the prediction intervals
for 20 randomly selected points from the test set.

For better visualization, Figure 6 shows all traffic noise observations
ordered according to the length of the corresponding PIs. It can be easily
seen that the overall width of PIs of the quantile GBM is wider than that of
the quantile LightGBM and QXGBoost. QGBoost has a slightly tighter width than
the quantile LightGBM. All three models performed relatively well at the
middle-level traffic noise and covered almost 90% of the observed traffic
noise in the testing dataset. Figure 7 shows the PIs from all three models
together for better comparison. Quantile GBM and LightGBM have a wider length
of PIs than QXGBoost for most of the observations (around 1300 out of 1649
observations). For the rest of the observations, the PIs width of QXGBoost is
wider as was reflected at the high index of figure 7, indicating that the
widths of the PIs from QXGBoost vary more than those of quantile GBM and
quantile LightGBM.

Figure 6: Ordered and centered PIs for the traffic noise data from all three
models. The red dots represent the observations that are covered by PIs. The
green dots are the observations that are not covered by PIs. The blue areas
show the ordered and centered PIs. Figure 7: The figure on the second row is
the ordered and centered PIs from all three models on the same plot.

##  7 Conclusion and Discussion

We developed QXGBoost, a new method that employs quantile regression as the
customized objective function for XGBoost for uncertainty analysis. In order
to have an objective function that can be differentiable everywhere and have
second derivatives, we adapted the concept developed in previous studies for
neural networks [31, 33] by introducing the Huber norm function to provide a
smooth approximation to the true quantile regression. To our knowledge, this
is the first study that combines quantile regression, Huber norm function, and
XGBoost to do uncertainty analysis for XGBoost. The results from our
experiments suggest that QXGBoost produces competitive or better quality PIs
compared to quantile GBM and LightGBM.

The primary goal of this study was to develop an easy-to-understand
implementation of uncertainty quantification for XGBoost. Even though popular
uncertainty quantification methods such as Bayesian or Delta methods for the
neural networks have been studied for years [6], there are some disadvantages
to using them. First, there is a learning curve for non-experts to understand
Bayesian or Delta methods. Second, their integration with XGBoost is very
challenging due to the fact that tree-based ML methods are fundamentally
different from neural network models. GBM was able to use the quantile
regression function directly because it does not require second-order
approximation for the objective function. LightGBM is similar to XGBoost in
terms of its objective optimization method. However, the success of quantile
LightGBM requires diving into the source code (c++) and exporting the tree
leaf index to the objective function, which would be difficult for people who
are not familiar with the source code and programming language. QXGBoost
solved the problem simply by defining an optimized objective function. Even
though QXGBoost has one additional parameter to tune, it provides a simple and
quick solution for performing uncertainty analysis for one of the most
powerful ML methods.

The additional threshold hyperparameter $\tau$ plays an important role during
model fitting. Based on our experiment, $\tau$ controls both the PICP and
PINAW of the PIs. Small $\tau$ (less than 1) usually leads to large PICP and
PINAW and large $\tau$s typically produce small PICP and PINAW, meaning more
constraints on the PIs. In practice, we prefer large PIs coverage with a
narrower width. Thus, both the threshold and XGBoost parameters need to be
tuned and balanced to achieve better model performance. However, as is shown
in figure 2, as $\tau$ gets larger (more than 2), the Huber smoothed
approximation function departs further away from the original quantile
regression function, and the poor approximation might lead to biased PIs. In
practice, we want to control the threshold parameter $\tau$ under 2.

There are limitations to our method. First, as we have seen from figure 5, for
extremely low and high traffic noise values, when there were few observations
for training, QXGBoost could not give informative PIs because less data led to
greater uncertainty in the predictions. This is true for other ML models
including GBM and LightGBM, which is why all three models have better PI
quality for mid-level than at low and high ends of traffic noise. We generally
see overestimation for low levels of traffic noise and underestimation for
high-level traffic noise. Second, we usually could not achieve the "hard"
coverage as the $\tau$ specifies for the proposed QXGBoost method. For
example, we have to uniformly pad the quantile intervals by 3% and 1% to get a
90% coverage using QXGBoost for the simulated and traffic noise data,
respectively. This limitation is not unique to our method. As is shown in
other studies [11, 37], additional padding was also needed for quantile random
forest to achieve desired quantile coverage. Given these limitations, future
efforts should be focused on optimizing the objective function of the XGBoost
so that the “hard” coverage can be achieved without additional padding.

Data uncertainty exists in every stage of research from data collection,
modeling, analysis, visualization, and thus drawing conclusions. XGBoost has
been widely used in all research areas from image classification to disease
diagnosis. The popularity of the XGBoost algorithm motivates the need for its
ability to do uncertainty quantification to better understand and capture
uncertainty during the modeling stage. Ignoring uncertainty can lead to biased
results and thus hinder the improvement of health policies and any other
science-based decision making. As the first study proposing an uncertainty
method for XGBoost, we fill in the knowledge gap of modelling uncertainty
using one of the state-of-art ML algorithms. One important innovation of
QXGBoost is its easy implementation and customization. The success of QXGBoost
provides potential applications for epidemiological studies, both in modeling
exposure uncertainty and assessing health effects. With the increased
availability and ubiquity of spatial and health data, uncertainty
quantification methods such as QXGBoost can help us better understand the
causal effects relationship between exposures and outcomes.

## Acknowledgments

The research described in this article was conducted under contract to the
Health Effects Institute (HEI), an organization jointly funded by the United
States Environmental Protection Agency (EPA) (Assistance Award No.
CR-83590201) and certain motor vehicle and engine manufacturers. The contents
of this article do not necessarily reflect the views of HEI, or its sponsors,
nor do they necessarily reflect the views and policies of the EPA or motor
vehicle and engine manufacturers.

## References

  * [1] Christopher Irrgang, Niklas Boers, Maike Sonnewald, Elizabeth A. Barnes, Christopher Kadow, Joanna Staneva, and Jan Saynisch-Wagner.  Towards neural Earth system modelling by integrating artificial intelligence in Earth system science.  Nature Machine Intelligence, 3(8):667–674, 2021.
  * [2] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger.  On calibration of modern neural networks.  In International Conference on Machine Learning, pages 1321–1330. PMLR, 2017.
  * [3] National Academies of Sciences, Engineering, and Medicine.  Machine Learning and Artificial Intelligence to Advance Earth System Science: Opportunities and Challenges: Proceedings of a Workshop.  The National Academies Press, Washington, DC, 2022.
  * [4] HM Dipu Kabir, Abbas Khosravi, Mohammad Anwar Hosen, and Saeid Nahavandi.  Neural network-based uncertainty quantification: A survey of methodologies and applications.  IEEE access, 6:36218–36234, 2018.
  * [5] Vahid Nourani, Nardin Jabbarian Paknezhad, and Hitoshi Tanaka.  Prediction interval estimation methods for artificial neural network (ann)-based modeling of the hydro-climatic processes, a review.  Sustainability, 13(4):1633, 2021.
  * [6] Abbas Khosravi, Saeid Nahavandi, Doug Creighton, and Amir F Atiya.  Comprehensive review of neural network-based prediction intervals and new advances.  IEEE Transactions on neural networks, 22(9):1341–1356, 2011.
  * [7] Enrico Zio.  A study of the bootstrap method for estimating the accuracy of artificial neural networks in predicting nuclear transient processes.  IEEE Transactions on Nuclear Science, 53(3):1460–1478, 2006.
  * [8] Rachid Errouissi, Julian Cardenas-Barrera, Julian Meng, Eduardo Castillo-Guerra, Xun Gong, and Liuchen Chang.  Bootstrap prediction interval estimation for wind speed forecasting.  In 2015 IEEE Energy Conversion Congress and Exposition (ECCE), pages 1919–1924. IEEE, 2015.
  * [9] Arni Magnusson, André E Punt, and Ray Hilborn.  Measuring uncertainty in fisheries stock assessment: the delta method, bootstrap, and mcmc.  Fish and Fisheries, 14(3):325–342, 2013.
  * [10] Liang-Cheng Zhang and AC Worthington.  Evaluating the accuracy of scope economies: comparisons among delta method, bootstrap, and bayesian approach.  In Australian Conference of Economists PhD Colloquium. Retrieved from http://www98. griffith. edu. au/dspace/handle/10072/69409, 2015.
  * [11] Quantile regression forests and bootstrap.  https://saattrupdan.github.io/2020-04-05-quantile-regression-forests/.  Accessed: 2023-02-22.
  * [12] James W Taylor.  A quantile regression neural network approach to estimating the conditional density of multiperiod returns.  Journal of Forecasting, 19(4):299–311, 2000.
  * [13] Yaoyao He and Haiyan Li.  Probability density forecasting of wind power using quantile regression neural network and kernel density estimation.  Energy conversion and management, 164:374–384, 2018.
  * [14] Yaoyao He, Yang Qin, Shuo Wang, Xu Wang, and Chao Wang.  Electricity consumption probability density forecasting method based on lasso-quantile regression neural network.  Applied energy, 233:565–575, 2019.
  * [15] Alex J Cannon.  Non-crossing nonlinear regression quantiles by monotone composite quantile regression neural network, with application to rainfall extremes.  Stochastic environmental research and risk assessment, 32(11):3207–3225, 2018.
  * [16] Wen-Tsao Pan, Chiung-En Huang, and Chiung-Lin Chiu.  Study on the performance evaluation of online teaching using the quantile regression analysis and artificial neural network.  The Journal of Supercomputing, 72(3):789–803, 2016.
  * [17] Qifa Xu, Kai Deng, Cuixia Jiang, Fang Sun, and Xue Huang.  Composite quantile regression neural network with applications.  Expert Systems with Applications, 76:129–139, 2017.
  * [18] Nicolai Meinshausen and Greg Ridgeway.  Quantile regression forests.  Journal of Machine Learning Research, 7(6), 2006.
  * [19] Happy Aprillia, Hong-Tzer Yang, and Chao-Ming Huang.  Statistical load forecasting using optimal quantile regression random forest and risk assessment index.  IEEE Transactions on Smart Grid, 12(2):1467–1480, 2020.
  * [20] Kévin Vaysse and Philippe Lagacherie.  Using quantile regression forest to estimate uncertainty of digital soil mapping products.  Geoderma, 291:55–64, 2017.
  * [21] T Francke, JA López-Tarazón, and B Schröder.  Estimation of suspended sediment concentration and yield using linear models, random forests and quantile regression forests.  Hydrological Processes: An International Journal, 22(25):4892–4904, 2008.
  * [22] Alexey Natekin and Alois Knoll.  Gradient boosting machines, a tutorial.  Frontiers in neurorobotics, 7:21, 2013.
  * [23] Mark Landry, Thomas P Erlinger, David Patschke, and Craig Varrichio.  Probabilistic gradient boosting machines for gefcom2014 wind forecasting.  International Journal of Forecasting, 32(3):1061–1066, 2016.
  * [24] Hadrien Verbois, Andrivo Rusydi, and Alexandre Thiery.  Probabilistic forecasting of day-ahead solar irradiance using quantile gradient boosting.  Solar Energy, 173:313–327, 2018.
  * [25] Jing Wei, Zhanqing Li, Rachel T Pinker, Jun Wang, Lin Sun, Wenhao Xue, Runze Li, and Maureen Cribb.  Himawari-8-derived diurnal variations in ground-level pm 2.5 pollution across china using the fast space-time light gradient boosting machine (lightgbm).  Atmospheric Chemistry and Physics, 21(10):7863–7880, 2021.
  * [26] Antoine Guitton and William W Symes.  Robust inversion of seismic data using the huber norm.  Geophysics, 68(4):1310–1319, 2003.
  * [27] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785–794, 2016.
  * [28] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Additive logistic regression: a statistical view of boosting (with discussion and a rejoinder by the authors).  The annals of statistics, 28(2):337–407, 2000.
  * [29] Roger Koenker and Gilbert Bassett Jr.  Regression quantiles.  Econometrica: journal of the Econometric Society, pages 33–50, 1978\.
  * [30] Keming Yu, Zudi Lu, and Julian Stander.  Quantile regression: applications and current research areas.  Journal of the Royal Statistical Society: Series D (The Statistician), 52(3):331–350, 2003.
  * [31] Colin Chen.  A finite smoothing algorithm for quantile regression.  Journal of Computational and Graphical Statistics, 16(1):136–164, 2007.
  * [32] Peter J Huber.  Robust regression: asymptotics, conjectures and monte carlo.  The annals of statistics, pages 799–821, 1973.
  * [33] Alex J Cannon.  Quantile regression neural networks: Implementation in r and application to precipitation downscaling.  Computers & geosciences, 37(9):1277–1284, 2011.
  * [34] Tim Pearce, Alexandra Brintrup, Mohamed Zaki, and Andy Neely.  High-quality prediction intervals for deep learning: A distribution-free, ensembled approach.  In International Conference on Machine Learning, pages 4075–4084. PMLR, 2018.
  * [35] Yuandu Lai, Yucheng Shi, Yahong Han, Yunfeng Shao, Meiyu Qi, and Bingshuai Li.  Exploring uncertainty in deep learning for construction of prediction intervals.  arXiv preprint arXiv:2104.12953, 2021.
  * [36] Xiaozhe Yin, Masoud Fallah-Shorshani, Rob McConnell, Scott Fruin, and Meredith Franklin.  Predicting fine spatial scale traffic noise using mobile measurements and machine learning.  Environmental Science & Technology, 54(20):12860–12869, 2020.
  * [37] Quantile regression forests.  https://scikit-garden.github.io/examples/QuantileRegressionForests/.  Accessed: 2023-02-22.