#  Conformalized Fairness via Quantile Regression

Meichen Liu1, Lei Ding1, Dengdeng Yu2, Wulong Liu3,
Linglong Kong1, Bei Jiang1∗
1Department of Mathematical and Statistical Sciences, University of Alberta
2Department of Mathematics, University of Texas at Arlington
3 Huawei Noah’s Ark Lab Canada
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Co-Corresponding Authors

###### Abstract

Algorithmic fairness has received increased attention in socially sensitive
domains. While rich literature on mean fairness has been established, research
on quantile fairness remains sparse but vital. To fulfill great needs and
advocate the significance of quantile fairness, we propose a novel framework
to learn a real-valued quantile function under the fairness requirement of
_Demographic Parity_ with respect to sensitive attributes, such as race or
gender, and thereby derive a reliable _fair_ prediction interval. Using
optimal transport and functional synchronization techniques, we establish
theoretical guarantees of distribution-free coverage and exact fairness for
the induced prediction interval constructed by fair quantiles. A hands-on
pipeline is provided to incorporate flexible quantile regressions with an
efficient fairness adjustment post-processing algorithm. We demonstrate the
superior empirical performance of this approach on several benchmark datasets.
Our results show the model’s ability to uncover the mechanism underlying the
fairness-accuracy trade-off in a wide range of societal and medical
applications.

##  1 Introduction

We are increasingly leaning on machine learning systems to tackle human
problems. A primary objective is to develop intelligent algorithms that can
automatically produce accurate decisions which also enjoy equitable
properties, as unintended social bias has been identified as a rising concern
in various fields [15, 20, 16].

As a means of providing quantitative measures of fairness, a number of metrics
have been proposed. These metrics can be categorized into three broad
categories: group fairness [4], individual fairness [28], and causality-based
fairness [34]. In contrast to causality-based fairness that requires domain
knowledge to develop a fair causal structure and individual fairness that
seeks equality only between similar individuals, group fairness does not
require any prior knowledge and seeks equality for groups as a whole [7].
Among the metrics defined for group fairness such as equalized odds [10, 32]
and predictive rate parity [11], demographic parity (DP) is generic since it
does not allow prediction results in aggregate to depend on sensitive
attributes [2, 24, 13, 46]. In particular, an algorithm is said to satisfy DP
if its prediction is independent of any given sensitive attribute [2].

There have been a number of studies on algorithmic fairness concerning DP [2,
12, 13, 24, 35, 46]. In the context regression analysis, much attention have
been paid on conditional mean inferences [2, 12, 13, 35], few are concerned
with conditional quantiles [44, 46]. As real-world data often exhibit
heterogeneity, contain extreme outliers, or do not meet satisfactory
distributional assumptions, like Gaussianity, a fairness discussion on
conditional quantiles may be more rational and essential since they are able
to provide a more complete understanding of the dependence structure between
response and explanatory variables [46], as well as better accommodate
asymmetry and extreme tail behavior [45]. It should also be noted that bias or
unfairness that arises in mean regression may also be propagated through
quantile regression, therefore it must be properly dealt with separately: a
graphic demonstration can be found in Figure 1. More intuitively, we may take
an example from a Spanish labor market study [19, 20]. The study found that in
Spain, also in line with other countries, the mean wage gap between men and
women is quite substantial: on average, women earn around 70 percent of what
men earn. While wage gaps are not uniform across all pay scales, they are
greater at higher quantiles than at lower quantiles. As biases and disparities
at different quantiles tend to be overshadowed by the mean behavior of the
entire population, we propose a novel framework to seek fair predictions at
different quantiles. It uses optimal transport techniques [3, 13] by
transforming _bias-affected_ distributions into an _only-fair_ Wasserstein-2
barycenter through a kernel-based functional synchronization method [9, 52],
in order to provide fair quantile estimators.

Figure 1: An illustration of quantile fairness: for a skewed and
heteroscedastic quantile estimation $\\{q_{\alpha,i}\\}_{i=1}^{N}$ affected by
the sensitive attribute $S\in\\{0,1\\}$, for example, the higher quantile of
the salary distribution, the optimal fair quantile prediction
$Q_{g_{\alpha}}(t),t\in(0,1)$ is derived through a convex combination of the
conditional quantile functions of $Q_{q_{\alpha}|S=0}$ and
$Q_{q_{\alpha}|S=1}$.

Since quantile fairness poses a number of theoretical challenges, no previous
literature has been able to provide any inference results such as prediction
intervals concerning quantile fairness. It is imperative to keep in mind that
fairness is only one of two legs of the primary goal of modern machine
learning algorithms, the other being accuracy. Building a reliable prediction
with valid confidence is a significant challenge that is encountered by many
machine learning algorithms [51]. Towards this end, we propose the
conformalized fair quantile prediction (CFQP) inspired by the works of Romano
et al. [37, 38]. Our analysis demonstrates, both mathematically and
experimentally, that CFQR provides finite sample, distribution-free validity,
DP fairness for different quantiles, and precise control of the miscoverage
rate, regardless of the underlying quantile algorithm.

Contributions and Outlines. In this paper, we propose a new quantile based
method with valid inference that enhances both accuracy and fairness while
maintaining a balance between the two. It is a novel framework that allows an
exact control of prediction miscoverages while ensuring quantile fairness
simultaneously. The main contributions are summarized as follows:

  1. i.

We successfully transform the problem of searching quantiles under DP fairness
to the construction of multi-marginal Wasserstein-2 barycenters via the
optimal transport theory [3, 13, 21]. We incorporate a novel kernel smoothing
step into the preceding method, which is particularly advantageous for
subgroups whose sample sizes are too small to obtain reliable quantile
function estimations.

  2. ii.

In Section 4, we propose a conformalized fair quantile regression prediction
interval (CFQP) inspired by the works of Romano et al. [37, 38]. It is
mathematically proved to achieve a distribution-free validity, demographic
parity on different quantiles, and an exact control of miscoverage rates,
regardless of the quantile algorithm used. The theoretical validity of
prediction interval constructed by CFQP and exact DP of the fair quantile
estimators are given in Section 5 and the supplement.

  3. iii.

The experimental results presented in Section 6 include a numerical comparison
of the proposed CFQP and fair quantile estimation with both state-of-the-art
conformal and fairness-oriented methods. By reducing the discriminatory bias
dramatically, our method outperforms the state-of-the-art methods while
maintaining reasonable short interval lengths.

Related works. Existing approaches for building a fair mean regression broadly
fall into three classes: pre-processing, in-processing and post-processing. In
particular, preprocessing methods focus on transforming the data to remove any
unwanted bias [6, 34, 50]; in-processing methods aim to build in fairness
constraints into the training step [2, 5, 28, 32]; post-processing methods
target to modify the trained predictor [12, 13, 31]. As few previous works
have focused on the quantile fairness of and fair prediction interval, the
most related are Yang et al. [46], where a different fairness measure was
used. While Agarwal et al. [2] mentioned that their reduction-based approach
can be adapted into quantile regression, Williamson and Menon [44] brought
forward a novel conditional variance at risk fairness measure aiming to
control the largest subgroup risk. For interval fairness measure, the approach
by Romano et al. [37] achieved equalized coverage among groups without
fairness on interval endpoints. Methodologically, integrating algorithmic
fairness with Wasserstein distance based barycenter problem has been studied
in [3, 12, 13, 21, 26]. Both in-processing [2, 26] and post-processing [12,
13] methods were proposed to solve classification and mean regression
problems. As a post-processing method, our work is distinct from above-
mentioned methods by constructing the DP-fairness for each population
quantile, and generating a fair prediction interval accordingly.

Notations. We denote by $[K]$ the set $\\{1,\ldots,K\\}$ for arbitrary integer
$K$. $|\mathcal{S}|$ represents the cardinality for a finite set
$\mathcal{S}$. ${E}$ and ${P}$ represent the expectation and probability and
$\mathbbm{1}\\{\cdot\\}$ is the indicator function. Let
$\\{Z_{n}\\}_{n=1}^{\infty}$ be a sequence of random variables, and
$\\{k_{n}\\}_{n=1}^{\infty}$ be a sequence of positive numbers, we say that
$Z_{n}=O_{p}(k_{n})$, if
$\lim_{T\rightarrow\infty}\limsup_{n\rightarrow\infty}{P}(|Z_{n}|>Tk_{n})=0$,
then $Z_{n}/k_{n}=O_{p}(1)$. To denote the equality in distribution of two
random variables $A$ and $B$, we write $A\overset{d}{=}B$.

##  2 Problem statement

Consider the regression problem where a “sensitive characteristic” $S$ is
available, which by the U.S. law [21, 37] can be enumerated as sex, race, age,
disability, etc. We observe the triplets
$\left({X}_{1},S_{1},Y_{1}\right),\dots,\left({X}_{n},S_{n},Y_{n}\right)$,
denote $({X}_{i},S_{i},Y_{i})$ by $Z_{i}$, $i=1,\dots,n$ and $Z_{i}$ is a
random variable in $\mathbb{R}^{p}\times[K]\times\mathbb{R}$. The aim is to
predict the unknown value of $Y_{n+1}$ at a test point $X_{n+1},S_{n+1}$. Let
${P}$ be the joint distribution of $Z$, we assume that all the samples
$\left\\{Z_{i}\right\\}_{i=1}^{n+1}$ are drawn exchangeable, where i.i.d. is a
special case.

Our goal is to construct a marginal distribution-free prediction band
$C\left(X_{n+1},S_{n+1}\right)\subseteq\mathbb{R}$ that is likely to cover the
unknown response $Y_{n+1}$ with finite-sample (nonasymptotic) validity.
Formally, given a desired miscoverage rate $\alpha$, the predicted interval
satisfies

${P}\left\\{Y_{n+1}\in C\left(X_{n+1},S_{n+1}\right)\right\\}\geq 1-\alpha$
(1)

for any joint distribution ${P}$ and any sample size $n$, while the left and
right endpoint of $C\left(X_{n+1},S_{n+1}\right)$ satisfies the fairness
constraint of Demographic Parity concerning the sensitive variable $S$.

Demographic Parity. We introduce the quantitative definition of DP in fair
regression and connect the DP-fairness with a quantile regressor $q_{\alpha}$.
The result that $q_{\alpha}$ can be projected to the fair counterparts using
optimal transport will be invoked later.

Given a fixed quantile level $\alpha$ (it may refer to
${\alpha_{\mathrm{lo}}}$ or ${\alpha_{\mathrm{hi}}}$ indicating the upper and
lower quantile estimates for the prediction band $C(X_{n+1},S_{n+1})$). Let
$q_{\alpha}:\mathbb{R}^{p}\times[K]\rightarrow\mathbb{R}$ represent an
arbitrary conditional quantile predictor. Denote by $\nu_{q_{\alpha}\mid s}$
the distribution of $(q_{\alpha}(X,S)\mid S=s)$, the Cumulative Distribution
Function (CDF) of $\nu_{q_{\alpha}\mid s}$ is given by

$F_{\nu_{q_{\alpha}\mid s}}(t)={P}(q_{\alpha}(X,S)\leq t\mid S=s).$ (2)

The quantile function $Q_{\nu_{q_{\alpha}\mid s}}=F^{-1}_{\nu_{q_{\alpha}\mid
s}}:[0,1]\rightarrow\mathbb{R}$ ,namely, the generalized inverse of
$F_{\nu_{q_{\alpha}\mid s}}$, can thus be defined as for all levels
$t\in(0,1]$,

$Q_{\nu_{q_{\alpha}\mid s}}(t)=\inf\\{y\in\mathbb{R}:F_{\nu_{q_{\alpha}\mid
s}}(y)\geq t\\}\ \text{with}\ Q_{\nu_{q_{\alpha}\mid
s}}(0)=Q_{\nu_{q_{\alpha}\mid s}}(0+).$ (3)

To simplify the notations, we will write $F_{q_{\alpha}\mid s}$ and
$Q_{q_{\alpha}\mid s}$ instead of $F_{\nu_{q_{\alpha}\mid s}}$ and
$Q_{\nu_{q_{\alpha}\mid s}}$ respectively, for any prediction rule
$q_{\alpha}$.

In the following, we introduce the definition of Demographic Parity (DP),
which is most commonly used in the context of fairness research [2, 12, 13,
24, 32].

######  Definition 1 (Demographic Parity).

An arbitrary prediction $g:\mathbb{R}^{d}\times[K]\rightarrow\mathbb{R}$
satisfies demographic parity under a distribution ${P}$ over $(X,S,Y)$, if
$g(X,S)$ is statistically independent of the sensitive attribute $S$.
Formally, for every $s,s^{\prime}\in[K]$,

$\sup_{t\in\mathbb{R}}\left|{P}(g(X,S)\leq t\mid S=s)-{P}\left(g(X,S)\leq
t\mid S=s^{\prime}\right)\right|=0.$

Demographic Parity (DP) requires the predictions to be independent of the
sensitive attribute, and it demands the Kolmogorov-Smirnov distance [29] (the
difference between CDFs measured in the $l_{\infty}$ norm) between $\nu_{g\mid
s}$ and $\nu_{g\mid s^{\prime}}$ to vanish for all categories $s,s^{\prime}$.

##  3 Quantile Regression and Conformal Prediction

In this section, we recall the CQR approach for finite sample, distribution-
free prediction interval inference. Quantile regression was proposed by
Koenker and Bassett [27] to estimate the $\alpha$-th quantile of the
conditional distribution of $Y$ given $\tilde{X}:=(X,S)$ for some quantile
level $\alpha\in(0,1)$, since then it has become more pervasive with various
applications, such as providing prediction intervals, detecting outliers, or
perceiving the entire distribution [30, 22]. Denote the conditional cumulative
distribution of $Y$ given $\tilde{X}$ by
$F(y\mid\tilde{X}=\tilde{x}):={P}\\{Y\leq y\mid\tilde{X}=\tilde{x}\\}$. The
$\alpha$-th conditional quantile prediction is defined as
$q_{\alpha}(\tilde{x}):=\inf\\{y\in\mathbb{R}:F(y\mid\tilde{X}=\tilde{x})\geq\alpha\\}.$
Quantile regression can be cast as an optimization problem [30, 49, 33, 42,
48], by minimizing the expected check loss function
$E(\rho_{\alpha})=E[\rho_{\alpha}(y,q)|\tilde{X}=\tilde{x}]$, where

$\rho_{\alpha}(y,q_{\alpha}(\tilde{x}))=\left\\{\begin{aligned}
\alpha|y-q_{\alpha}(\tilde{x})|\quad&\text{if}\ y\geq
q_{\alpha}(\tilde{x}),\\\ (1-\alpha)|y-q_{\alpha}(\tilde{x})|\quad&\text{if}\
y<q_{\alpha}(\tilde{x}).\end{aligned}\right.$ (4)

Quantile regression offers a principled way of judging the reliability of
predictions by building a prediction interval for the new observation
$(\tilde{X}_{n+1},Y_{n+1})$. In contrast to asymptopia, Romano et al. [37, 38]
brought forward the conformalized quantile regression (CQR) by combining the
merits of robust quantile regression with conformal prediction, thus finite
sample validity in Eq. (1) is guaranteed. Inspired by the split conformal
method, a split CQR likewise starts with splitting the data into a proper
training set and a calibration set, indexed by
$\mathcal{I}_{1},\mathcal{I}_{2}$ respectively. Given any quantile regression
algorithm $\mathcal{Q}$, we then fit two conditional quantile functions
$\hat{q}_{\alpha_{\mathrm{lo}}}$ and $\hat{q}_{\alpha_{\mathrm{hi}}}$ on the
proper training set:
$\left\\{\hat{q}_{\alpha_{\mathrm{lo}}},\hat{q}_{\alpha_{\mathrm{hi}}}\right\\}\leftarrow\mathcal{Q}\left(\left\\{\left(\tilde{X}_{i},Y_{i}\right):i\in\mathcal{I}_{1}\right\\}\right).$
The conformity scores are calculated to quantify the error made by the plug-in
prediction interval
$\hat{C}(\tilde{x})=[\hat{q}_{\alpha_{\mathrm{lo}}}(\tilde{x}),\hat{q}_{\alpha_{\mathrm{hi}}}(\tilde{x})]$.
We evaluate the scores on the calibration set as
$E_{k}:=\max\left\\{\hat{q}_{\alpha_{\mathrm{lo}}}(\tilde{X}_{k})-Y_{k},Y_{k}-\hat{q}_{\alpha_{\mathrm{hi}}}(\tilde{X}_{k})\right\\}$
for each $k\in\mathcal{I}_{2}$, where both undercoverage and overcoverage of
the interval are taken into consideration [38]. Given a new input data
$\tilde{X}_{n+1}$, we construct the prediction interval for $Y_{n+1}$ as
$C\left(\tilde{X}_{n+1}\right)=\left[\hat{q}_{\alpha_{\mathrm{lo}}}\left(\tilde{X}_{n+1}\right)-Q_{1-\alpha}\left(E,\mathcal{I}_{2}\right),\hat{q}_{\alpha_{\mathrm{hi}}}\left(\tilde{X}_{n+1}\right)+Q_{1-\alpha}\left(E,\mathcal{I}_{2}\right)\right],$
where $Q_{1-\alpha}(E,\mathcal{I}_{2}):=(1-\alpha)(1+1/|\mathcal{I}_{2}|)$-th
empirical quantile of $\\{E_{k}:k\in\mathcal{I}_{2}\\}$ conformalizes the
plug-in prediction interval. Note that the constructed interval
$C(\tilde{X}_{n+1})$ could be highly influenced by the sensitive variable $S$.

##  4 Conformal fair quantile prediction (CFQP)

We formally describe our proposed conformal fair prediction (CFQP) framework
for constructing DP fairness constrained prediction intervals in this section.
A kernel smoothing quantile function is introduced during the functional
synchronization, which can improve the estimation when some subgroups are too
small to give reliable sample quantile function estimations.

######  Definition 2 (Wasserstein-2 distance).

Let $\mu$ and $\nu$ be two univariate probability measures with finite second
moments. The squared Wasserstein-2 distance between $\mu$ and $\nu$ is defined
as

$\mathcal{W}_{2}^{2}(\mu,\nu)=\inf\left\\{\int_{\mathbb{R}\times\mathbb{R}}|x-y|^{2}d\gamma(x,y),\gamma\in\Gamma_{\mu,\nu}\right\\}$

where $\Gamma_{\mu,\nu}$ is the set of probability measures (couplings) on
$\mathbb{R}\times\mathbb{R}$ having $\mu$ and $\nu$ as marginals.

######  Proposition 1 (Fair optimal prediction [13]).

Assume, for each $s\in[K]$, that the univariate measure $\nu_{q_{\alpha}\mid
s}$ has a density and let $p_{s}={P}(S=s)$. Then,

$\min_{g_{\alpha}\text{ is fair
}}{E}\left(q_{\alpha}(X,S)-g_{\alpha}(X,S)\right)^{2}=\min_{\nu}\sum_{s\in[K]}p_{s}\mathcal{W}_{2}^{2}\left(\nu_{q_{\alpha}\mid
s},\nu\right).$ (5)

Moreover, if $g_{\alpha}$ and $\nu$ solve the l.h.s. and the r.h.s. problems
respectively, then $\nu=\nu_{g_{\alpha}}$ and specifically,

$g_{\alpha}(x,s)=\sum_{s^{\prime}\in[K]}p_{s^{\prime}}Q_{q_{\alpha}\mid
s^{\prime}}\circ F_{q_{\alpha}\mid s}\circ q_{\alpha}(x,s).$ (6)

Proposition 1 implies that the optimal fair quantile predictor for an input
$(x,s)$ is obtained by a nonlinear transformation of the vector
$[q_{\alpha}(x,s)]_{s=1}^{K}$ linking to a Wasserstein barycenter problem[3,
13]. The explicit closed form solution comes from [3, 13, 18], which relies on
the classical characterization of optimal coupling in one dimension of the
Wasserstein-2 distance. A rigorous proof is given in [13, 26]. It shows that a
minimizer $g_{\alpha}$ of the $L_{2}$-risk can be used to construct $\nu$ and
vice-versa, given $\nu$, there is a explicit expression Eq. (6) for the multi-
marginal Wasserstein barycenter [3].

First, we provide a sketch of the CFQP approach. We start with splitting the
whole training data into a proper training set $\mathcal{I}_{1}$ and a
calibration set $\mathcal{I}_{2}$, then fit an arbitrary quantile regression
algorithm $\mathcal{Q}$ on $\mathcal{I}_{1}$,
$\left\\{\hat{q}_{\alpha_{\mathrm{lo}}},\hat{q}_{\alpha_{\mathrm{hi}}}\right\\}\leftarrow\mathcal{Q}\left(\left\\{\left(\tilde{X}_{i},Y_{i}\right):i\in\mathcal{I}_{1}\right\\}\right).$
We apply the fitted quantile algorithm $\mathcal{Q}$ on the calibration set
$\mathcal{I}_{2}$ to obtain the predicted
$\\{\hat{q}_{\alpha_{\mathrm{lo}}}(\tilde{X}_{i}),\hat{q}_{\alpha_{\mathrm{hi}}}(\tilde{X}_{i})\\}_{i\in\mathcal{I}_{2}}$.
Since the quantile estimates for $\mathcal{I}_{2}$ will be used for
conformalization, it is essential to transform them into fair ones, i.e.
$\hat{g}_{\alpha,i},\forall i\in\mathcal{I}_{2}$ (Eq. (9)), through Algorithm
2. Finally, for a test point $\tilde{X}_{n+1}$, we will predict two quantile
estimates $\hat{q}_{\alpha}(x,s)$ affected by the sensitive variable $S$ by
$\mathcal{Q}$, then apply the functional synchronization (details in Algorithm
2) and calibration (Algorithm 1) steps in turn to generate the fair constraint
prediction interval $C(\tilde{X}_{n+1})$ for $Y_{n+1}$.

Next, we explicate in detail how to remove the effect of the sensitive
variable for the predicted quantile estimates. By Proposition 1, the optimal
fair quantiles take the form of Eq. (6). Therefore, we propose an empirical
optimal fair quantile estimator $\hat{g}_{\alpha}$ that relies on the plug-in
principle. In particular, Eq.(6) indicates that for each quantile level
$\alpha$ and each category $s\in[K]$, we only need estimators for the
regression function $\hat{q}_{\alpha}$, the proportion $\hat{p}_{s}$, the
cumulative distribution function $F_{\hat{q}_{\alpha}\mid s}$ and the quantile
function $Q_{\hat{q}_{\alpha}\mid s}$.

Note that we can empirically estimate the CDF and quantile function for each
sensitive group in the calibration set $\mathcal{I}_{2}$ separately. Hence for
each quantile level $\alpha$, let $N_{s}:=|\mathcal{I}_{2}^{s}|$, and the
quantile estimators
$(\hat{q}_{1}^{s},\hat{q}_{2}^{s},\dots,\hat{q}_{N_{s}}^{s})$111$\hat{q}_{i}^{s}$
depends on the quantile level $\alpha$, we suppress $\alpha$ for notational
simplification. are calculated through the fitted quantile regression
$\mathcal{Q}$ with training set $\mathcal{I}_{1}$. We define the augmented
random variable for each data point in $\mathcal{I}_{2}$,

$\tilde{q}_{i}^{s}:=\hat{q}_{i}^{s}+U_{i}^{s}([-\sigma,\sigma])\quad\forall
i\in\mathcal{I}_{2}^{s},s\in[K],$

where $U_{i}^{s}$ are i.i.d. random variables, uniformly distributed on
$[-\sigma,\sigma]$ for some small positive $\sigma$, and independent from all
the previously introduced random variables. It serves as a smoothing random
variable, for the random variables $\tilde{q}_{i}^{s}$ are i.i.d. continuous
for any ${P}(Y|\tilde{X})$ and $\mathcal{Q}$. Otherwise, the original
$\hat{q}_{i}^{s}$ might have atoms resulting in a non-zero probability to
observe ties in the group $\\{\hat{q}_{i}^{s}\\}$ for $s=1,\dots,K$. This
trick, also called jittering [8, 13] is often used for data visualization for
tie-breaking. Using the above quantities, we build the CDF and quantile
function estimators for each subgroup $s^{\prime}\in[K]$ as follows,

$\displaystyle\hat{F}_{q_{\alpha}|s^{\prime}}(t)=N_{s}^{-1}\sum_{i=1}^{N_{s}}\mathbbm{1}\left\\{\tilde{q}_{i}^{s^{\prime}}\leq
t\right\\},\quad\text{for all }t\in\mathbb{R},$ (7)
$\displaystyle\hat{Q}_{2,q_{\alpha}|s^{\prime}}(t)=\int_{0}^{1}\hat{F}^{-1}_{q_{\alpha}|s^{\prime}}(v)K_{h}(t-v)dv,\quad
t\in(0,1).$ (8)

The smoothed kernel estimator Eq.(8) was firstly proposed by Cheng and Parzen
[9], where $K_{h}(\cdot)=K(\cdot/h)/h$ is a kernel function chosen as a
probability density function that is symmetric around zero with bandwidth
parameter $h>0$.

If the quantile functions ${Q}_{2,q_{\alpha}|s^{\prime}}$ is differentiable,
the derivative
$Q^{\prime}_{s^{\prime}}(t):={Q}_{2,q_{\alpha}|s^{\prime}}^{\prime}(t)$ for
$t\in(0,1)$ is the quantile density function [9, 52]. We hereby give an
estimation bound for Eq. (8) using kernel smoothing. For this purpose, we
invoke the conditions (A1) - (A3) that are needed for deducing the following
proposition. They can also be found from [52] and are included in the
supplementary material.

######  Proposition 2.

Under conditions (A1), (A2), and (A3), we have

$\sup_{s^{\prime}}\sup_{t\in[0,1]}\left|\hat{Q}_{2,q_{\alpha}|s^{\prime}}(t)-{Q}_{q_{\alpha}|s^{\prime}}(t)\right|=O_{p}\left(N^{-\nicefrac{{1}}{{2}}}\right),\quad
s^{\prime}=1,\ldots,K.$

The motivation for including a smoothing step is twofold: First, smoothing the
quantile function eliminates the troublesomeness in defining arbitrary
quantiles from the empirical one when the sample sizes of subgroups are small.
Second, the proposed kernel smoothing improves second-order efficiency by
alleviating the relative deficiency [17, 52].

######  Remark 1.

One can utilize various kernels such as the Gaussian or Epanechnikov kernel
with adaptive bandwidth for better practical performance. Other smoothing
methods such as splines or local linear fitting can likewise be applied with
equal effectiveness.

Consequently, for each quantile level $\alpha$, the functional synchronized
quantile estimator is

$\hat{g}_{\alpha,i}=\sum_{s^{\prime}=1}^{K}\hat{p}_{s^{\prime}}\hat{Q}_{2,q_{\alpha}|s^{\prime}}\circ\hat{F}_{q_{\alpha}|s}\circ\tilde{q}_{i}^{s},\quad\forall
i\in\mathcal{I}_{2}.$ (9)

The proposed estimator can be deemed as the empirical counterpart with
additional randomization of the explicit fair optimal formula Eq.(6).

To conformalize the adjusted fair quantiles Eq (9), we need to compute the
conformity scores $E_{i}$ for each $i\in\mathcal{I}_{2}$ that quantify the
error made by the plug-in fair prediction interval
$\hat{C}^{g}(\tilde{x})=[\hat{g}_{\alpha_{\mathrm{lo}}}(\tilde{x}),\hat{g}_{\alpha_{\mathrm{hi}}}(\tilde{x})]$.
The scores are evaluated on the calibration set as

$E_{i}:=\max\\{\hat{g}_{\alpha_{\mathrm{lo}},i}-Y_{i},Y_{i}-\hat{g}_{\alpha_{\mathrm{hi}},i}\\}.$
(10)

At the last stage, for a new data point $\tilde{X}_{n+1}=(x,s)$, and
$\alpha\in\\{\alpha_{\mathrm{lo}},\alpha_{\mathrm{hi}}\\}$, by defining

$\tilde{q}_{1,i}^{s}=\hat{q}_{i}^{s}+U_{i}^{s}([-\sigma,\sigma])\quad\forall
i\in\mathcal{I}_{1}^{s}\quad\text{ and
}\quad\tilde{q}_{\alpha}(x,s)=\hat{q}_{\alpha}(x,s)+U([-\sigma,\sigma]).$

We use the empirical CDF of training set 222Still, $\hat{q}_{i}^{s}$ depends
on quantile level $\alpha$.

$\textstyle{\hat{F}_{1,q_{\alpha}|s}(t):=\frac{1}{|\mathcal{I}_{1}^{s}|+1}\left(\sum_{i=1}^{|\mathcal{I}_{1}^{s}|}\mathbbm{1}\left\\{\tilde{q}_{1,i}^{s}<t\right\\}+U([0,1])\left(1+\sum_{i=1}^{|\mathcal{I}_{1}^{s}|}\mathbbm{1}\left\\{\tilde{q}_{1,i}^{s}=t\right\\}\right)\right)}$
(11)

to estimate the location
$\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s)$. Thus the fair quantile
estimator is built as follows

$\hat{g}_{\alpha}(x,s)=\sum_{s^{\prime}=1}^{K}\hat{p}_{s^{\prime}}\hat{Q}_{2,q_{\alpha}|s^{\prime}}\circ\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s),\forall\alpha\in\\{\alpha_{\mathrm{lo}},\alpha_{\mathrm{hi}}\\}.$
(12)

The fair prediction interval for $Y_{n+1}$ is constructed as

$C(\tilde{X}_{n+1})=[\hat{g}_{\alpha_{\mathrm{lo}}}(x,s)-Q_{1-\alpha}(E,\mathcal{I}_{2}),\hat{g}_{\alpha_{\mathrm{hi}}}(x,s)+Q_{1-\alpha}(E,\mathcal{I}_{2})],$
(13)

where $Q_{1-\alpha}(E,\mathcal{I}_{2}):=(1-\alpha)(1+1/|\mathcal{I}_{2}|)$-th
empirical quantile of $\\{E_{i}:i\in\mathcal{I}_{2}\\}$ will adjust the plug-
in fair prediction interval. We present the pseudo-codes of CFQP as well as
the construction of $\hat{g}_{\alpha}$ for Eq. 9 in Algorithm 1, 2
respectively.

Algorithm 1 Split Conformal Fair Prediction (CFQP)

0: $\mathcal{D}=\\{(X_{i},S_{i},Y_{i})\\}_{i=1}^{n}$; miscoverage level
$\alpha\in(0,1)$; quantile regression algorithm $\mathcal{Q}$.

1: Randomly split $[n]$ into disjoint proper training and calibration indices
$\mathcal{I}_{1}$, $\mathcal{I}_{2}$.

2: Fit two conditional quantile functions on the training set
$\\{\hat{q}_{\alpha_{\mathrm{lo}}},\hat{q}_{\alpha_{\mathrm{hi}}}\\}\leftarrow\mathcal{Q}(\\{(X_{i},S_{i},Y_{i}),i\in\mathcal{I}_{1}\\})$.

3: Call functional Synchronization (Algorithm 2) to calculate
$\\{\hat{g}_{\alpha_{\mathrm{lo}}},\hat{g}_{\alpha_{\mathrm{hi}}}\\}$ for each
$i\in\mathcal{I}_{2}$.

4: Compute
$E_{i}\leftarrow\max\left\\{\hat{g}_{\alpha_{\mathrm{lo}}}\left(X_{i}\right)-Y_{i},Y_{i}-\hat{g}_{\alpha_{\mathrm{hi}}}\left(X_{i}\right)\right\\}$
for $\forall\ i\in\mathcal{I}_{2}$.

5: Compute
$Q_{1-\alpha}(E,\mathcal{I}_{2})\leftarrow(1-\alpha)(1+\nicefrac{{1}}{{|\mathcal{I}_{2}|}})$-th
empirical quantile of $\\{E_{i}:i\in\mathcal{I}_{2}\\}$.

6: For a new test point $(x,s)$, compute
$\\{\hat{g}_{\alpha_{\mathrm{lo}}}(x,s),\hat{g}_{\alpha_{\mathrm{hi}}}(x,s)\\}$
through Algorithm 2

6: Fair prediction interval
$C(x,s)=[\hat{g}_{\alpha_{\mathrm{lo}}}(x,s)-Q_{1-\alpha}(E,\mathcal{I}_{2}),\hat{g}_{\alpha_{\mathrm{hi}}}(x,s)+Q_{1-\alpha}(E,\mathcal{I}_{2})]$
for $(X_{n+1},S_{n+1})=(x,s)$.

Algorithm 2 Functional Synchronization

0: Calibration set $\\{(X_{i},S_{i})\\}_{i\in\mathcal{I}_{2}}$ or new point
$(x,s)$; base quantile estimator $\mathcal{Q}$;slack parameter $\sigma$;
training set $\\{(X_{i},S_{i})\\}_{i\in\mathcal{I}_{1}}$;

1: if Calibration set $\\{(X_{i},S_{i})\\}_{i\in\mathcal{I}_{2}}$ then

2: for $\alpha\in\\{\alpha_{\mathrm{lo}},\alpha_{\mathrm{hi}}\\}$ do

3:
$\\{\tilde{q}_{\alpha}(X_{i},S_{i})\\}\leftarrow\\{q_{\alpha}(X_{i},S_{i})+U([-\sigma,\sigma])\\}_{i\in\mathcal{I}_{2}}$
$\triangleright\ U([-\sigma,\sigma])$ are used for tie-breaking

4: for $s^{\prime}\in[K]$ do

5: Compute $\hat{F}_{q_{\alpha}|s^{\prime}}(t)$, and
$\hat{F}^{-1}_{2,q_{\alpha}|s^{\prime}}(t)$ by Eq. (7) and (8).

6: Obtain
$\hat{g}_{\alpha}(X_{i},S_{i})\leftarrow\sum_{s^{\prime}=1}^{K}\hat{p}_{s^{\prime}}\hat{F}^{-1}_{2,q_{\alpha}|s^{\prime}}\circ\hat{F}_{q_{\alpha}|s^{\prime}}\circ\tilde{q}_{\alpha}(X_{i},S_{i}),\
\forall i\in\mathcal{I}_{2}$

7: end for

8: end for

9: else if New test point $(x,s)$ then

10: for $\alpha\in\\{\alpha_{\mathrm{lo}},\alpha_{\mathrm{hi}}\\}$ do

11:
$\\{\tilde{q}^{s}_{1,\alpha}\\}\leftarrow\\{\hat{q}_{\alpha}^{s}+U([-\sigma,\sigma])\\}_{i\in\mathcal{I}_{1}^{s}}\
\text{and}\ \tilde{q}_{\alpha}(x,s)\leftarrow
q_{\alpha}(x,s)+U([-\sigma,\sigma])$

12: Compute
$\hat{g}_{\alpha}(x,s)\leftarrow\sum_{s^{\prime}=1}^{K}\hat{p}_{s^{\prime}}\hat{F}^{-1}_{2,q_{\alpha}|s^{\prime}}\circ\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}({x},s)$
by Eq. (8) and (7)

13: end for

14: end if

14: fair quantile prediction $\hat{g}_{\alpha}$ for calibration set or new
test point $(x,s)$.

##  5 Theoretical results

We provide a statistical analysis of the proposed algorithm with coverage and
DP-fairness guarantees in this part.

######  Theorem 1 (Prediction coverage guarantee).

If $(\tilde{X}_{i},Y_{i}),i=1,\dots,n+1$ are exchangeable, then the prediction
interval $C(\tilde{X}_{n+1})$ constructed by the split CFQP algorithm
satisfies

${P}\\{Y_{n+1}\in C(\tilde{X}_{n+1})\\}\geq 1-\alpha.$

Moreover, if the conformity scores $E_{i}$ are almost surely distinct, the
prediction interval is nearly exactly calibrated,

${P}\\{Y_{n+1}\in C(\tilde{X}_{n+1})\\}\leq 1-\alpha+1/(|\mathcal{I}_{2}|+1).$

######  Remark 2.

Corollary 1 in the supplementary material gives an extension for the
conformalization step which allows coverage errors to be spread arbitrarily
over the left and right tails. Controlling the left and right tails
independently yields a stronger coverage guarantee.

######  Theorem 2 (Demographic parity guarantee).

For any joint distribution ${P}$ of $({X},S,Y)$, any $\sigma>0$, as well as
the base quantile estimator
$\hat{q}_{\alpha}:\mathbb{R}^{p}\times[K]\rightarrow\mathbb{R}$ constructed on
labeled data, the estimator $\hat{g}_{\alpha}$ defined in Eq. (12) satisfies

$(\hat{g}_{\alpha}({X},S)\mid
S=s)\overset{d}{=}\left(\hat{g}_{\alpha}({X},S)\mid
S=s^{\prime}\right)\quad\forall s,s^{\prime}\in[K].$

Quantile DP guarantees provided by Theorem 2 are derived directly from
distribution-free properties on rank and order statistics in Lemma E.1.,
Theorem 7.2, of Chzhen and Schreuder [12]. Further information can be found in
their papers and the references they provide. We extend the estimator
$\hat{g}$ for quantile regression based on the estimators developed in their
work with solid theoretical foundations. In contrast to the seminal work of
Chzhen and Schreuder [12], we regard the densities and quantile functions as
functional data, namely, as samples of stochastic processes. Typically, this
approach is used by economists when dealing with the densities of income
distribution across populations. As a result of introducing the kernel
smoothing step, potential performance improvements are demonstrated in the
supplemental material.

##  6 Experiments

To evaluate our proposed method 333We utilize the local linear fitting
smoothing method in the experiments., we report the performance of post-
processing fairness adjustment on quantiles through four benchmark datasets:
Law School (LAW), Community&Crime (CRIME), MEPS 2016 (MEPS), Government Salary
(GOV). A detailed description of these datasets can be found in the
supplementary material. The code for reproducing our results is avaiable at
https://github.com/Lei-Ding07/Conformal_Quantile_Fairness.

| LAW | CRIME
---|---|---
Coverage |  Length |  KS(lo) |  KS(hi) |  Coverage |  Length |  KS(lo) |  KS(hi)
Ln-CQR |  90.16$\pm$0.47 |  0.46$\pm$.004 |  0.39$\pm$0.03 |  0.11$\pm$0.02 |  90.22$\pm$1.88 |  1.30$\pm$0.05 |  0.62$\pm$0.06 |  0.53$\pm$0.06
Ln-CFQP |  90.02$\pm$0.51 |  0.46$\pm$.004 |  0.02$\pm$0.01 |  0.02$\pm$0.01 |  90.44$\pm$1.84 |  1.64$\pm$0.05 |  0.11$\pm$0.03 |  0.12$\pm$0.04
RF-CQR |  90.25$\pm$0.55 |  0.39$\pm$.005 |  0.20$\pm$0.02 |  0.15$\pm$0.02 |  90.27$\pm$1.66 |  1.15$\pm$0.03 |  0.64$\pm$0.05 |  0.59$\pm$0.05
RF-CFQP |  90.11$\pm$0.48 |  0.38$\pm$.004 |  0.02$\pm$.008 |  0.02$\pm$.009 |  90.34$\pm$1.84 |  1.54$\pm$0.04 |  0.12$\pm$0.04 |  0.12$\pm$0.03
NN-CQR |  90.00$\pm$0.50 |  0.40$\pm$0.02 |  0.41$\pm$0.07 |  0.18$\pm$0.05 |  90.01$\pm$1.89 |  1.16$\pm$0.05 |  0.70$\pm$0.05 |  0.63$\pm$0.06
NN-CFQP |  90.01$\pm$0.51 |  0.39$\pm$0.01 |  0.02$\pm$.009 |  0.03$\pm$.009 |  89.95$\pm$1.62 |  1.54$\pm$0.12 |  0.12$\pm$0.04 |  0.12$\pm$0.03
| MEPS | GOV
Coverage |  Length |  KS (lo) |  KS(hi) |  Coverage |  Length |  KS (lo) |  KS(hi)
Ln-CQR |  89.92$\pm$0.66 |  0.66$\pm$0.01 |  0.09$\pm$0.03 |  0.33$\pm$0.05 |  90.00$\pm$0.19 |  0.79$\pm$.002 |  0.26$\pm$.014 |  0.44$\pm$0.02
Ln-CFQP |  89.99$\pm$0.69 |  0.66$\pm$0.01 |  0.03$\pm$0.01 |  0.03$\pm$0.01 |  90.02$\pm$0.19 |  0.78$\pm$.002 |  0.05$\pm$0.01 |  0.04$\pm$0.01
RF-CQR |  90.07$\pm$0.65 |  0.38$\pm$.009 |  0.19$\pm$0.02 |  0.30$\pm$0.03 |  90.03$\pm$0.17 |  0.61$\pm$.002 |  0.29$\pm$0.01 |  0.28$\pm$0.02
RF-CFQP |  90.38$\pm$0.60 |  0.39$\pm$0.01 |  0.02$\pm$0.01 |  0.03$\pm$0.01 |  90.03$\pm$0.17 |  0.62$\pm$.002 |  0.05$\pm$0.01 |  0.04$\pm$0.01
NN-CQR |  89.95$\pm$0.68 |  0.37$\pm$0.04 |  0.24$\pm$0.09 |  0.37$\pm$ 0.06 |  90.01$\pm$0.19 |  0.58$\pm$0.01 |  0.28$\pm$0.03 |  0.32$\pm$0.04
NN-CFQP |  89.97$\pm$0.61 |  0.37$\pm$0.04 |  0.03$\pm$0.01 |  0.04$\pm$0.01 |  90.01$\pm$0.18 |  0.59$\pm$0.01 |  0.05$\pm$0.01 |  0.05$\pm$0.01

Table 1: Results reported on test set of 200 repeated experiments with
$\alpha=0.1$. CQR refers to the conformalized quantile regression in [38]. Ln,
RF, and NN denote the linear, random forest, and neural network quantile
regression models. Our methods are shown in bold.

Figure 2:  Results for estimating the lower ($\alpha_{\mathrm{lo}}$) and upper
($\alpha_{\mathrm{hi}}$) quantiles using some state-of-the-art DP-fairness
requirement methods on all the datasets. ‘Unfair’, ‘Chzhen’, and ‘Agarwal’
stand for the linear quantile model without fairness adjustment, barycenter
method [13] and reduction-based algorithm [2] respectively. We present the MAE
and KS of lower and upper quantile estimations. A Linear quantile model is
implemented in this comparison.

We measure the violation of DP-fairness of the quantiles required by
Definition 1 through the empirical Kolmogorov-Smirnov (KS) distance. The value
represents the disparity between groups
$\mathcal{Z}^{s}=\\{(X,S,Y)\in\mathcal{Z}:S=s\\}$ for all $s\in[K]$,

$\textstyle{\mathrm{KS}(g_{\alpha})=\underset{s,s^{\prime}\in[K]}{\max}\
\underset{t\in\mathbb{R}}{\sup}\left\lvert\frac{1}{|\mathcal{Z}^{s}|}\underset{(X,S,Y)\in\mathcal{Z}^{s}}{\sum}\mathbbm{1}{\left\\{g_{\alpha}({X,S})\leq
t\right\\}}-\frac{1}{|\mathcal{Z}^{s^{\prime}}|}\underset{(X,S,Y)\in\mathcal{Z}^{s^{\prime}}}{\sum}\mathbbm{1}{\left\\{g_{\alpha}({X,S})\leq
t\right\\}}\right\rvert}.$

Experiment results. In Table 1, we report the average performance of the
proposed CFQP over 200 randomly training-test splits as well as the baseline
model CQR by the coverage rate, length of prediction interval, and the KS
distance of the interval endpoint. We split the training data into proper
training and calibration sets of equal sizes. Throughout the experiments, the
nominal miscoverage rate is fixed to $\alpha=0.1$. Among pre-existing quantile
algorithms, we select three leading variants: linear model[27], random forests
[30] and neural networks [39]. Overall, our CFQP likewise CQR constructs
prediction bands attaining desirable coverage around 90%, as claimed in
Theorem 1. Random forest based approaches tend to be slightly more
conservative than the other two w.r.t the coverage rate among all four
datasets.

In the KS column concerning the DP fairness of interval endpoints, our CFQP
method greatly reduces the discriminatory bias (quantified by KS) by 70% up to
90% compared to that of CQR. In addition, the lengths of the prediction
intervals mostly remain the same except for the Crime dataset, probably due to
its inherent high discriminatory bias between sensitive groups. In addition,
CFQP is more robust according to the standard errors over experiment
repetitions.

Figure 2 presents the comparison of our post-processing fairness adjustment
procedure on quantiles $\hat{g}_{\alpha}$ using the test set
$\mathcal{Z}_{test}=\\{(X_{i},S_{i},Y_{i})\\}_{i=1}^{n_{test}}$ with some
state-of-the-art fairness algorithms. Since most of the algorithms are
targeted at mean prediction, there is no direct comparison with our quantile
fairness method; we accordingly modified the existing methods into quantile
versions for comparison. A detailed description can be found in the
supplementary material.

The points in Figure 2 represents the mean of 200 repeated experiments with
$x$-axis as KS distance and $y$-axis as Mean Absolute Error(MAE), where
$\operatorname{MAE}(g_{\alpha})=\nicefrac{{1}}{{n_{test}}}\sum_{\mathcal{Z}_{test}}|Y-g_{\alpha}(X,S)|$
measures the prediction error of quantiles, the bars are the standard error on
both axes. The optimal points should locate at the bottom left corner of the
graph, where smaller KS distance and smaller MAE are achieved. In each
subplot, our method consistently performs better with the smallest KS distance
while keeping the MAE equal or slightly higher than the others or the unfair
version. Note that due to the highly right skewness of real datasets, the MAE
of the upper quantile estimation is larger than that of the lower quantile for
all quantile approaches. Additional ablation studies for computational time
analysis and testing the kernel smoothing approach are incorporated in the
supplementary material.

##  7 Conclusion and future work

Conformal fair quantile regression is a novel approach for creating fair
prediction intervals that attain valid coverage and reach independence between
sensitive attributes while making minimal modifications to the quantile
endpoints simultaneously. It becomes superior within heteroskedastic and/or
asymmetric datasets and robust to outliers.

Our method is supported by rigorous distribution-free coverage and exact DP-
fairness guarantees, as proved in theoretical parts. We conducted several real
data examples demonstrating the effectiveness of our method in achieving exact
coverage while imposing DP-fairness in practice. The method outperforms
several state-of-the-art approaches by comparison.

A limitation in our numerical experiments is that we simply utilize the local
linear smoothing method in defining quantile functions of the subgroups; we
believe incorporating flexible kernel smoothing approaches [47, 52] would
improve the experimental performances.

As potential future works, it would be valuable to introduce a DP relaxation
framework based on an unfairness measure in a similar manner as [12, 44],
allowing controlling the level of unfairness in quantile estimates. We also
expect to extend the scope to other potential fairness metrics which is
dependent on the underlying response like equalizing quantile loss across
groups by incorporating a fairness penalty term in training, or the fairness
metric defined for conditional variance-at-risk. In addition, it is worthwhile
exploring the quantile fairness in other areas of AI, such as NLP[16, 25],
Computer Vision, Recommendation systems, etc.

## Acknowledgements

This work was supported by the Economic and Social Research Council (ESRC
ES/T012382/1) and the Social Sciences and Humanities Research Council (SSHRC
2003-2019-0003) under the scheme of the Canada-UK Artificial Intelligence
Initiative. The project title is BIAS: Responsible AI for Gender and Ethnic
Labour Market Equality. We thank Yang Hu, Xiaojun Du and Matthew Pietrosanu
for their helpful discussion and valuable input and all the constructive
suggestions and comments from the reviewers.

## References

  * [1] Medical expenditure panel survey.  URL https://meps.ahrq.gov/mepsweb/.
  * Agarwal et al. [2019] A. Agarwal, M. Dudík, and Z. S. Wu.  Fair regression: Quantitative definitions and reduction-based algorithms.  In _International Conference on Machine Learning_ , page 120–129. PMLR, 2019.
  * Agueh and Carlier [2011] M. Agueh and G. Carlier.  Barycenters in the wasserstein space.  _SIAM Journal on Mathematical Analysis_ , 43(2):904–924, 2011.  doi: 10.1137/100805741.
  * Barocas et al. [2017] S. Barocas, M. Hardt, and A. Narayanan.  Fairness in machine learning.  _Nips tutorial_ , 1:2, 2017.
  * Berk et al. [2017] R. Berk, H. Heidari, S. Jabbari, M. Joseph, M. Kearns, J. Morgenstern, S. Neel, and A. Roth.  A convex framework for fair regression.  _arXiv preprint arXiv:1706.02409_ , 2017.
  * Calmon et al. [2017] F. Calmon, D. Wei, B. Vinzamuri, K. Natesan Ramamurthy, and K. R. Varshney.  Optimized pre-processing for discrimination prevention.  _Advances in neural information processing systems_ , 30, 2017.
  * Castelnovo et al. [2022] A. Castelnovo, R. Crupi, G. Greco, D. Regoli, I. G. Penco, and A. C. Cosentini.  A clarification of the nuances in the fairness metrics landscape.  _Scientific Reports_ , 12(1):1–21, 2022.
  * Chambers et al. [2018] J. M. Chambers, W. S. Cleveland, B. Kleiner, and P. A. Tukey.  _Graphical methods for data analysis_.  Chapman and Hall/CRC, 2018.
  * Cheng and Parzen [1997] C. Cheng and E. Parzen.  Unified estimators of smooth quantile and quantile density functions.  _Journal of statistical planning and inference_ , 59(2):291–307, 1997.
  * Cho et al. [2020] J. Cho, G. Hwang, and C. Suh.  A fair classifier using kernel density estimation.  _Advances in Neural Information Processing Systems_ , 33:15088–15099, 2020.
  * Chouldechova [2017] A. Chouldechova.  Fair prediction with disparate impact: A study of bias in recidivism prediction instruments.  _Big data_ , 5(2):153–163, 2017.
  * Chzhen and Schreuder [2022] E. Chzhen and N. Schreuder.  A minimax framework for quantifying risk-fairness trade-off in regression.  _The Annals of Statistics_ , 50(4):2416–2442, 2022.
  * Chzhen et al. [2020] E. Chzhen, C. Denis, M. Hebiri, L. Oneto, and M. Pontil.  Fair regression with wasserstein barycenters.  _Advances in Neural Information Processing Systems_ , 33:7321–7331, 2020.
  * Csorgo and Revesz [1978] M. Csorgo and P. Revesz.  Strong approximations of the quantile process.  _The Annals of Statistics_ , 6(4):882–894, 1978\.  ISSN 00905364.
  * Czarnowska et al. [2021] P. Czarnowska, Y. Vyas, and K. Shah.  Quantifying social biases in nlp: A generalization and empirical comparison of extrinsic fairness metrics.  _Transactions of the Association for Computational Linguistics_ , 9:1249–1267, 2021.
  * Ding et al. [2022] L. Ding, D. Yu, J. Xie, W. Guo, S. Hu, M. Liu, L. Kong, H. Dai, Y. Bao, and B. Jiang.  Word embeddings via causal inference: Gender bias reducing and semantic information preserving.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, pages 11864–11872, 2022.
  * Falk [1984] M. Falk.  Relative deficiency of kernel type estimators of quantiles.  _The Annals of Statistics_ , pages 261–268, 1984.
  * Gallón et al. [2013] S. Gallón, J.-M. Loubes, and E. Maza.  Statistical properties of the quantile normalization method for density curve alignment.  _Mathematical biosciences_ , 242(2):129–142, 2013.  doi: 10.1016/j.mbs.2012.12.007.
  * García et al. [2001] J. García, P. J. Hernández, and A. López-Nicolás.  How wide is the gap? an investigation of gender wage differences using quantile regression.  _Empirical Economics_ , 26(1):149–167, Mar 2001\.  doi: 10.1007/s001810000050.
  * Gardeazabal and Ugidos [2005] J. Gardeazabal and A. Ugidos.  Gender wage discrimination at quantiles.  _Journal of Population Economics_ , 18(1):165–179, Mar 2005.  doi: 10.1007/s00148-003-0172-z.
  * Gouic et al. [2020] T. L. Gouic, J.-M. Loubes, and P. Rigollet.  Projection to fairness in statistical learning.  _arXiv:2005.11720 [cs, math, stat]_ , Jun 2020.  arXiv: 2005.11720.
  * Han et al. [2019] P. Han, L. Kong, J. Zhao, and X. Zhou.  A general framework for quantile estimation with incomplete data.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 81(2):305–333, 2019.
  * Harris et al. [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, Sept. 2020.  doi: 10.1038/s41586-020-2649-2.
  * Holzer and Neumark [2006] H. J. Holzer and D. Neumark.  Affirmative action: What do we know?  _Journal of policy analysis and management_ , 25(2):463–490, 2006.
  * Hu et al. [2022] S. Hu, J. A. Al-Ani, K. D. Hughes, N. Denier, A. Konnikov, L. Ding, J. Xie, Y. Hu, M. Tarafdar, B. Jiang, et al.  Balancing gender bias in job advertisements with text-level bias mitigation.  _Frontiers in big Data_ , 5, 2022.
  * Jiang et al. [2020] R. Jiang, A. Pacchiano, T. Stepleton, H. Jiang, and S. Chiappa.  Wasserstein fair classification.  In _Proceedings of The 35th Uncertainty in Artificial Intelligence Conference_ , page 862–872. PMLR, Aug 2020.
  * Koenker and Bassett [1978] R. Koenker and G. Bassett.  Regression quantiles.  _Econometrica: journal of the Econometric Society_ , pages 33–50, 1978.
  * Kusner et al. [2017] M. J. Kusner, J. Loftus, C. Russell, and R. Silva.  Counterfactual fairness.  _Advances in neural information processing systems_ , 30, 2017.
  * Lehmann et al. [2005] E. L. Lehmann, J. P. Romano, and G. Casella.  _Testing statistical hypotheses_ , volume 3.  Springer, 2005.
  * Meinshausen and Ridgeway [2006] N. Meinshausen and G. Ridgeway.  Quantile regression forests.  _Journal of Machine Learning Research_ , 7(6), 2006.
  * Milli et al. [2019] S. Milli, J. Miller, A. D. Dragan, and M. Hardt.  The social cost of strategic classification.  In _Proceedings of the Conference on Fairness, Accountability, and Transparency_ , pages 230–239, 2019.
  * Oneto et al. [2019] L. Oneto, M. Donini, and M. Pontil.  General fair empirical risk minimization.  _arXiv:1901.10080 [cs, stat]_ , Dec 2019.  arXiv: 1901.10080.
  * Pietrosanu et al. [2021] M. Pietrosanu, J. Gao, L. Kong, B. Jiang, and D. Niu.  Advanced algorithms for penalized quantile and composite quantile regression.  _Computational Statistics_ , 36(1):333–346, Mar 2021.  ISSN 0943-4062, 1613-9658.  doi: 10.1007/s00180-020-01010-1.
  * Plečko and Meinshausen [2020] D. Plečko and N. Meinshausen.  Fair data adaptation with quantile preservation.  _Journal of Machine Learning Research_ , 21(242):1–44, 2020.
  * Plečko et al. [2021] D. Plečko, N. Bennett, and N. Meinshausen.  fairadapt: Causal reasoning for fair data pre-processing.  _arXiv preprint arXiv:2110.10200_ , 2021.
  * Redmond and Baveja [2002] M. Redmond and A. Baveja.  A data-driven software tool for enabling cooperative information sharing among police departments.  _European Journal of Operational Research_ , 141(3):660–678, 2002.
  * Romano et al. [2019a] Y. Romano, R. Barber, C. Sabatti, and E. Candès.  With malice towards none: Assessing uncertainty via equalized coverage.  _ArXiv_ , 2019a.  doi: 10.1162/99608f92.03f00592.
  * Romano et al. [2019b] Y. Romano, E. Patterson, and E. Candes.  Conformalized quantile regression.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019b.
  * Taylor [2000] J. W. Taylor.  A quantile regression neural network approach to estimating the conditional density of multiperiod returns.  _Journal of Forecasting_ , 19(4):299–311, 2000\.
  * Vovk et al. [2005] V. Vovk, A. Gammerman, and G. Shafer.  _Algorithmic learning in a random world_.  Springer, New York, 2005.  ISBN 978-0-387-00152-4.
  * Vovk et al. [2009] V. Vovk, I. Nouretdinov, and A. Gammerman.  On-line predictive linear regression.  _The Annals of Statistics_ , pages 1566–1590, 2009.
  * Wang et al. [2019] Y. Wang, L. Kong, B. Jiang, X. Zhou, S. Yu, L. Zhang, and G. Heo.  Wavelet-based lasso in functional linear quantile regression.  _Journal of Statistical Computation and Simulation_ , 89(6):1111–1130, 2019.
  * Wightman [1998] L. F. Wightman.  Lsac national longitudinal bar passage study. lsac research report series.  1998\.
  * Williamson and Menon [2019] R. Williamson and A. Menon.  Fairness risk measures.  In _International Conference on Machine Learning_ , pages 6786–6797. PMLR, 2019.
  * Xiao and Koenker [2009] Z. Xiao and R. Koenker.  Conditional quantile estimation for generalized autoregressive conditional heteroscedasticity models.  _Journal of the American Statistical Association_ , 104(488):1696–1712, 2009.
  * Yang et al. [2019] D. Yang, J. Lafferty, and D. Pollard.  Fair quantile regression.  _arXiv preprint arXiv:1907.08646_ , 2019.
  * Yang [1985] S.-S. Yang.  A smooth nonparametric estimator of a quantile function.  _Journal of the American Statistical Association_ , 80(392):1004–1011, 1985.  ISSN 0162-1459.  doi: 10.2307/2288567.
  * Yu et al. [2016] D. Yu, L. Kong, and I. Mizera.  Partial functional linear quantile regression for neuroimaging data analysis.  _Neurocomputing_ , 195:74–87, 2016.
  * Yu et al. [2019] D. Yu, L. Zhang, I. Mizera, B. Jiang, and L. Kong.  Sparse wavelet estimation in quantile regression with multiple functional predictors.  _Computational statistics & data analysis_, 136:12–29, 2019.
  * Zemel et al. [2013] R. Zemel, Y. Wu, K. Swersky, T. Pitassi, and C. Dwork.  Learning fair representations.  In _International conference on machine learning_ , pages 325–333. PMLR, 2013.
  * Zeni et al. [2020] G. Zeni, M. Fontana, and S. Vantini.  Conformal prediction: a unified review of theory and new challenges.  _arXiv preprint arXiv:2005.07972_ , 2020.
  * Zhang and Müller [2011] Z. Zhang and H.-G. Müller.  Functional density synchronization.  _Computational Statistics & Data Analysis_, 55(7):2234–2249, 2011.  doi: 10.1016/j.csda.2011.01.007.

Supplementary Material

##  Appendix A Conditions for Proposition 2.

_Condition (A1)._ The quantile density functions
$Q^{\prime}_{s^{\prime}},s^{\prime}=1,\dots,K$, are twice continuously
differentiable in $(0,1)$, and satisfy
$\inf_{t\in[0,1]}Q^{\prime}_{s^{\prime}}(t)\geq c_{0}>0$, for a constant
$c_{0}>0$, $\forall s^{\prime}$. There is a $\gamma>0$ such that
$\sup_{u\in[0,1]}u(1-u)|J_{s^{\prime}}(u)|\leq\gamma$ for all $s^{\prime}$,
where $J_{s^{\prime}}(u):=[\mathrm{d}\log
Q^{\prime}_{s^{\prime}}(u)/\mathrm{d}u]$.

Condition (A2). There exists $0<L_{0}<\infty$, such that
$\sup_{u\in[0,1]}|\int_{0}^{1}Q^{\prime}_{s^{\prime}}(t)K_{h}(u-t)\mathrm{d}t|\leq
L_{0},\forall k$.

Condition (A3). The kernel functions $K$ is probability density functions
which is symmetric around $0$. For any function $f$ that is at least twice
continuously differentiable in $(0,1)$, it holds that for a $\rho>1/2$,
$\lim\sup_{N\rightarrow\infty}h^{2}N^{\rho}<\infty$, and
$\sup_{u\in[a,b]\subset(0,1)}\left|f(u)-\int_{0}^{1}f(t)K_{h}(u-t)\mathrm{d}t\right|=O\left(N^{-\rho}\right),$
where $N$ is defined to ensure that the numbers of measurement asymptotically
increases in the same way across the groups, we assume that there exists a
sequence $N=N(n)$ with $N\rightarrow\infty,$ as $n\rightarrow\infty$, such
that $\nicefrac{{N_{s^{\prime}}}}{{N}}\rightarrow\tau_{s^{\prime}}$ for
positive constants, and $0<c_{0}\leq\inf_{1\leq s^{\prime}\leq
K}\tau_{s^{\prime}}\leq\sup_{1\leq s^{\prime}\leq
K}\tau_{s^{\prime}}<C_{0}<\infty,s^{\prime}=1,\dots,K$.

######  Remark 3.

Conditions (A1) - (A3) guarantee the existence of a strong approximation of
the empirical quantile process by a sequence of weighted Brownian bridges as
established in [14]. Condition (A3) posited on kernel functions assures that
the integral transform $\hat{F}^{-1}\mapsto\hat{Q}$ possesses good
approximation properties for smooth functions, and it is shown that (A3) holds
for any difference kernel
$\mathrm{d}_{t}K_{h}(u,t)=h^{-1}k\left((u-t)/h\right)\mathrm{d}t$ with a
vanishing bandwidth $h$. for example, the gaussian density
$K(u)=\exp(-{u^{2}}/{2h^{2}})$ [9, 47, 52] or the triangular density function
$K(u)=(1-|u|/h)I(|u|/h\leq 1)$ with a vanishing bandwidth $h_{n}$.

##  Appendix B Proofs related to CFQP

In this section, we prove the validity of the CFQP prediction intervals
described in Section 4 of the main paper. First, we recap some results on
distribution-free order statistics from Romano et al. [37].

The quantile function $Q$ of a random variable $Z$, with cumulative
distribution function $F(z):={P}\\{Z\leq z\\}$, is defined by the equivalence

$Q(\alpha)\leq z\quad\text{ if and only if }\quad\alpha\leq F(z)$

for all $\alpha\in(0,1)$ and $z\in\mathbb{R}$. And less standardly, the right
quantile function $R$ of the random variable $Z$ is defined by the equivalence
$F^{-}(z)\leq\alpha\text{ if and only if }z\leq R(\alpha),$ where
$F^{-}(z):=F(z-)={P}\\{Z<z\\}$. The quantile functions have the explicit
formulas

$Q(\alpha)=\inf\\{z\in\mathbb{R}:\alpha\leq F(z)\\},\quad
R(\alpha)=\sup\left\\{z\in\mathbb{R}:F^{-}(z)\leq\alpha\right\\}.$

As a special case, the empirical quantile function $\hat{Q}_{n}$ of random
variables $Z_{1},\ldots,Z_{n}$ is the quantile function with respect to the
empirical
$\operatorname{CDF}\hat{F}_{n}(z):=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\\{Z_{i}\leq
z\\}.$ Likewise, the right empirical quantile function $\hat{R}_{n}$ of
$Z_{1},\ldots,Z_{n}$ is the right quantile function with respect to
$\hat{F}_{n}^{-}(z)=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\\{Z_{i}<z\\}$. They
have the explicit formulas

$\hat{Q}_{n}(\alpha)=Z_{(\lceil\alpha
n\rceil)},\quad\hat{R}_{n}(\alpha)=Z_{(\lfloor\alpha n\rfloor+1)},$

where $Z_{(k)}$ denotes the $k$ th smallest value in $Z_{1},\ldots,Z_{n}$.
Several variants of the following lemmas are showed in Vovk et al. [40],Vovk
et al. [41] and Romano et al. [38]. We restate them here for clarification.

######  Lemma 1 (Quantiles and exchangeability).

Suppose $Z_{1},\ldots,Z_{n}$ are exchangeable random variables. For any
$\alpha\in(0,1)$,

${P}\left\\{Z_{n}\leq\hat{Q}_{n}(\alpha)\right\\}\geq\alpha.$

Moreover, if the random variables $Z_{1},\ldots,Z_{n}$ are almost surely
distinct, then also

${P}\left\\{Z_{n}\leq\hat{Q}_{n}(\alpha)\right\\}\leq\alpha+\frac{1}{n}.$

Here the probabilities are taken over all the variables $Z_{1},\ldots,Z_{n}$.

######  Lemma 2 (Inflation of quantiles).

Suppose $Z_{1},\ldots,Z_{n+1}$ are exchangeable random variables. For any
$\alpha\in(0,1)$,

${P}\left\\{Z_{n+1}\leq\hat{Q}_{n}\left(\left(1+\frac{1}{n}\right)\alpha\right)\right\\}\geq\alpha.$

Moreover, if the random variables $Z_{1},\ldots,Z_{n+1}$ are almost surely
distinct, then also

${P}\left\\{Z_{n+1}\leq\hat{Q}_{n}\left(\left(1+\frac{1}{n}\right)\alpha\right)\right\\}\leq\alpha+\frac{1}{n+1}.$

Applying the previous auxiliary lemmas, we can prove the validity of the CFQP
prediction intervals described in Section 4.

Theorem 1. If $(\tilde{X}_{i},Y_{i}),i=1,\dots,n+1$ are exchangeable, then the
prediction interval $C(\tilde{X}_{n+1})$ constructed by the split CFQP
algorithm satisfies

${P}\\{Y_{n+1}\in C(\tilde{X}_{n+1})\\}\geq 1-\alpha.$

Moreover, if the conformity scores $E_{i}$ are almost surely distinct, the
prediction interval is nearly exactly calibrated,

${P}\\{Y_{n+1}\in C(\tilde{X}_{n+1})\\}\leq 1-\alpha+1/(|\mathcal{I}_{2}|+1).$

###### Proof of Theorem 1..

Conditionally on the proper training set. Denote by $E_{n+1}$ the conformity
score

$E_{i}:=\max\left\\{\hat{g}_{\alpha_{\mathrm{lo}}}\left(\tilde{X}_{i}\right)-Y_{i},Y_{i}-\hat{g}_{\alpha_{\mathrm{hi}}}\left(\tilde{X}_{i}\right)\right\\}$

at the test point $\left(\tilde{X}_{n+1},Y_{n+1}\right)$. By the construction
of the prediction interval, we have

$Y_{n+1}\in C\left(\tilde{X}_{n+1}\right)\quad\text{ if and only if }\quad
E_{n+1}\leq Q_{1-\alpha}\left(E,\mathcal{I}_{2}\right)\text{, }$

and thus,

${P}\left\\{Y_{n+1}\in
C\left(\tilde{X}_{n+1}\right)\mid\left(\tilde{X}_{i},Y_{i}\right):i\in\mathcal{I}_{1}\right\\}={P}\left\\{E_{n+1}\leq
Q_{1-\alpha}\left(E,\mathcal{I}_{2}\right)\mid\left(\tilde{X}_{i},Y_{i}\right):i\in\mathcal{I}_{1}\right\\}$
(14)

Since the original pairs $\left(\tilde{X}_{i},Y_{i}\right)$ are exchangeable,
so are the calibration variables $E_{i}$ for $i\in\mathcal{I}_{2}$ and
$i=n+1$. Thus, by Lemma 2 on inflated empirical quantiles,

${P}\left\\{E_{n+1}\leq
Q_{1-\alpha}\left(E,\mathcal{I}_{2}\right)\mid\left(\tilde{X}_{i},Y_{i}\right):i\in\mathcal{I}_{1}\right\\}\geq
1-\alpha,$ (15)

and, under the additional assumption that the $E_{i}$ ’s are almost surely
distinct,

${P}\left\\{E_{n+1}\leq
Q_{1-\alpha}\left(E,\mathcal{I}_{2}\right)\mid\left(\tilde{X}_{i},Y_{i}\right):i\in\mathcal{I}_{1}\right\\}\leq
1-\alpha+\frac{1}{\left|\mathcal{I}_{2}\right|+1}$ (16)

The exact coverage result is derived by taking expectations over the proper
training set in 14, 15, and 16. ∎

Next, we likewise prove the validity of the extended CFQP prediction intervals
that control the left and right tails independently.

Corollary 1. Define the prediction interval

$C(\tilde{X}_{n+1}):=\left[\hat{g}_{\alpha_{\mathrm{lo}}}\left(\tilde{X}_{n+1}\right)-Q_{1-\alpha_{\mathrm{lo}}}\left(E_{\mathrm{lo}},\mathcal{I}_{2}\right),\hat{g}_{\alpha_{\mathrm{hi}}}\left(\tilde{X}_{n+1}\right)+Q_{1-\alpha_{\mathrm{hi}}}\left(E_{\mathrm{hi}},\mathcal{I}_{2}\right)\right]$

where $Q_{1-\alpha_{\mathrm{lo}}}\left(E_{\mathrm{lo}},\mathcal{I}_{2}\right)$
is the $\left(1-\alpha_{\mathrm{lo}}\right)$-th empirical quantile of
$\left\\{\hat{g}_{\alpha_{\mathrm{lo}},i}-Y_{i}:i\in\mathcal{I}_{2}\right\\}$
and $Q_{1-\alpha_{\mathrm{hi}}}\left(E_{\mathrm{hi}},\mathcal{I}_{2}\right)$
is the $\left(1-\alpha_{\mathrm{hi}}\right)$-th empirical quantile of
$\left\\{Y_{i}-\hat{g}_{\alpha_{\mathrm{hi}},i}:i\in\mathcal{I}_{2}\right\\}$.
If the samples $\left(\tilde{X}_{n+1},Y_{i}\right),i=1,\ldots,n+1$ are
exchangeable, then

$\displaystyle{P}\left\\{Y_{n+1}\geq\hat{g}_{\alpha_{\mathrm{lo}}}\left(\tilde{X}_{n+1}\right)-Q_{1-\alpha_{\mathrm{lo}}}\left(E_{\mathrm{lo}},\mathcal{I}_{2}\right)\right\\}\geq
1-\alpha_{\mathrm{lo}},$ (17) and
$\displaystyle{P}\left\\{Y_{n+1}\leq\hat{g}_{\alpha_{\mathrm{hi}}}\left(\tilde{X}_{n+1}\right)+Q_{1-\alpha_{\mathrm{hi}}}\left(E_{\mathrm{hi}},\mathcal{I}_{2}\right)\right\\}\geq
1-\alpha_{\mathrm{hi}}.$ (18)

Consequently, we also have ${P}\left\\{Y_{n+1}\in
C\left(\tilde{X}_{n+1}\right)\right\\}\geq 1-\alpha$ assuming
$\alpha=\alpha_{\mathrm{lo}}+\alpha_{\mathrm{hi}}$.

###### Proof.

The two events inside the probabilities (17) as well as (18) are equivalent to
$\hat{g}_{\alpha_{\mathrm{lo}}}\left(X_{n+1}\right)-Y_{n+1}\leq$
$Q_{1-\alpha_{\mathrm{lo}}}\left(E_{\mathrm{lo}},\mathcal{I}_{2}\right)$ and
$Y_{n+1}-\hat{g}_{\alpha_{\mathrm{hi}}}\left(X_{n+1}\right)\leq
Q_{1-\alpha_{\mathrm{hi}}}\left(E_{\mathrm{hi}},\mathcal{I}_{2}\right)$,
respectively. The results are derives by applying Lemma 2 twice, in the same
manner as in the proof of Theorem 1. ∎

##  Appendix C Proofs related to DP

###### Proof of Proposition 2.

By Theorem 2.1(2) in Cheng and Parzen [9], we have

$\sup_{t\in[0,1]}\left|\hat{Q}_{2,q_{\alpha}|s^{\prime}}(t)-{Q}_{q_{\alpha}|s^{\prime}}(t)\right|=O_{p}\left(N^{-1/2}+N^{-\rho}\right)=O_{p}\left(N^{-1/2}\right),$
(19)

for each $k$, as $\rho>1/2$. according to Assumption (A1) , conditions
(Q1)-(Q3) in Cheng and Parzen [9] are satisfied, and since we choose the
kernel function $K_{h}$ with the properties in (A3), conditions (K1)-(K3) in
Cheng and Parzen [9] are ensured. Moreover, the extension from
$[a,b]\subset(0,1)$ to $[0,1]$ is made possible by Assumption (A1). At last,
assumption (A2) and the fact that the bound in equation (2.7) in Cheng and
Parzen [9] is a universal bound which does not depend on $s^{\prime}$ allow
the extension from Eq. 19 to

$\sup_{s^{\prime}}\sup_{t\in[0,1]}\left|\hat{Q}_{2,q_{\alpha}|s^{\prime}}(t)-{Q}_{q_{\alpha}|s^{\prime}}(t)\right|=O_{p}\left(N^{-\nicefrac{{1}}{{2}}}\right),\quad
s^{\prime}=1,\ldots,K.$

∎

Before showing the exact DP guarantee, we utilize Lemma E.1. (stated in Lemma
3) from Chzhen and Schreuder [12] and references therein, where rigorous
proofs are given.

######  Lemma 3.

Let $V_{1},\ldots,V_{n},V_{n+1},n\geq 1$ be exchangeable real-valued random
variables and $U$ distributed uniformly on $[0,1]$ be independent from
$V_{1},\ldots,V_{n},V_{n+1}$, then the constructed location statistic

$T\left(V_{1},\ldots,V_{n},V_{n+1},U\right)=\frac{1}{n+1}\left(\sum_{i=1}^{n}\mathbbm{1}\left\\{V_{i}<V_{n+1}\right\\}+U\cdot\left(1+\sum_{i=1}^{n}\mathbbm{1}\left\\{V_{i}=V_{n+1}\right\\}\right)\right)$

is distributed uniformly on $[0,1]$.

The proof of Theorem 2 is a direct adaptation of quantile demographic parity
guarantee from Chzhen and Schreuder [12] for the mean regression.

###### Proof of Theorem 2..

To prove the claimed DP guarantee for fixed quantile level
$\alpha\in\\{\alpha_{\mathrm{lo}},\alpha_{\mathrm{hi}}\\}$, we will show that
the Kolmogorov-Smirnov distance between $\nu_{\hat{g}_{\alpha}\mid s}$ and
$\nu_{\hat{g}_{\alpha}\mid s^{\prime}}$ equals to zero for any $s\neq
s^{\prime}\in[K]$.

Note that, according to the formulation of $\hat{g}_{\alpha}$ in Eq. (12), we
have for any $(x,s)\in\mathbb{R}^{p}\times[K]$,

$\hat{g}_{\alpha}(x,s)=\sum_{s^{\prime}=1}^{K}\hat{p}_{s^{\prime}}\hat{Q}_{2,q_{\alpha}|s^{\prime}}\circ\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s),\forall\alpha\in\\{\alpha_{\mathrm{lo}},\alpha_{\mathrm{hi}}\\}.$

Denote by
$\hat{Q}(t)=\sum_{s^{\prime}=1}^{K}\hat{p}_{s^{\prime}}\hat{Q}_{2,q_{\alpha}|s^{\prime}}$.
Note that we use the training set to estimate the location statistic for the
new test point, $\hat{Q}(t)$ is independent from
$\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s)$ for each $s\in[K]$.

Since the test point belongs to group $S=s$ for some fixed $s\in[K]$ and, for
all $i=1,\ldots,|\mathcal{I}^{s}_{1}|$, set $V_{i}=\tilde{q}_{1,i}^{s}$ with
$V_{N_{s}+1}\overset{d}{=}(\tilde{q}_{\alpha}(x,s))$ independent from
$\left(V_{i}\right)_{i=1,\ldots,|\mathcal{I}^{s}_{1}|}$. Since the random
variables $V_{1},\ldots,V_{N_{s}},V_{N_{s}+1}$ are exchangeable (more ideally,
independent), Lemma 3 implies that for all $s\in[K]$, the location statistic
$\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s)$ is distributed
uniformly on $[0,1]$. Thus for all $s,s^{\prime}\in[K]$, we have

$\displaystyle\mathrm{KS}\left(\nu_{\hat{g}_{\alpha}\mid
s},\nu_{\hat{g}_{\alpha}\mid s^{\prime}}\right)$ $\displaystyle=$
$\displaystyle\sup_{t\in\mathbb{R}}\left|{P}(\hat{g}_{\alpha}\leq t\mid
S=s)-{P}\left(\hat{g}_{\alpha}\leq t\mid S=s^{\prime}\right)\right|$
$\displaystyle=$
$\displaystyle\sup_{t\in\mathbb{R}}\left|{P}\left(\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s)\leq\hat{Q}^{-1}(t)\mid
S=s\right)-{P}\left(\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s)\leq\hat{Q}^{-1}(t)\mid
S=s^{\prime}\right)\right|$ $\displaystyle=$
$\displaystyle\sup_{t\in\mathbb{R}}\left|{E}\left[\hat{Q}^{-1}(t)\mid
S=s\right]-{E}\left[\hat{Q}^{-1}(t)\mid S=s^{\prime}\right]\right|=0.$

The first equality uses the definition of Eq. (9); the second uses the fact
that $\hat{Q}$ is monotone by construction [12]; finally since the
independence of $\hat{Q}$ is independent from
$\hat{F}_{1,q_{\alpha}|s}\circ\tilde{q}_{\alpha}(x,s)$ conditionally on $S=s$
for any $s\in[K]$, also $\hat{Q}$ remains independent from $S$. The exact DP
is concluded. ∎

##  Appendix D Experiments

###  D.1 Data Description & Pre-processing

  1. 1.

Law School (LAW)[43]: The dataset contains 20,649 examples aiming to predict
students’ GPA based on their information and capacities, with gender as the
sensitive attribute (male vs. female).

  2. 2.

Community&Crime (CRIME)[36]: This dataset contains socio-economic, law
enforcement, and crime data about communities in the US with 1,994 examples.
The task is to predict the number of violent crimes per 100,000 population
with race as the sensitive attribute.

  3. 3.

MEPS 2016 (MEPS)[1, 38]: The Medical Expenditure Panel Survey 2016 dataset
contains information on individuals and their utilization of medical services.
The goal is to predict the health care system utilization score of each
individual by their features including age, marital status, race, poverty
status, health status, health insurance type, and more. There are 15,656
examples on $p=41$ features with race as the sensitive attribute (nonwhite vs.
white).

  4. 4.

Government Salary (GOV)[34, 35]: The government salary dataset (available in R
package "fairadapt") is collected from the 2018 American Community Survey by
the US Census Bureau. The yearly salary for over 200,000 examples is the
response variable, and employee race (7 categories) is identified as the
sensitive attribute.

In order to smoothly run the regression model, several preprocessing steps are
utilized before running the model. For example, in the CRIME dataset, we
impute the missing value by mean, and the race feature is created by the
maximum value of the four races. The clean data are uploaded to the GitHub
repo for future research.

###  D.2 Adaptation of other algorithms in [2, 13]

First, our approach is built upon the one proposed by Chzhen et al. [13], we
incorporated the kernel smoothing procedure in quantile estimation and applied
the local linear smoothing method provided in _NumPy_ [23] which shows its
functionality when there are subgroups of small sample sizes, especially for
the CRIME dataset. It is also feasible to compute Eq. (8) via some global
kernels, such as the gaussian density $K(u)=\exp(-{u^{2}}/{2h^{2}})$ [9, 47,
52] or the triangular density function $K(u)=(1-|u|/h)I(|u|/h\leq 1)$ with a
vanishing bandwidth $h_{n}$. For the bandwidth selection, there is a publicly
available R package _lokern_ with the global bandwidth choice.

Second, we adjusted the reduction-based approach in Agarwal et al. [2] in
several respects: we rescale and discretize the responses, and modify their
algorithm by replacing the loss function $l$ by $\rho_{\alpha}$ of Eq. (4) in
our paper. It is worth mentioning that the reduction-based approach is
sensitive to the hyperparameters: discretization parameter $N$ and slack
$\hat{\varepsilon}_{a}$ in training. We used logistic regression and SVM
classifiers in tuning.

###  D.3 Additional experiment results

We show the additional results using quantile random forest and quantile
neural network for the comparison of our post-processing fairness adjustment
procedure in figures 3 and 4. The three quantile models present consistent
results, our post-processing based upon kernel smoothing outperforms the other
two approaches, which is reflected in more KS value reduction. We also
conducted a brief ablation study by testing the experimental results of
removing either or both of the smoothing strategies in our CFQP approach.
Results are presented in table 2. We found incorporating the jittering and
kernel smoothing methods works better when the subgroups are unbalanced and
there exist subgroups of small sample sizes, especially for the CRIME dataset.
Finally, the computational time analysis is also incorporated in figure 5.

Figure 3:  Results using Quantile random forest for estimating the lower
($\alpha_{\mathrm{lo}}$) and upper ($\alpha_{\mathrm{hi}}$) quantiles using
some state-of-the-art DP-fairness requirement methods on all the datasets.
‘Unfair’, ‘Chzhen’, and ‘Agarwal’ stand for the quantile model without
fairness adjustment, barycenter method [13] and reduction-based algorithm [2]
respectively. We present the MAE and KS of lower quantile estimation, as well
as upper quantile estimation. We set ‘n estimator’ to 50.

Figure 4:  Results using Quantile neural network model for estimating the
lower ($\alpha_{\mathrm{lo}}$) and upper ($\alpha_{\mathrm{hi}}$) quantiles
using some state-of-the-art DP-fairness methods on all the datasets. ‘Unfair’,
‘Chzhen’, and ‘Agarwal’ stand for the quantile model without fairness
adjustment, barycenter method [13] and reduction-based algorithm [2]
respectively. We present the MAE and KS of lower quantile estimation, as well
as upper quantile estimation.

| LAW
---|---
Coverage | Length | KS(lo) | KS(hi)
Ln-CFQP (with both smoothing) |  90.02$\pm$0.51 |  0.46$\pm$.004 |  0.02$\pm$0.01 |  0.02$\pm$0.01
Ln-CFQP (with jittering) |  89.99$\pm$0.49 |  0.39$\pm$.002 |  0.03$\pm$.008 |  0.03$\pm$0.01
Ln-CFQP (with kernel smoothing) |  89.99$\pm$0.50 |  0.38$\pm$.002 |  0.03$\pm$0.009 |  0.03$\pm$0.01
Ln-CFQP (without smoothing) |  89.96$\pm$0.49 |  0.40$\pm$.002 |  0.04$\pm$0.01 |  0.03$\pm$0.01
| CRIME
Coverage | Length | KS(lo) | KS(hi)
Ln-CFQP (with both smoothing) |  90.44$\pm$1.84 |  1.64$\pm$0.05 |  0.11$\pm$0.03 |  0.12$\pm$0.04
Ln-CFQP (with jittering) |  90.55$\pm$1.35 |  1.69$\pm$0.05 |  0.31$\pm$0.11 |  0.33$\pm$0.12
Ln-CFQP (with kernel smoothing) |  90.58$\pm$1.36 |  1.69$\pm$0.05 |  0.25$\pm$0.13 |  0.29$\pm$0.12
Ln-CFQP (without smoothing) |  90.53$\pm$1.27 |  1.69$\pm$0.05 |  0.36$\pm$0.12 |  0.40$\pm$0.12

Table 2: Ablation test results when removing either or both of the smoothing
strategies. Our methods are shown in bold.

In the current experiment, the kernel we used is the local linear one in
defining the quantile functions of subgroups instead of the global kernels.
When calculating the $\tau$-th quantile $x_{\tau}$ of $q_{\alpha}$ using local
linear smoothing, we choose a constant distance size $h$ (kernel radius) and
compute a weighted average for all data points that are closer to $x_{\tau}$
(the closer to $x_{\tau}$ points get higher weights). The time complexity for
computing the local kernel smoother (Eq.(8)) is $\mathcal{O}(1)$, while if we
applied the global kernels in Eq.(8), $\mathcal{O}(n)$ time would cost.

For the CFQP method, The steps determining the time complexity of Algorithms 1
and 2 reside in the following two parts:

  1. 1.

The for-loop where we perform a post-processing which takes
$\sum_{s^{\prime}\in[K]}\mathcal{O}\left(N_{s^{\prime}}\log
N_{s^{\prime}}\right)$ time, as we need to sort the grouped samples;

  2. 2.

The evaluation of $\hat{g}_{\alpha}$ on a new point $(x,s)$ is performed in
$\max_{s^{\prime}\in[K]}\mathcal{O}(\log N_{s})$ time as it involves locating
$\hat{g}_{\alpha}$ in a sorted array.

Figure 5: Empirical running times of methods used for experimental
comparisons. We utilized the linear quantile model on the MEPS dataset. For
methods of "Chzhen" and "Agarwal", normalizing factors are applied to present
the graph.


displaying a rather blue UV/optical slope, most noticeable at $\lambda_{\rm
rest}$ $\lesssim$ 3000Å. The origin of this excess blue light is likely
emission from the accretion disk of the central SMBH, which has been scattered
into our line of sight by gas or dust. Two of the sources in the study of
Assef et al. (2020) are in the Hot DOG sample analyzed in this work
(W0116–0505 and W0220+0137). The contribution of this component to the rest-
frame near-IR emission is, however, small ($\lesssim$ 10% at $H$-band; see
their Fig. 1) and therefore has a limited impact on the estimation of the
stellar masses of Hot DOGs given the significantly larger uncertainties
associated to the SED modeling assumptions.

Figure 16: Results from the SED fitting. The photometric data of the sources
are shown as black open circles. The AGN component is shown as a pink, solid
line, while the E, Sbc and Im galaxy templates are displayed as red, green and
blue dashed lines, respectively. The SED of the host galaxy (E + Sbc + Im) is
shown as a yellow solid line, while the total SED of the source (including the
AGN component) is displayed as a black solid line, with the synthetic
photometry obtained from it (filter-convolved, color-corrected) shown as black
open diamonds. The $H$-band flux of the host, used to obtain its stellar mass
via scaling with a [M/L]H ratio (see text), is shown as a yellowe open square.

After the best SED fit is found, the three galaxy templates are added together
to obtain the rest-frame $H$-band flux density of the host, free from AGN
contamination. In order to estimate the uncertainty associated with this
value, we run 1000 Monte Carlo simulations for each Hot DOG, sampling from the
errors of the photometric data-points assuming they follow a normal
distribution, and calculate the $H$-band flux density for each instance. The
uncertainty is obtained as the standard deviation of the simulated results.
Because the galaxy (and AGN) templates are purely empirical, the SED fitting
procedure does not provide information about the SFR history (SFRH) or any
other intrinsic physical property of the stellar populations of the source.
Therefore, in order to obtain the stellar masses of the host galaxies, an
additional step is required; that is, a [M/L]H ratio. To derive it, we use the
code ezgal (Mancone & Gonzalez, 2012) and the Galaxev stellar population SED
models (Bruzual & Charlot, 2003), which include a contribution from thermally
pulsing-asymptotic giant branch (TP-AGB) stars, with a metallicity Z = 0.008
($\equiv$ 0.4 $Z_{\odot}$, where $Z_{\odot}$ is the metallicity of the Sun)
and a Chabrier intial mass function (IMF) (Chabrier, 2003). We consider
exponentially declining SFRHs that depend only on two parameters: the time of
formation, $t_{\rm f}$, defined with respect to the current age of the galaxy,
$t$, where $\Delta t$ = $t-t_{\rm f}$, and the characteristic decaying time-
scale of the SFR, $\tau$. The SFRH is then parametrized as SFR($t^{\prime})$
$\propto$ e${}^{-(t^{\prime}-t_{\rm f})/\tau}$ with $t^{\prime}$ $\in$
[$t_{\rm f}$, $t$]. The integral of the SFRH provides the $M_{\star}$ of the
model together with the resulting SED, from which the [M/L]H is derived. In
turn, the scaling with the data yields the normalization factor with which to
obtain an absolute value of SFR($t$).

Because it is extremely difficult to infer any information regarding the SFRH
of galaxies based on sparse photometric data, for this work we decide to
consider four simple limiting cases that are selected to represent reasonable
lower and upper boundaries for $\Delta t$ and $\tau$, given the redshifts of
the Hot DOGs, which correspond to between $t$ $\simeq$ 1.3 and 2.2 Gyr after
the Big Bang. Specifically, we assume that star formation began $\Delta t$ =
{500, 100} Myr before the age of each galaxy, $t$, and two time-scales, of
$\tau$ = {0.1, 1} Gyr. Note that the combination of $\Delta t$ = 500 Myr and
$\tau$ = 0.1 Gyr closely mimics a past starburst event, while the combination
of $\Delta t$ = 100 Myr and $\tau$ = 1 Gyr is equivalent to having an almost
constant SFR. The $M_{\star}$ values derived based on these four combination
of parameters are used in Section 6 to place the Hot DOGs along the x-axis in
the star-formation main sequence (Figure 13). In turn, the SFRs obtained from
these SFRHs and scalings are not used in the calculation of the average SFR of
each source in order to avoid the intrinsic correlation existing between SFR
and $M_{\star}$ for each assumed SFRH.

## References

  * Anglés-Alcázar et al. (2017a) Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, F., & Hopkins, P. F. 2017a, MNRAS, 464, 2840, doi: 10.1093/mnras/stw2565
  * Anglés-Alcázar et al. (2017b) Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017b, MNRAS, 472, L109, doi: 10.1093/mnrasl/slx161
  * Angles-Alcazar et al. (2020) Angles-Alcazar, D., Quataert, E., Hopkins, P., et al. 2020, arXiv e-prints, arXiv:2008.12303.  https://arxiv.org/abs/2008.12303
  * Appleton et al. (2013) Appleton, P. N., Guillard, P., Boulanger, F., et al. 2013, ApJ, 777, 66, doi: 10.1088/0004-637X/777/1/66
  * Appleton et al. (2017) Appleton, P. N., Guillard, P., Togi, A., et al. 2017, ApJ, 836, 76, doi: 10.3847/1538-4357/836/1/76
  * Appleton et al. (2018) Appleton, P. N., Diaz-Santos, T., Fadda, D., et al. 2018, ApJ, 869, 61, doi: 10.3847/1538-4357/aaed2a
  * Armus et al. (2009) Armus, L., Mazzarella, J. M., Evans, A. S., et al. 2009, PASP, 121, 559, doi: 10.1086/600092
  * Assef et al. (2010) Assef, R. J., Kochanek, C. S., Brodwin, M., et al. 2010, ApJ, 713, 970, doi: 10.1088/0004-637X/713/2/970
  * Assef et al. (2015) Assef, R. J., Eisenhardt, P. R. M., Stern, D., et al. 2015, ApJ, 804, 27, doi: 10.1088/0004-637X/804/1/27
  * Assef et al. (2020) Assef, R. J., Brightman, M., Walton, D. J., et al. 2020, ApJ, 897, 112, doi: 10.3847/1538-4357/ab9814
  * Bañados et al. (2015) Bañados, E., Decarli, R., Walter, F., et al. 2015, ApJ, 805, L8, doi: 10.1088/2041-8205/805/1/L8
  * Barcos-Muñoz et al. (2017) Barcos-Muñoz, L., Leroy, A. K., Evans, A. S., et al. 2017, ApJ, 843, 117, doi: 10.3847/1538-4357/aa789a
  * Béthermin et al. (2020) Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2, doi: 10.1051/0004-6361/202037649
  * Bischetti et al. (2019) Bischetti, M., Maiolino, R., Carniani, S., et al. 2019, A&A, 630, A59, doi: 10.1051/0004-6361/201833557
  * Bischetti et al. (2017) Bischetti, M., Piconcelli, E., Vietri, G., et al. 2017, A&A, 598, A122, doi: 10.1051/0004-6361/201629301
  * Bischetti et al. (2018) Bischetti, M., Piconcelli, E., Feruglio, C., et al. 2018, A&A, 617, A82, doi: 10.1051/0004-6361/201833249
  * Bischetti et al. (2021) Bischetti, M., Feruglio, C., Piconcelli, E., et al. 2021, A&A, 645, A33, doi: 10.1051/0004-6361/202039057
  * Brauher et al. (2008) Brauher, J. R., Dale, D. A., & Helou, G. 2008, ApJS, 178, 280, doi: 10.1086/590249
  * Brisbin et al. (2015) Brisbin, D., Ferkinhoff, C., Nikola, T., et al. 2015, ApJ, 799, 13, doi: 10.1088/0004-637X/799/1/13
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Capak et al. (2015) Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455, doi: 10.1038/nature14500
  * Carniani et al. (2013) Carniani, S., Marconi, A., Biggs, A., et al. 2013, A&A, 559, A29, doi: 10.1051/0004-6361/201322320
  * Cassata et al. (2020) Cassata, P., Morselli, L., Faisst, A., et al. 2020, A&A, 643, A6, doi: 10.1051/0004-6361/202037517
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Chu et al. (2017) Chu, J. K., Sanders, D. B., Larson, K. L., et al. 2017, ApJS, 229, 25, doi: 10.3847/1538-4365/aa5d15
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11, doi: 10.1111/j.1365-2966.2005.09675.x
  * Cutri et al. (2012) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2012, Explanatory Supplement to the WISE All-Sky Data Release Products, Tech. rep.
  * Davidzon et al. (2017) Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70, doi: 10.1051/0004-6361/201730419
  * Decarli et al. (2016) Decarli, R., Walter, F., Aravena, M., et al. 2016, ApJ, 833, 70, doi: 10.3847/1538-4357/833/1/70
  * Decarli et al. (2018) Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97, doi: 10.3847/1538-4357/aaa5aa
  * Decarli et al. (2019) Decarli, R., Dotti, M., Bañados, E., et al. 2019, ApJ, 880, 157, doi: 10.3847/1538-4357/ab297f
  * Dekel et al. (2009) Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451, doi: 10.1038/nature07648
  * Di Teodoro & Fraternali (2015) Di Teodoro, E. M., & Fraternali, F. 2015, MNRAS, 451, 3021, doi: 10.1093/mnras/stv1213
  * Díaz-Santos et al. (2010) Díaz-Santos, T., Charmandaris, V., Armus, L., et al. 2010, ApJ, 723, 993, doi: 10.1088/0004-637X/723/2/993
  * Díaz-Santos et al. (2011) —. 2011, ApJ, 741, 32, doi: 10.1088/0004-637X/741/1/32
  * Díaz-Santos et al. (2013) Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68, doi: 10.1088/0004-637X/774/1/68
  * Díaz-Santos et al. (2016) Díaz-Santos, T., Assef, R. J., Blain, A. W., et al. 2016, ApJ, 816, L6, doi: 10.3847/2041-8205/816/1/L6
  * Díaz-Santos et al. (2017) Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32, doi: 10.3847/1538-4357/aa81d7
  * Díaz-Santos et al. (2018) Díaz-Santos, T., Assef, R. J., Blain, A. W., et al. 2018, Science, 362, 1034, doi: 10.1126/science.aap7605
  * Eisenhardt et al. (2012) Eisenhardt, P. R. M., Wu, J., Tsai, C.-W., et al. 2012, ApJ, 755, 173, doi: 10.1088/0004-637X/755/2/173
  * Eisenhardt et al. (in prep.) Eisenhardt, P. R. M., et al. in prep.
  * Elvis et al. (1994) Elvis, M., Wilkes, B. J., McDowell, J. C., et al. 1994, ApJS, 95, 1, doi: 10.1086/192093
  * Falgarone et al. (2017) Falgarone, E., Zwaan, M. A., Godard, B., et al. 2017, Nature, 548, 430, doi: 10.1038/nature23298
  * Fan et al. (2017) Fan, L., Jones, S. F., Han, Y., & Knudsen, K. K. 2017, PASP, 129, 124101, doi: 10.1088/1538-3873/aa8e91
  * Fan et al. (2018) Fan, L., Knudsen, K. K., Fogasy, J., & Drouart, G. 2018, ApJ, 856, L5, doi: 10.3847/2041-8213/aab496
  * Farrah et al. (2013) Farrah, D., Lebouteiller, V., Spoon, H. W. W., et al. 2013, ApJ, 776, 38, doi: 10.1088/0004-637X/776/1/38
  * Faucher-Giguère & Quataert (2012) Faucher-Giguère, C.-A., & Quataert, E. 2012, MNRAS, 425, 605, doi: 10.1111/j.1365-2966.2012.21512.x
  * Finnerty et al. (2020) Finnerty, L., Larson, K., Soifer, B. T., et al. 2020, ApJ, 905, 16, doi: 10.3847/1538-4357/abc3bf
  * Förster Schreiber et al. (2006) Förster Schreiber, N. M., Genzel, R., Lehnert, M. D., et al. 2006, ApJ, 645, 1062, doi: 10.1086/504403
  * Förster Schreiber et al. (2009) Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364, doi: 10.1088/0004-637X/706/2/1364
  * Fujimoto et al. (2020) Fujimoto, S., Silverman, J. D., Bethermin, M., et al. 2020, ApJ, 900, 1, doi: 10.3847/1538-4357/ab94b3
  * Fumagalli et al. (2011) Fumagalli, M., Prochaska, J. X., Kasen, D., et al. 2011, MNRAS, 418, 1796, doi: 10.1111/j.1365-2966.2011.19599.x
  * Gaspari et al. (2013) Gaspari, M., Ruszkowski, M., & Oh, S. P. 2013, MNRAS, 432, 3401, doi: 10.1093/mnras/stt692
  * Gillman et al. (2019) Gillman, S., Swinbank, A. M., Tiley, A. L., et al. 2019, MNRAS, 486, 175, doi: 10.1093/mnras/stz765
  * Gnerucci et al. (2011) Gnerucci, A., Marconi, A., Cresci, G., et al. 2011, A&A, 528, A88, doi: 10.1051/0004-6361/201015465
  * González-López et al. (in prep.) González-López, J., et al. in prep.
  * Graciá-Carpio et al. (2011) Graciá-Carpio, J., Sturm, E., Hailey-Dunsheath, S., et al. 2011, ApJ, 728, L7, doi: 10.1088/2041-8205/728/1/L7
  * Griffith et al. (2012) Griffith, R. L., Kirkpatrick, J. D., Eisenhardt, P. R. M., et al. 2012, AJ, 144, 148, doi: 10.1088/0004-6256/144/5/148
  * Guillard et al. (2015) Guillard, P., Boulanger, F., Lehnert, M. D., et al. 2015, A&A, 574, A32, doi: 10.1051/0004-6361/201423612
  * Haan et al. (2011) Haan, S., Surace, J. A., Armus, L., et al. 2011, AJ, 141, 100, doi: 10.1088/0004-6256/141/3/100
  * Harikane et al. (2018) Harikane, Y., Ouchi, M., Shibuya, T., et al. 2018, ApJ, 859, 84, doi: 10.3847/1538-4357/aabd80
  * Häring & Rix (2004) Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89, doi: 10.1086/383567
  * Herrero-Illana et al. (2019) Herrero-Illana, R., Privon, G. C., Evans, A. S., et al. 2019, A&A, 628, A71, doi: 10.1051/0004-6361/201834088
  * Hodge et al. (2019) Hodge, J. A., Smail, I., Walter, F., et al. 2019, ApJ, 876, 130, doi: 10.3847/1538-4357/ab1846
  * Hopkins (2012) Hopkins, P. F. 2012, MNRAS, 420, L8, doi: 10.1111/j.1745-3933.2011.01179.x
  * Hopkins et al. (2008) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356, doi: 10.1086/524362
  * Hopkins et al. (2005) Hopkins, P. F., Hernquist, L., Martini, P., et al. 2005, ApJ, 625, L71, doi: 10.1086/431146
  * Howell et al. (2010) Howell, J. H., Armus, L., Mazzarella, J. M., et al. 2010, ApJ, 715, 572, doi: 10.1088/0004-637X/715/1/572
  * Inami et al. (2018) Inami, H., Armus, L., Matsuhara, H., et al. 2018, A&A, 617, A130, doi: 10.1051/0004-6361/201833053
  * Iwasawa et al. (2011) Iwasawa, K., Sanders, D. B., Teng, S. H., et al. 2011, A&A, 529, A106, doi: 10.1051/0004-6361/201015264
  * Izumi et al. (2019) Izumi, T., Onoue, M., Matsuoka, Y., et al. 2019, PASJ, 71, 111, doi: 10.1093/pasj/psz096
  * Jones et al. (2020) Jones, G. C., Béthermin, M., Fudamoto, Y., et al. 2020, MNRAS, 491, L18, doi: 10.1093/mnrasl/slz154
  * Jones et al. (2014) Jones, S. F., Blain, A. W., Stern, D., et al. 2014, MNRAS, 443, 146, doi: 10.1093/mnras/stu1157
  * Jones et al. (2017) Jones, S. F., Blain, A. W., Assef, R. J., et al. 2017, MNRAS, 469, 4565, doi: 10.1093/mnras/stx1141
  * Jun et al. (2020a) Jun, H. D., Assef, R. J., Carroll, C. M., et al. 2020a, arXiv e-prints, arXiv:2010.15460.  https://arxiv.org/abs/2010.15460
  * Jun et al. (2020b) Jun, H. D., Assef, R. J., Bauer, F. E., et al. 2020b, ApJ, 888, 110, doi: 10.3847/1538-4357/ab5e7b
  * Kim et al. (2013) Kim, D.-C., Evans, A. S., Vavilkin, T., et al. 2013, ApJ, 768, 102, doi: 10.1088/0004-637X/768/2/102
  * Kimball et al. (2015) Kimball, A. E., Lacy, M., Lonsdale, C. J., & Macquart, J.-P. 2015, MNRAS, 452, 88, doi: 10.1093/mnras/stv1160
  * Knudsen et al. (2016) Knudsen, K. K., Richard, J., Kneib, J.-P., et al. 2016, MNRAS, 462, L6, doi: 10.1093/mnrasl/slw114
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * LaMassa et al. (2017) LaMassa, S. M., Glikman, E., Brusa, M., et al. 2017, ApJ, 847, 100, doi: 10.3847/1538-4357/aa87b5
  * Langer & Pineda (2015) Langer, W. D., & Pineda, J. L. 2015, A&A, 580, A5, doi: 10.1051/0004-6361/201525950
  * Le Fèvre et al. (2020) Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1, doi: 10.1051/0004-6361/201936965
  * Leung et al. (2019) Leung, T. K. D., Riechers, D. A., Baker, A. J., et al. 2019, ApJ, 871, 85, doi: 10.3847/1538-4357/aaf860
  * Linden et al. (2017) Linden, S. T., Evans, A. S., Rich, J., et al. 2017, ArXiv e-prints.  https://arxiv.org/abs/1705.03370
  * Lupi et al. (2019) Lupi, A., Volonteri, M., Decarli, R., et al. 2019, MNRAS, 488, 4004, doi: 10.1093/mnras/stz1959
  * Lutz et al. (2016) Lutz, D., Berta, S., Contursi, A., et al. 2016, A&A, 591, A136, doi: 10.1051/0004-6361/201527706
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Maiolino et al. (2005) Maiolino, R., Cox, P., Caselli, P., et al. 2005, A&A, 440, L51, doi: 10.1051/0004-6361:200500165
  * Maiolino et al. (2015) Maiolino, R., Carniani, S., Fontana, A., et al. 2015, MNRAS, 452, 54, doi: 10.1093/mnras/stv1194
  * Malhotra et al. (2001) Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766, doi: 10.1086/323046
  * Mancone & Gonzalez (2012) Mancone, C. L., & Gonzalez, A. H. 2012, PASP, 124, 606, doi: 10.1086/666502
  * Mazzucchelli et al. (2017) Mazzucchelli, C., Bañados, E., Venemans, B. P., et al. 2017, ApJ, 849, 91, doi: 10.3847/1538-4357/aa9185
  * McConnell & Ma (2013) McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184, doi: 10.1088/0004-637X/764/2/184
  * McGreer et al. (2014) McGreer, I. D., Fan, X., Strauss, M. A., et al. 2014, AJ, 148, 73, doi: 10.1088/0004-6256/148/4/73
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Miller et al. (2018) Miller, T. B., Chapman, S. C., Aravena, M., et al. 2018, Nature, 556, 469, doi: 10.1038/s41586-018-0025-2
  * Mullaney et al. (2011) Mullaney, J. R., Alexander, D. M., Goulding, A. D., & Hickox, R. C. 2011, MNRAS, 414, 1082, doi: 10.1111/j.1365-2966.2011.18448.x
  * Murphy et al. (2009) Murphy, E. J., Chary, R., Alexander, D. M., et al. 2009, ApJ, 698, 1380, doi: 10.1088/0004-637X/698/2/1380
  * Murphy et al. (2011) Murphy, E. J., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 732, 126, doi: 10.1088/0004-637X/732/2/126
  * Muzzin et al. (2012) Muzzin, A., Wilson, G., Yee, H. K. C., et al. 2012, ApJ, 746, 188, doi: 10.1088/0004-637X/746/2/188
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Neeleman et al. (2019) Neeleman, M., Bañados, E., Walter, F., et al. 2019, ApJ, 882, 10, doi: 10.3847/1538-4357/ab2ed3
  * Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 160, doi: 10.1086/590483
  * Nguyen et al. (2020) Nguyen, N. H., Lira, P., Trakhtenbrot, B., et al. 2020, arXiv e-prints, arXiv:2003.00525.  https://arxiv.org/abs/2003.00525
  * Ota et al. (2014) Ota, K., Walter, F., Ohta, K., et al. 2014, ApJ, 792, 34, doi: 10.1088/0004-637X/792/1/34
  * Ouchi et al. (2013) Ouchi, M., Ellis, R., Ono, Y., et al. 2013, ApJ, 778, 102, doi: 10.1088/0004-637X/778/2/102
  * Penney et al. (2019) Penney, J. I., Blain, A. W., Wylezalek, D., et al. 2019, MNRAS, 483, 514, doi: 10.1093/mnras/sty3128
  * Penney et al. (2020) Penney, J. I., Blain, A. W., Assef, R. J., et al. 2020, MNRAS, 496, 1565, doi: 10.1093/mnras/staa1582
  * Perrotta et al. (2019) Perrotta, S., Hamann, F., Zakamska, N. L., et al. 2019, MNRAS, 488, 4126, doi: 10.1093/mnras/stz1993
  * Petitjean et al. (1996) Petitjean, P., Pécontal, E., Valls-Gabaud, D., & Chariot, S. 1996, Nature, 380, 411, doi: 10.1038/380411a0
  * Price et al. (2019) Price, S. H., Kriek, M., Barro, G., et al. 2019, arXiv e-prints, arXiv:1902.09554.  https://arxiv.org/abs/1902.09554
  * Rafelski et al. (2012) Rafelski, M., Wolfe, A. M., Prochaska, J. X., Neeleman, M., & Mendez, A. J. 2012, ApJ, 755, 89, doi: 10.1088/0004-637X/755/2/89
  * Rich et al. (2015) Rich, J. A., Kewley, L. J., & Dopita, M. A. 2015, ApJS, 221, 28, doi: 10.1088/0067-0049/221/2/28
  * Riguccini et al. (2015) Riguccini, L., Le Floc’h, E., Mullaney, J. R., et al. 2015, MNRAS, 452, 470, doi: 10.1093/mnras/stv1297
  * Riguccini et al. (2019) Riguccini, L. A., Treister, E., Menéndez-Delmestre, K., et al. 2019, AJ, 157, 233, doi: 10.3847/1538-3881/ab16cd
  * Sanders et al. (1988) Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988, ApJ, 325, 74, doi: 10.1086/165983
  * Schawinski et al. (2015) Schawinski, K., Koss, M., Berney, S., & Sartori, L. F. 2015, MNRAS, 451, 2517, doi: 10.1093/mnras/stv1136
  * Schreiber et al. (2015) Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74, doi: 10.1051/0004-6361/201425017
  * Scoville et al. (2015) Scoville, N., Sheth, K., Walter, F., et al. 2015, ApJ, 800, 70, doi: 10.1088/0004-637X/800/1/70
  * Scoville et al. (2017) Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150, doi: 10.3847/1538-4357/aa61a0
  * Smirnova-Pinchukova et al. (2019) Smirnova-Pinchukova, I., Husemann, B., Busch, G., et al. 2019, A&A, 626, L3, doi: 10.1051/0004-6361/201935577
  * Smith et al. (2017) Smith, J. D. T., Croxall, K., Draine, B., et al. 2017, ApJ, 834, 5, doi: 10.3847/1538-4357/834/1/5
  * Stacey et al. (2010) Stacey, G. J., Hailey-Dunsheath, S., Ferkinhoff, C., et al. 2010, ApJ, 724, 957, doi: 10.1088/0004-637X/724/2/957
  * Stern (2015) Stern, D. 2015, ApJ, 807, 129, doi: 10.1088/0004-637X/807/2/129
  * Stern et al. (2014) Stern, D., Lansbury, G. B., Assef, R. J., et al. 2014, ApJ, 794, 102, doi: 10.1088/0004-637X/794/2/102
  * Stierwalt et al. (2013) Stierwalt, S., Armus, L., Surace, J. A., et al. 2013, ApJS, 206, 1, doi: 10.1088/0067-0049/206/1/1
  * Stierwalt et al. (2014) Stierwalt, S., Armus, L., Charmandaris, V., et al. 2014, ApJ, 790, 124, doi: 10.1088/0004-637X/790/2/124
  * Stott et al. (2016) Stott, J. P., Swinbank, A. M., Johnson, H. L., et al. 2016, MNRAS, 457, 1888, doi: 10.1093/mnras/stw129
  * Tacconi et al. (2018) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179, doi: 10.3847/1538-4357/aaa4b4
  * Tadaki et al. (2018) Tadaki, K., Iono, D., Yun, M. S., et al. 2018, Nature, 560, 613, doi: 10.1038/s41586-018-0443-1
  * Tan et al. (2019) Tan, Q.-H., Gao, Y., Kohno, K., et al. 2019, ApJ, 887, 24, doi: 10.3847/1538-4357/ab50be
  * Temple et al. (2019) Temple, M. J., Banerji, M., Hewett, P. C., et al. 2019, MNRAS, 487, 2594, doi: 10.1093/mnras/stz1420
  * Torres-Albà et al. (2018) Torres-Albà, N., Iwasawa, K., Díaz-Santos, T., et al. 2018, A&A, 620, A140, doi: 10.1051/0004-6361/201834105
  * Trainor & Steidel (2013) Trainor, R., & Steidel, C. C. 2013, ApJ, 775, L3, doi: 10.1088/2041-8205/775/1/L3
  * Trakhtenbrot et al. (2017) Trakhtenbrot, B., Lira, P., Netzer, H., et al. 2017, ApJ, 836, 8, doi: 10.3847/1538-4357/836/1/8
  * Tsai et al. (2015) Tsai, C.-W., Eisenhardt, P. R. M., Wu, J., et al. 2015, ApJ, 805, 90, doi: 10.1088/0004-637X/805/2/90
  * Tsai et al. (2018) Tsai, C.-W., Eisenhardt, P. R. M., Jun, H. D., et al. 2018, ApJ, 868, 15, doi: 10.3847/1538-4357/aae698
  * Veilleux et al. (2002) Veilleux, S., Kim, D., & Sanders, D. B. 2002, ApJS, 143, 315, doi: 10.1086/343844
  * Veilleux et al. (2009) Veilleux, S., Rupke, D. S. N., Kim, D., et al. 2009, ApJS, 182, 628, doi: 10.1088/0067-0049/182/2/628
  * Venemans et al. (2016) Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37, doi: 10.3847/0004-637X/816/1/37
  * Venemans et al. (2012) Venemans, B. P., McMahon, R. G., Walter, F., et al. 2012, ApJ, 751, L25, doi: 10.1088/2041-8205/751/2/L25
  * Venemans et al. (2017) Venemans, B. P., Walter, F., Decarli, R., et al. 2017, ApJ, 837, 146, doi: 10.3847/1538-4357/aa62ac
  * Venemans et al. (2018) Venemans, B. P., Decarli, R., Walter, F., et al. 2018, ApJ, 866, 159, doi: 10.3847/1538-4357/aadf35
  * Venemans et al. (2020) Venemans, B. P., Walter, F., Neeleman, M., et al. 2020, ApJ, 904, 130, doi: 10.3847/1538-4357/abc563
  * Vietri et al. (2018) Vietri, G., Piconcelli, E., Bischetti, M., et al. 2018, A&A, 617, A81, doi: 10.1051/0004-6361/201732335
  * Walter et al. (2009) Walter, F., Riechers, D., Cox, P., et al. 2009, Nature, 457, 699, doi: 10.1038/nature07681
  * Walter et al. (2020) Walter, F., Carilli, C., Neeleman, M., et al. 2020, arXiv e-prints, arXiv:2009.11126.  https://arxiv.org/abs/2009.11126
  * Wang et al. (2009) Wang, J.-G., Dong, X.-B., Wang, T.-G., et al. 2009, ApJ, 707, 1334, doi: 10.1088/0004-637X/707/2/1334
  * Wang et al. (2013) Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44, doi: 10.1088/0004-637X/773/1/44
  * Wang et al. (2016) Wang, R., Wu, X.-B., Neri, R., et al. 2016, ApJ, 830, 53, doi: 10.3847/0004-637X/830/1/53
  * Whiting (2012) Whiting, M. T. 2012, MNRAS, 421, 3242, doi: 10.1111/j.1365-2966.2012.20548.x
  * Willott et al. (2015) Willott, C. J., Bergeron, J., & Omont, A. 2015, ApJ, 801, 123, doi: 10.1088/0004-637X/801/2/123
  * Willott et al. (2017) —. 2017, ApJ, 850, 108, doi: 10.3847/1538-4357/aa921b
  * Willott et al. (2013) Willott, C. J., Omont, A., & Bergeron, J. 2013, ApJ, 770, 13, doi: 10.1088/0004-637X/770/1/13
  * Wisnioski et al. (2015) Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209, doi: 10.1088/0004-637X/799/2/209
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Wu et al. (2012) Wu, J., Tsai, C.-W., Sayers, J., et al. 2012, ApJ, 756, 96, doi: 10.1088/0004-637X/756/1/96
  * Wu et al. (2014) Wu, J., Bussmann, R. S., Tsai, C.-W., et al. 2014, ApJ, 793, 8, doi: 10.1088/0004-637X/793/1/8
  * Wu et al. (2018) Wu, J., Jun, H. D., Assef, R. J., et al. 2018, ApJ, 852, 96, doi: 10.3847/1538-4357/aa9ff3
  * Zakamska et al. (2016) Zakamska, N. L., Hamann, F., Pâris, I., et al. 2016, MNRAS, 459, 3144, doi: 10.1093/mnras/stw718
  * Zanella et al. (2018) Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976, doi: 10.1093/mnras/sty2394