# Differentially Private Boxplots

Kelly Ramsay
Department of Mathematics and Statistics
York University
Toronto, Ontario M3J 1P3
<EMAIL_ADDRESS>
&Jairo Diaz-Rodriguez
Department of Mathematics and Statistics
York University
Toronto, Ontario M3J 1P3
<EMAIL_ADDRESS>

###### Abstract

Despite the potential of differentially private data visualization to
harmonize data analysis and privacy, research in this area remains relatively
underdeveloped. Boxplots are a widely popular visualization used for
summarizing a dataset and for comparison of multiple datasets.
Consequentially, we introduce a differentially private boxplot. We evaluate
its effectiveness for displaying location, scale, skewness and tails of a
given empirical distribution. In our theoretical exposition, we show that the
location and scale of the boxplot are estimated with optimal sample
complexity, and the skewness and tails are estimated consistently. In
simulations, we show that this boxplot performs similarly to a non-private
boxplot, and it outperforms a boxplot naively constructed from existing
differentially private quantile algorithms. Additionally, we conduct a real
data analysis of Airbnb listings, which shows that comparable analysis can be
achieved through differentially private boxplot visualization.

##  1 Introduction

It is now well established that differential privacy (Dwork et al.,, 2006) is
a powerful framework for protecting the privacy of individuals’ data. As a
result, a plethora of differentially private data analysis tools have been
developed over the last several decades (Dwork,, 2008; Ji et al.,, 2014; Liu
et al.,, 2023). However, one area that has been considerably underdeveloped is
that of differentially private visualization. This is despite the fact that
data visualization is a key tool in exploratory analysis, which is an
essential component of data analysis. In fact, in a recent study of data
analysts and practitioners, Garrido et al., (2023) found that “most analysts
employed aggregations and visualizations to fulfill their use case,” rather
than machine learning models.

Differentially private versions of several popular visualizations have been
considered in the literature. There has been substantial study of the
differentially private histogram (Hay et al.,, 2010; Li et al.,, 2010; Acs et
al.,, 2012; Kellaris and Papadopoulos,, 2013; Xu et al.,, 2013; Zhang et al.,,
2014; Budiu et al.,, 2022). Visualizations based on clustering (Hongde et
al.,, 2014), for mobility data (He et al.,, 2016), heatmaps (Zhang et al.,,
2016) and scatterplots (Panavas et al.,, 2024) have also been considered.

Surprisingly, despite its widespread use, a differentially private boxplot has
not been directly studied. Boxplots, invented by Tukey et al., (1977), are
used for visualizing the characteristics of a univariate sample, as well as
for visualizing the relationship between a continuous variable and a
categorical variable. Despite its simplicity, many key distributional
characteristics can be evaluated from the boxplot: namely, location, scale,
skewness, and tails. This has made it popular among practitioners.

In addition, its simplicity is beneficial from a privacy standpoint. The
boxplot only requires estimating a few summary statistics in order to convey
the distributional characteristics of a univariate sample, making the boxplot
efficient in its use of privacy budget. For instance, by contrast, a histogram
can also be used to convey the same information. However, a histogram is
generally a consistent estimate of the population density, which intuitively,
should then require more noise to ensure privacy.

Motivated by these observations, we develop and study a differentially private
boxplot. Our contributions are as follows: We introduce a differentially
private boxplot and evaluate its ability to convey location, scale, skewness
and tails both theoretically and empirically. We show that the inner quantiles
used for the location and scale of the boxplot are estimated optimally
(Theorem 2). We then show that tails and skewness are estimated consistently
(Theorem 5). Through extensive simulations, we show that the proposed boxplot
is often more accurate than one naively constructed via popular, existing
differentially private quantile methods (see Section 5). We also conduct a
differentially private exploratory data analysis, demonstrating the utility of
the proposed boxplots, as well as their limitations (see Section 6). As a by-
product of our work, we develop some new results concerning private quantiles.
We introduce a minimax lower bound on private quantile estimation, and show
that when the number of quantiles is small, the quantile algorithms JointExp
(Gillenwater et al.,, 2021) and PrivateQuantile (Smith,, 2011) attain this
lower bound, up to logarithmic factors. We also theoretically confirm the
empirical observation of Durfee, (2023), which says that the unbounded
quantile algorithm performs better than existing algorithms for estimating
extreme quantiles. We show that JointExp is inconsistent for extreme quantiles
(Lemma 9), where the unbounded algorithm of Durfee, (2023) is consistent
(Lemma 8). We also present simulations comparing four differentially private
quantile algorithms, those of Smith, (2011); Gillenwater et al., (2021);
Kaplan et al., (2022) and Durfee, (2023).

Before turning to our main results, we mention some other related works.
Nanayakkara et al., (2022) developed a plot to visualize privacy utility
trade-offs. Zhou et al., (2022) developed DPVisCreator, a visualization
system, which relies on publishing a synthetic dataset. Visualizations in
other privacy models were considered by Dasgupta and Kosara, (2011); Dasgupta
et al., (2013, 2019), and recently reviewed by Bhattacharjee et al., (2020).
In addition, a related area is that of differentially private quantile
estimation (Smith,, 2011; Gillenwater et al.,, 2021; Kaplan et al.,, 2022;
Alabi et al.,, 2022; Durfee,, 2023; Lalanne et al., 2023a, ; Lalanne et al.,
2023b, ). One could apply any one of these algorithms to create a
differentially private boxplot. We show in Section 5 that this approach is
often less accurate than the proposed methodology.

##  2 Preliminaries

First, we briefly review the boxplot. We define a dataset of size
$n\in\mathbb{N}$ to be a set of real numbers of size $n$. Let
$\mathcal{D}_{n}$ be the set of datasets of size $n$ and let
$\mathbf{X}_{n}=\\{X_{1},\ldots,X_{n}\\}$ be a dataset size $n$. The boxplot
consists of a box with a line drawn through it, with two whiskers emanating
from the lower and upper bounds of the box. Specific points are sometimes
indicated above or below the whiskers. The central line inside the box marks
the median of the dataset. The box itself is constructed from the quartiles,
detailing the middle 50% of $\mathbf{X}_{n}$. The lower (upper) whisker is
constructed from the larger (smaller) of the minimum (maximum) of
$\mathbf{X}_{n}$ and the lower (upper) quartile of $\mathbf{X}_{n}$ minus
(plus) 1.5 times the interquartile range of $\mathbf{X}_{n}$. Lastly, any
points in $\mathbf{X}_{n}$ falling outside of the whiskers are added to the
plot. See Figure 2 for an example. The boxplot describes various
characteristics of $\mathbf{X}_{n}$. The median line approximates the
_location_ of $\mathbf{X}_{n}$. The box itself approximates the spread or
_scale_ of $\mathbf{X}_{n}$. The _skewness_ of $\mathbf{X}_{n}$ can be
approximated as follows: A median line placed in the center of the box, along
with whiskers that are similar in magnitude, suggests a symmetric
distribution. A median line placed away from the center of the box or
asymmetric whiskers indicates skewness. Additionally, the _tails_ can be
assessed by the number of points beyond the whiskers. Many such points signify
outliers or heavy tails, highlighting data points that deviate significantly
from the rest.

Next, we introduce differential privacy. First, we say that
$\mathbf{Y}_{n}\in\mathcal{D}_{n}$, is adjacent to $\mathbf{X}_{n}$ if
$\mathbf{X}_{n}$ and $\mathbf{Y}_{n}$ differ in exactly one point. Let
$\mathcal{A}_{n}$ be the set of pairs of adjacent datasets of size $n$. Next,
denote by $P_{\mathbf{X}_{n}}$ a probability measure which depends on
$\mathbf{X}_{n}$. For a space $S$, let $\mathscr{B}(S)$ denote the Borel sets
of $S$ and let $\mathcal{M}_{1}(S)$ be the set of Borel probability measures
over $S$. Formally, $P_{\mathbf{X}_{n}}$ is a map from $\mathbb{R}^{n}$ to
$\mathcal{M}_{1}(\mathbb{R}^{d})$, for some $d\in\mathbb{N}$. We can now
define differential privacy.

######  Definition 1.

The quantity $\theta\sim P_{\mathbf{X}_{n}}$ is $\epsilon$-differentially
private if

$\sup_{(\mathbf{X}_{n},\mathbf{Y}_{n})\in\mathcal{A}_{n}}\sup_{B\in\mathscr{B}(\mathbb{R}^{d})}\frac{P_{\mathbf{X}_{n}}(B)}{P_{\mathbf{Y}_{n}}(B)}\leq
e^{\epsilon}.$ (1)

Here, $\theta$ is a $d$-dimensional differentially private quantity and
$\epsilon$ is the privacy budget, where values of $\epsilon$ which are closer
to zero enforce stricter privacy constraints.

The proposed differentially private boxplot makes use of existing
differentially private algorithms. Specifically, the Laplace mechanism, and
algorithms for computing differentially private quantiles. The Laplace
mechanism (Dwork et al.,, 2006) is a fundamental mechanism for constructing
differentially private statistics. Let $Z$ be a standard Laplace random
variable. Specifically, for a statistic $T\colon\mathcal{D}_{n}\to\mathbb{R}$,
define its global sensitivity to be
$\sup_{(\mathbf{X}_{n},\mathbf{Y}_{n})\in\mathcal{A}_{n}}|T(\mathbf{X}_{n})-T(\mathbf{Y}_{n})|$.
If $T$ has global sensitivity bounded by 1, then
$\tilde{T}=T(\mathbf{X}_{n})+Z/\epsilon$ is a differentially private quantity.
This is known as the Laplace mechanism.

Of course, a private boxplot relies on differentially private quantiles. The
proposed private boxplot relies on the unbounded algorithm of Durfee, (2023)
and the JointExp algorithm of Gillenwater et al., (2021). These algorithms are
the main focus of our theoretical exposition. We do, however, consider the
performance of the naive boxplot constructed via the algorithms of Smith,
(2011), Kaplan et al., (2022) and Durfee, (2023) in our simulation study, see
Section 5. We also discuss the extent to which our theoretical results extend
to these other algorithms, see Remark 4.

We now give a brief overview of JointExp and unbounded. JointExp works as
follows: Suppose we wish to estimate $m\in\mathbb{N}$ quantiles of levels
$0<q_{1}<\ldots<q_{m}<1$. For an integer $k\in\mathbb{N}$, let $[k]$ denote
the set $\\{1,\ldots,k\\}$. Define the $j$th quantile generated by JointExp
for a given $j\in[m]$ to be $\tilde{\xi}_{q_{j}}$, and let
$\tilde{\xi}_{q}=(\tilde{\xi}_{q_{1}},\ldots,\tilde{\xi}_{q_{m}})$. For a
given $\nu\in\mathcal{M}_{1}(\mathbb{R})$, let $F_{\nu}$ be the cumulative
distribution function of $\nu$ and, if it exists, let $f_{\nu}$ be its
probability density function. Now, define the following utility function:
$\phi(x)=-\sum_{j=1}^{m+1}|F_{\hat{\nu}}(x_{j})-F_{\hat{\nu}}(x_{j-1})-(q_{j}-q_{j-1})|,$
where $x_{0}=-\infty<x_{1}<\ldots<x_{m}<x_{m+1}=\infty$, $q_{0}=0$ and
$q_{m+1}=1$. Then $\tilde{\xi}_{q}$ has law $Q_{\mathbf{X}_{n}}$, where
$f_{Q_{\mathbf{X}_{n}}}\propto\exp(-n\phi(x)/2\epsilon)\mathbbm{1}\left\\{x\in[a,b]\right\\}.$
It was shown by Gillenwater et al., (2021) that $\tilde{\xi}_{q}$ is
$\epsilon$-differentially private. Furthermore, Lalanne et al., 2023b  showed
that this algorithm is consistent for continuous distributions. If one
suspects that the population distribution contains atoms, then we recommend
that one uses the jittering modification proposed by Lalanne et al., 2023b .
We show in Section 4 that this algorithm has optimal sample complexity for
estimating a single quantile, and therefore, for the scale and location of the
boxplot. (Here we are only interested in estimating three quantiles, and so,
how sample complexity scales in $m$ in not of particular concern.)

The unbounded algorithm produces a private quantile given a lower bound on the
data. We present a slightly modified version, which performs better for
extreme quantiles. For a given quantile $q\in(1/2,1)$ and $a>0$, let
$V_{0},V_{1},\ldots$ be independent, standard exponential random variables and
let $\beta>1$. Define $i^{*}=\inf\\{i\in\mathbb{N}\colon
F_{\hat{\nu}}(\beta^{i}-a-1)+\frac{2}{n\epsilon}V_{i}\geq
q+V_{0}\frac{2}{n\epsilon}\\}$. The quantile estimate generated by the
unbounded algorithm is given by: $\tilde{\psi}_{q}=\beta^{i^{*}}-a-1$. If
$q\in(0,1/2)$, then, we apply the above procedure to the dataset
$-\mathbf{X}_{n}=\\{-X_{1},\ldots,-X_{n}\\}$, with input parameters $1-q$ and
$a=-b$. It follows from (Durfee,, 2023) that this algorithm is
$\epsilon$-differentially private. In the proposed differentially private
boxplot algorithm, we use unbounded to estimate the minimum and maximum of the
dataset. For a measure $\nu\in\mathcal{M}_{1}(\mathbb{R})$, let
$\xi_{p,\nu}=\inf\\{x\in\mathbb{R}\colon p\leq F_{\nu}(x)\\}$ be the
associated $p$th quantile, $\min(\nu)=\inf\\{x\colon F_{\nu}(x)>0\\}$ and
$\max(\nu)=\inf\\{x\colon F_{\nu}(x)=1\\}$. We show that if $\mathbf{X}_{n}$
is an independent sample from $\nu$, then when $q\to 1$ as $n\to\infty$,
$\tilde{\psi}_{q}$ is weakly consistent for $\max(\nu)$, and when $q\to 0$ as
$n\to\infty$, $\tilde{\psi}_{q}$ is weakly consistent for $\min(\nu)$, see
Lemma 8. By contrast, we show that JointExp is inconsistent for $q=O(1/n)$,
see Lemma 9. Instead of exponential noise, one could also use Laplace or
Gumbel noise (Durfee,, 2023). We found this to make little difference in
simulation, and so we only present the version which incorporates exponential
noise.

##  3 A differentially private boxplot

We can now introduce the differentially private boxplot. Let $q=(1/4,1/2,3/4)$
and $\tilde{\operatorname{IQR}}=\tilde{\xi}_{3/4}-\tilde{\xi}_{1/4}$. Given
lower and upper bounds on the data $a$ and $b$111If one does not wish to
supply input bounds for the data, one can use the “fully unbounded” version of
unbounded (Durfee,, 2023). However, in simulation, (see Section 5), we have
observed that the procedure is still accurate, even when the input bounds are
very loose., we first generate the minimum and maximum estimates using the
unbounded algorithm, denoted by $\tilde{\psi}_{c/\sqrt{n}}$ and
$\tilde{\psi}_{1-c/\sqrt{n}}$. Next, we run JointExp with input bounds
$(\tilde{\psi}_{c/\sqrt{n}},\tilde{\psi}_{1-c/\sqrt{n}})$ to generate
$\tilde{\xi}_{q}$. The box bounds and center line are then constructed from
$(\tilde{\xi}_{1/4},\tilde{\xi}_{1/2},\tilde{\xi}_{3/4})$. Letting
$\tilde{\ell}_{1}=\tilde{\xi}_{1/4}-\tilde{\operatorname{IQR}}$,
$\tilde{u}_{1}=\tilde{\xi}_{3/4}+\tilde{\operatorname{IQR}}$ and $\lambda_{n}$
be a sequence decreasing to 0, the lower and upper whisker are defined as

$\tilde{\ell}=\begin{cases}\tilde{\ell}_{1}&\tilde{\xi}_{c/\sqrt{n}}\leq\lambda_{n}|\tilde{\ell}_{1}|+\tilde{\ell}_{1}\\\
\tilde{\psi}_{c/\sqrt{n}}&\tilde{\xi}_{c/\sqrt{n}}>\lambda_{n}|\tilde{\ell}_{1}|+\tilde{\ell}_{1}\end{cases}\qquad\text{and}\qquad\tilde{u}=\begin{cases}\tilde{u}_{1}&\tilde{\xi}_{1-c/\sqrt{n}}\geq\tilde{u}_{1}-\lambda_{n}|\tilde{u}_{1}|\\\
\tilde{\psi}_{1-c/\sqrt{n}}&\tilde{\xi}_{1-c/\sqrt{n}}<\tilde{u}_{1}-\lambda_{n}|\tilde{u}_{1}|\end{cases},$

respectively. To elaborate, the role of the minimum and the maximum of the
dataset in a traditional boxplot is played by the estimated extreme quantiles
$\tilde{\psi}_{c/\sqrt{n}}$ and $\tilde{\psi}_{1-c/\sqrt{n}}$. Given that
estimated extreme quantiles are more variable than the estimated inner
quantiles, we set a given whisker equal to its associated extreme quantile
only when that quantile is at least $\lambda_{n}\times$100% shorter than the
quantity based on the noisy interquartile range. We take
$\lambda_{n}=n^{-1/4}$ and c=$1/20$. Several others were tested in simulation,
see Appendix B.2. Note that we use the unbounded algorithm because JointExp
(and, by consequence, PrivateQuantile and ApproxQuantile) are inconsistent for
estimating quantiles that satisfy $q_{n}=O(1/n)$, see Lemma 9.

A traditional boxplot also plots the observations which lie above and below
the upper and lower whiskers, respectively. We cannot release such data points
under the constraints of differential privacy. Instead, we plot a noisy
version of the number of points above $\tilde{u}$, denoted $\tilde{o}_{u}$ and
below $\tilde{\ell}$, denoted $\tilde{o}_{l}$. These are generated via the
Laplace mechanism, where it is easy to see that the count of observations
above or below a threshold has global sensitivity 1. Lastly, we attribute less
privacy budget for computing the outliers, as we deem these values to be of
less interest than the box itself. We assign each noisy outlyingness number a
privacy budget of $\epsilon/16$. Together, these seven values make up the
differentially private boxplot. See Figure 3 for an example of the proposed
private boxplot. The algorithm, which we call DPBoxplot is summarized in
Algorithm 1.

######  Remark 1 (On the tuning parameters).

Our algorithm contains several tuning parameters, which were chosen via
simulation. That is, we chose a $O(n^{-1/2})$ quantile, rather than a
$O(n^{-1})$ quantile, the constant $c=1/20$ and $\lambda_{n}=n^{-1/4}$, which
showed the best performance in our simulation study. See Appendix B.2 for
simulation results concerning the tuning parameters.

input : $\ \mathbf{X}_{n},\epsilon,a,b,\lambda_{n}$

// Construct the private quantiles

$\tilde{\psi}_{1-c/\sqrt{n}}\leftarrow$unbounded($q,a,b,3\epsilon/16$);

$\tilde{\psi}_{c/\sqrt{n}}\leftarrow$unbounded($q,a,b,3\epsilon/16$);

$\tilde{\xi}_{q}\leftarrow$JointExp($q,a,b,\epsilon/2$);

$\tilde{\xi}_{1/4}\leftarrow\tilde{\xi}_{1/4}\wedge\tilde{\xi}_{1/2}$;

$\tilde{\xi}_{3/4}\leftarrow\tilde{\xi}_{3/4}\vee\tilde{\xi}_{1/2}$;

// Construct the whiskers and outlyingness counts

$\tilde{\ell}\leftarrow\tilde{\xi}_{1/4}-1.5(\tilde{\xi}_{3/4}-\tilde{\xi}_{1/4})$;

$\tilde{u}\leftarrow\tilde{\xi}_{3/4}+1.5(\tilde{\xi}_{3/4}-\tilde{\xi}_{1/4})$;

if  _ $\tilde{\psi}_{c/\sqrt{n}}>\lambda_{n}|\tilde{\ell}|+\tilde{\ell}$_ then

$\tilde{\ell}\leftarrow\tilde{\psi}_{c/\sqrt{n}}$;

$\tilde{o}_{\ell}\leftarrow 0$;

else

$\tilde{o}_{\ell}\leftarrow\sum_{i=1}^{n}{\mathbbm{1}\left\\{x_{i}<\tilde{\ell}\right\\}}+\operatorname*{{\rm
Laplace}}(0,\frac{1}{\epsilon/16})$;

end if

if  _ $\tilde{\psi}_{1-c/\sqrt{n}}<\tilde{u}-\lambda_{n}|\tilde{u}|$_ then

$\tilde{u}\leftarrow\tilde{\psi}_{1-c/\sqrt{n}}$;

$\tilde{o}_{u}\leftarrow 0$;

else

$\tilde{o}_{u}\leftarrow\sum_{i=1}^{n}{\mathbbm{1}\left\\{x_{i}>\tilde{u}\right\\}}+\operatorname*{{\rm
Laplace}}(0,\frac{1}{\epsilon/16})$;

end if

return
$\tilde{B}(\mathbf{X}_{n},\epsilon)=(\tilde{o}_{\ell},\tilde{\ell},\tilde{\xi}_{1/4},\tilde{\xi}_{1/2},\tilde{\xi}_{3/4},\tilde{u},\tilde{o}_{u})$

Algorithm 1 DPBoxplot

##  4 Theoretical results

In this section, we derive several results concerning the different elements
of the private boxplot. We now assume that the dataset $\mathbf{X}_{n}$
consists of $n$ independent, identically distributed random variables drawn
from a population measure $\nu\in\mathcal{M}_{1}(\mathbb{R})$. We first prove
an upper bound on the sample complexity of the inner quantiles generated from
JointExp. We then prove a minimax lower bound for privately estimating a
quantile from a population measure lying in a general set, which matches the
upper bound. We then prove that the whiskers and outlyingness numbers and are
weakly consistent for their population counterparts.

We now define a boxplot rigorously, as a function of a measure on the set of
real numbers, which is convenient mathematically. The non-private boxplot
constructed from the data is taken to be the boxplot computed on the empirical
measure of the dataset $\mathbf{X}_{n}$, denoted by $\hat{\nu}$. For a measure
$\nu\in\mathcal{M}_{1}(\mathbb{R})$, let
$\xi_{p,\nu}=\inf\\{x\in\mathbb{R}\colon p\leq F_{\nu}(x)\\}$ be the
associated $p$th quantile, $\min(\nu)=\inf\\{x\colon F_{\nu}(x)>0\\}$,
$\max(\nu)=\inf\\{x\colon F_{\nu}(x)=1\\}$ and
$\operatorname{IQR}(\nu)=\xi_{3/4,\nu}-\xi_{1/4,\nu}$. For a measure
$\nu\in\mathcal{M}_{1}(\mathbb{R})$, let
$\operatorname{IQR}(\nu)=\xi_{3/4,\nu}-\xi_{1/4,\nu}$. For $x,y\in\mathbb{R}$,
we let $x\wedge y$ ($x\vee y$) be the minimum (maximum) of $x$ and $y$. Now,
letting $\ell_{1,\nu}=\xi_{1/4,\nu}-\operatorname{IQR}(\nu)$ and
$u_{1,\nu}=\xi_{3/4,\nu}+\operatorname{IQR}(\nu)$, the whiskers are defined as
$\ell_{\nu}=\ell_{1,\nu}\vee\min(\nu)$ and $u_{\nu}=u_{1,\nu}\wedge\max(\nu),$
Lastly, we can define
$o_{\ell,\nu}=F_{\nu}(\ell_{\nu})-\nu(\\{X=\ell_{\nu}\\})$,
$o_{u,\nu}=1-F_{\nu}(u_{\nu})$. The ‘population’ boxplot is the following
seven number summary
$B(\nu)=(o_{\ell,\nu},\ell_{\nu},\xi_{1/4,\nu},\xi_{1/2,\nu},\xi_{3/4,\nu},u_{\nu},o_{u,\nu})$.
Similarly, we can define the differentially private boxplot
$\tilde{B}(\mathbf{X}_{n},\epsilon)=\tilde{B}(\hat{\nu},\epsilon)=(\tilde{o}_{\ell},\tilde{\ell},\tilde{\xi}_{1/4},\tilde{\xi}_{1/2},\tilde{\xi}_{3/4},\tilde{u},\tilde{o}_{u})$.

For $L,r>0$ and $q\in(0,1)$, let $\mathcal{G}_{L,r,q}$ be the set of
absolutely continuous measures $\nu\in\mathcal{M}_{1}(\mathbb{R})$ such that
$\inf_{\xi_{q}-r\leq x\leq\xi_{q}+r}f_{\nu}(x)\geq L>0$. Our first result
gives an upper bound on the sample complexity of the private quantiles
generated via JointExp. For $0<q_{1}<\ldots,<q_{m}<1$ and
$\nu\in\mathcal{M}_{1}(\mathbb{R})$, let $q=(q_{1},\ldots,q_{m})$ and define
$\xi_{q,\nu}=(\xi_{q_{1},\nu},\ldots,\xi_{q_{m},\nu})$.

######  Theorem 2.

For $-\infty<a<b<\infty$ and $q=(q_{1},\ldots,q_{m})$ so that
$0<q_{1}<\ldots<q_{m}<1$, if $\nu\in\cap_{i=1}^{m}G_{L,r,q_{i}}$ is such that
$\sup_{x\in[a,b]}f_{\nu}(x)\leq K<\infty$ and
$a\leq\xi_{q_{1},\nu}\leq\ldots\xi_{q_{m},\nu}\leq b$, then there exists a
universal constant $C>0$ such that for all $0<\gamma<1$ and $t\leq r$, it
holds that $\lVert\xi_{q,\tilde{\nu}}-\xi_{q,\nu}\rVert\leq t$ with
probability $1-\gamma$, provided that

$n\geq C\left(m^{5/2}\frac{\log(1/\gamma)\vee
m\log(m/ctL)}{t^{2}L^{2}}\bigvee\
m^{2}\frac{\log\left(1/\gamma\right)\vee\log\left(\frac{K(b-a)}{q_{1}\wedge(1-q_{m})\wedge
tL}\right)}{tL\epsilon}\right).$

The proof of Theorem 2 relies on a general bound for the exponential mechanism
detailed by Ramsay et al., (2024), and so we defer it to Appendix A. Note that
the upper bound given in Theorem 2 has suboptimal scaling in $m$, the
ApproxQuantile algorithm of Kaplan et al., (2022) obtains logarithmic scaling
in $m$. However, in this context, $m=3$, and we are not so concerned about
this fact.

Next, we present a minimax lower bound for estimating a single quantile
subject to differential privacy, which in turn implies a lower bound for
estimating $m$ quantiles. Again, since our main interest is in boxplot
estimation, we are not concerned about optimal scaling in $m$. Next, for a set
$S$, let $\tilde{\mathcal{M}}_{1,\epsilon,n}(S)$ be the set of maps from
$\hat{\nu}$ to $\mathcal{M}_{1}(S)$ which satisfy (1). Note
$\tilde{\mathcal{M}}_{1,\epsilon,n}(S)$ is just a formal way of writing the
set of all differentially private estimators which lie in $S$. Let
$h_{q}(x)=x\sqrt{(-\log x-[\Phi^{-1}\left(q\right)]^{2}/2)/\pi}$,
$C_{q}=\operatorname*{{\rm argmax}}_{x>0}h_{q}(x)$ and denote the minimax risk
for differentially private quantile estimation by
$\mathcal{R}(\epsilon,L,r,q)=\inf_{T_{\epsilon}\sim\tilde{\mathcal{M}}_{1,\epsilon,n}(\mathbb{R})}\sup_{\nu\in\mathcal{G}_{L,r,q}}\mathbb{E}|T_{\epsilon}(\hat{\nu})-\xi_{q,\nu}|$.

######  Theorem 3.

For all $n\geq 1$, $q\in(0,1)$, $L,r>0$ such that $rL\leq\sup_{x>0}h_{q}(x)$,
it holds that

$n=\Omega\left(\frac{C_{q}^{2}}{L^{2}\tau^{2}}\vee\frac{C_{q}}{L\tau\epsilon}\right),$
(2)

samples are required for $\mathcal{R}(\epsilon,L,r,q)\leq t$ to hold.

Note that Theorem 3 implies a lower bound on the sample complexity for
estimating $m$ quantiles. That is, for all $q=(q_{1},\ldots,q_{m})$, $n\geq
1$, $L,r>0$ with $rL\leq\inf_{j\in[m]}C_{q_{j}}$, it holds that (2) samples
are also required for
$\inf_{T_{\epsilon}\sim\tilde{\mathcal{M}}_{1,\epsilon,n}(\mathbb{R}^{m})}\sup_{\nu\in\mathcal{G}_{L,r,q}}\mathbb{E}\lVert
T_{\epsilon}(\hat{\nu})-\xi_{q,\nu}\rVert\leq t$. Applying this lower bound in
conjunction with Theorem 2 yields that the scale and location of the proposed
private boxplot are estimated optimally, up to logarithmic factors. Note that
Tzamos and Vlatakis-Gkaragkounis, (2020) give a minimax lower bound on
differentially private median estimation. Their result is similar, though the
proof does not directly extend to estimating an arbitrary quantile, since it
relies on a set of uniform distributions who have that the property that the
median coincides with the mean. The proof of Theorem 3 is an elementary
application of the differentially private version of Fano’s inequality
(Acharya et al.,, 2021), it is deferred to Appendix A.

######  Remark 4.

Our theoretical exposition is based on JointExp, but any algorithm whose
quantiles satisfy the same sample complexity bound given by Theorem 2 would
result in optimal estimates of location and scale of the boxplot. For
instance, PrivateQuantile of Smith, (2011) is a special case of JointExp, when
only one quantile is being estimated, and therefore, the same results apply to
PrivateQuantile. In addition, Lalanne et al., 2023a  give concentration
results for the quantiles generated by the algorithm of Kaplan et al., (2022),
would also obtain the minimax lower bound given in Theorem 3. Lastly, recall
that JointExp is inconsistent for estimating extreme quantiles (see Lemma 9).
The proof of Lemma 9 also implies that PrivateQuantile and ApproxQuantile are
inconsistent for estimating extreme quantiles.

Next, we show that the whiskers and the outlyingness numbers, when
appropriately normalized, are weakly consistent for their population
counterparts. We first require a condition on the population measure.

######  Condition 1.

There exists $K,r,L>0$ such that
$\nu\in\mathcal{G}_{L,r,1/4}\cap\mathcal{G}_{L,r,1/2}\cap\mathcal{G}_{L,r,3/4}$
and $\sup_{x\in[a,b]}f_{\nu}(x)\leq K<\infty$.

Condition 1 requires that the population distribution have a bounded density,
which is also bounded below in a neighborhood of the median and the quartiles.
For instance, it is easily seen that Condition 1 is satisfied by the Gaussian,
student $t$, Gamma, Beta, and uniform distribution families.

######  Theorem 5.

For $-\infty<a<b<\infty$, if $\lambda_{n}\to 0$, Condition 1 holds and
$a\leq\xi_{1/4,\nu}<\xi_{3/4,\nu}\leq b$, then, for all $\beta>1$,
$q_{n}\rightarrow 0$ such that $nq_{n}\rightarrow\infty$ it holds that
$\tilde{\ell}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}\ell_{\nu}$,
$\tilde{u}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}u_{\nu}$
$\tilde{o}_{\ell}/n\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}o_{\ell,\nu}$
and
$\tilde{o}_{u}/n\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}o_{u,\nu}$.

Theorem 5 implies that for large sample sizes, the skewness and tails will be
correctly portrayed by the differentially private boxplot. Together, Theorems
2 and 5 imply that, given a large enough sample size, DPBoxplot will correctly
represent the location, scale, skewness and tail of the underlying
distribution. Thus, making it statistically valid for one to use the private
boxplot to describe and compare multiple samples in a differentially private
exploratory data analysis.

##  5 Simulations

We now assess the empirical performance of the differentially private boxplot
across the four key metrics: scale, location, skewness and tails, with a
simulation study. We consider the error between the private boxplot and its
non-private population counterpart. The errors in each of the component are
quantified as follows: $|\tilde{\xi}_{1/2}-\xi_{1/2,\hat{\nu}}|$ (location),
$|\tilde{\operatorname{IQR}}-\operatorname{IQR}(\hat{\nu})|$ (scale),
$|\tilde{\ell}-\ell_{\hat{\nu}}|+|\tilde{u}-u_{\hat{\nu}}|$ (skewness),
$|\tilde{o}_{\ell}-no_{\ell,\hat{\nu}}|+|\tilde{o}_{u}-no_{u,\hat{\nu}}|$
(tails).

We tested our proposed algorithm against the “naive private boxplot,” where a
single differentially private quantile algorithm is used to generate all
quantiles that make up the differentially private boxplot. We consider the
following quantile estimation methods: PrivateQuantile (Smith,, 2011),
unbounded (Durfee,, 2023), ApproxQuantile (Kaplan et al.,, 2022) and JointExp
(Gillenwater et al.,, 2021). For all methods, we used the default values of
all parameters according to their respective authors. For the input bounds on
the data, we set $a=-50$ and $b=50$, whenever input bounds were required. In
addition to the nonprivate boxplot, we also compute the population, or oracle
boxplot. For comparison, we also present the error in each component between
the non-private boxplot and the oracle boxplot.

Data were generated from the five distinct distributions parametrized to have
mean 0 and variance 1: normal distribution (normal), skew normal distribution
(skew), uniform distribution (uniform), beta distribution (beta), and
distribution of data from 2019 NY airbnb listing prices (airbnb)222These data
were sampled without replacement.. The results from the beta and airbnb
distributions did not convey any additional insights, and so we defer these to
Appendix B.1 We considered multiple sample sizes $n$ and privacy budgets
$\epsilon$. Each scenario was simulated 1000 times. (For more details, see
Appendix B.1).

Figure 1: Average error in each metric with 95% confidence intervals between
the different differentially private boxplots and their population
counterparts (y-axis) with $\epsilon\in\\{1,5\\}$ (line style) and increasing
$n$ (x-axis). The DPBoxplot algorithm is presented, along with the naive
method paired with each of the algorithms JointExp, ApproxQuantile, unbounded
and PrivateQuantile (line color). The different generated distributions are
represented row-wise, and the different metrics: location, scale, skewness and
tails are presented column-wise. The oracle line is the distance between the
non-private boxplot and its population counterpart.

Figure 1 provides a comprehensive summary of the simulation results. Each
column of figures corresponds to a distinct boxplot component, while each row
pertains to a different generating distribution. The sample size is depicted
along the $x$-axis, while the average error is represented on the $y$-axis.
The line color within the figures denotes the method used to generate the
boxplot and the line style corresponds to the privacy budget $\epsilon$.

Consistent with Theorems 2, 3, and 5, the error for the DPBoxplot is
converging to the oracle in all aspects under all simulated distributions with
increasing sample sizes. We also observe a slower convergence in skewness and
tails compared to location and scale, which is expected as distributional
extremes are more difficult to estimate privately. This shows that the error
between the proposed differentially private boxplot and its population
counterpart is similar to that of the error between the non-private boxplot
and its population counterpart. That is, the sampling error is larger than
that of the error attributed to privatization.

Observe that the naive boxplots do not exhibit the same behavior. In
particular, the inconsistency of the whiskers for the naive methods based on
JointExp, ApproxQuantile, and PrivateQuantile leads to very poor performance
in the skew metric for the skew and uniform distributions. This is due to the
inconsistency of all three of these algorithms for estimating extreme
quantiles. On the other hand, though it does well at estimating the extreme
quantiles, the unbounded algorithm exhibits poor performance in terms of
estimating scale and location, relative to the other algorithms. We observe
that, in general, DPBoxplot method exhibits consistently better behavior,
except notable with small sample sizes in the normal distribution. This
happens because at small sample sizes, the unbounded algorithm tends to
underestimate the magnitude of the minimum and maximum of the dataset, while
the other algorithms overestimate this magnitude.

We also assessed the ability to compare multiple distributions through
simulation, concluding that DPBoxplot also performs well under this setting.
For brevity, this study is defered to Appendix B.3

##  6 Case study

We now conduct an exploratory data analysis via boxplots, within the framework
of differential privacy. To maintain the authenticity of the case study,
privacy budgets are proportionally distributed across all visualizations
explored in the process of addressing each question. To elaborate, out of the
total budget ($\epsilon=1$), each visualization is assigned a privacy budget
equal to the number of boxplots in the visualization, divided by the number of
boxplots on all generated visualizations. The purpose for this allocation is
to provide more budget for visualizations where the data is partitioned more
times and hence sample sizes are smaller. Code for replicating these findings
has been made publicly accessible .

We analyze a dataset containing Airbnb listing prices and associated metrics
within New York City (NYC) in 2019 (Kaggle,, 2019). After removing listings
priced above 500 US dollars (USD) and requiring minimum nights of stay fewer
than 10, this dataset has $n=40738$ observations and $d=4$ explanatory
variables of business interest. We only consider listings priced below 500
USD, and so we set $a=0$ and $b=500$. We address two distinct business
inquiries:

Figure 2: Private and non-private boxplots for inquiry 1. The boxplots display
listing prices against borough, room type and room type separated by borough
in columns 1, 2 and 3, respectively. Private boxplots are obtained with a
total aggregated budget of $\epsilon=1$. Figure 3: Private and non-private
boxplots for inquiry 2. The boxplots display listing prices against minimum
nights and room type separated by minimum nights in columns 1 and 2,
respectively. Private boxplots are obtained with a total aggregated budget of
$\epsilon=1$.

#### Inquiry 1: Do discernible patterns emerge in Airbnb listing prices across
various boroughs in New York City and differing room types?

The dataset encompasses five distinct boroughs within New York City, namely
Bronx, Brooklyn, Manhattan, Queens, and Staten Island, alongside three offered
room types: Entire Home, Private Room, and Shared Room. We present three
visualizations: boxplots of prices by borough, by room type, and a combined
visualization of prices by borough and room type. These generate 5, 3, and 15
boxplots, respectively, totaling 23.

The differentially private boxplots are displayed in Figure 2, juxtaposed with
their non-private counterparts. Only looking at the differentially private
boxplots, Figure 2 reveals that prices predominantly lie in the bottom end of
the range [0,500], with a right skew and heavy right tail observed across all
boroughs and room types, except for in the Bronx, and for shared rooms. The
distribution of prices in the Bronx still exhibits a right skew, but has light
tails. The distribution of prices for shared rooms appears to be symmetric,
with light tails. Notably, prices appear elevated for Manhattan. As expected,
entire homes are priced higher than shared spaces. The right-most plot, which
displays prices by both borough and room type, affirms previous observations,
except shared rooms in Staten Island and the Bronx. Staten Island and the
Bronx exhibits higher variability for shared rooms.

It is natural to compare the patterns observed on the differentially private
boxplots to those observed in the non-private boxplots. As previously
mentioned, Figure 2 juxtaposes the differentially private boxplots with the
non-private boxplots. Many conclusions derived from the private visualizations
persist in the non-private domain. However, there are several disparities,
which are outlined are as follows: In reality, both the listing prices for
shared rooms and properties in the Bronx do have a heavy, right tail. The
second disparity occurs in the third plot, where the non-private plots does
not indicate higher prices and variability for properties in Staten Island and
the Bronx in shared rooms. These discrepancies can be explained by small
sample sizes and our conservative choice of privacy budget. For instance, the
number of shared room listings in Staten Island is only 9. Low sample sizes
underscore potential for greater discrepancies between private and non-private
boxplot representations. It may be of interest to allocate more privacy budget
to visualizations generated from smaller subsets of the dataset.

#### Inquiry 2: Are there observable trends in Airbnb listing prices
concerning minimum nights required for reservation and the types of rooms
offered?

We create a new categorical variable called minimum nights, in which a listing
is assigned “low” if it has less than or equal to three minimum nights, and
“high” otherwise. For this inquiry, we generate two visualizations: Boxplots
of prices by minimum nights, and prices by combinations of minimum nights and
room offered. These visualizations require 2, and 6 boxplots, respectively,
totaling 8.

The differentially private boxplots are displayed in Figure 3, juxtaposed with
their non-private counterparts. Again, we first analyze the patterns in the
differentially private boxplots, and compare to those observed in the private
boxplots afterward. The differentially private boxplots indicate the emergence
of a phenomenon akin to Simpson’s paradox. Specifically, a preliminary
examination suggests that listings requiring a higher minimum number of nights
are priced more steeply than their counterparts with lower minimum
requirements. However, this trend disappears when the data is divided by room
type. Listings for entire rooms exhibit no significant price differential
based on minimum night requirement, while private rooms slightly favor lower
minimum nights in terms of price. In shared rooms the median price does not
seem to be significantly different.

A comparative analysis with non-private boxplots reaffirms these observations;
we observe relatively consistent patterns between the private and non-private
boxplots. The primary visual disparity pertains to the positioning of the
lower whiskers on the boxplots. This underscores the recognized challenge
associated with differentially private estimation of extreme quantiles.
However, this discrepancy does not materially impede the analytical value of
our findings.

Overall, the patterns observed are generally consistent between the private
and non-private boxplots. The principal disparities are attributable to sample
size, a factor that we encourage practitioners to consider when conducting
differentially private exploratory analyses. This case study demonstrates the
potential of differentially private, exploratory data analysis and confirms
the efficacy of the differentially private boxplot, in accordance with our
theoretical and simulated results.

## References

  * Acharya et al., (2021) Acharya, J., Sun, Z., and Zhang, H. (2021).  Differentially private Assouad, Fano, and Le Cam.  In Feldman, V., Ligett, K., and Sabato, S., editors, Proceedings of the 32nd International Conference on Algorithmic Learning Theory, volume 132 of Proceedings of Machine Learning Research, pages 48–78. PMLR.
  * Acs et al., (2012) Acs, G., Castelluccia, C., and Chen, R. (2012).  Differentially private histogram publishing through lossy compression.  In 2012 IEEE 12th International Conference on Data Mining, pages 1–10. IEEE.
  * Alabi et al., (2022) Alabi, D., Ben-Eliezer, O., and Chaturvedi, A. (2022).  Bounded space differentially private quantiles.
  * Bhattacharjee et al., (2020) Bhattacharjee, K., Chen, M., and Dasgupta, A. (2020).  Privacy-preserving data visualization: Reflections on the state of the art and research opportunities.  Computer Graphics Forum, 39:675–692.
  * Budiu et al., (2022) Budiu, M., Thaker, P., Gopalan, P., Wieder, U., and Zaharia, M. (2022).  Overlook: Differentially private exploratory visualization for big data.  Journal of Privacy and Confidentiality, 12.
  * Dasgupta et al., (2013) Dasgupta, A., Chen, M., and Kosara, R. (2013).  Measuring privacy and utility in privacy-preserving visualization.  Computer Graphics Forum, 32:35–47.
  * Dasgupta and Kosara, (2011) Dasgupta, A. and Kosara, R. (2011).  Adaptive privacy-preserving visualization using parallel coordinates.  IEEE Transactions on Visualization and Computer Graphics, 17(12):2241–2248.
  * Dasgupta et al., (2019) Dasgupta, A., Kosara, R., and Chen, M. (2019).  Guess me if you can: A visual uncertainty model for transparent evaluation of disclosure risks in privacy-preserving data visualization.  In 2019 IEEE Symposium on Visualization for Cyber Security (VizSec), pages 1–10.
  * Durfee, (2023) Durfee, D. (2023).  Unbounded differentially private quantile and maximum estimation.  In Oh, A., Neumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S., editors, Advances in Neural Information Processing Systems, volume 36, pages 77691–77712. Curran Associates, Inc.
  * Dwork, (2008) Dwork, C. (2008).  Differential privacy: A survey of results.  In International conference on theory and applications of models of computation, pages 1–19. Springer.
  * Dwork et al., (2006) Dwork, C., McSherry, F., Nissim, K., and Smith, A. (2006).  Calibrating noise to sensitivity in private data analysis.  In Theory of Cryptography Conference, pages 265–284.
  * Garrido et al., (2023) Garrido, G. M., Liu, X., Matthes, F., and Song, D. (2023).  Lessons learned: Surveying the practicality of differential privacy in the industry.  Proceedings on Privacy Enhancing Technologies, 2023:151–170.
  * Gillenwater et al., (2021) Gillenwater, J., Joseph, M., and Kulesza, A. (2021).  Differentially private quantiles.  In Meila, M. and Zhang, T., editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 3713–3722. PMLR.
  * Hay et al., (2010) Hay, M., Rastogi, V., Miklau, G., and Suciu, D. (2010).  Boosting the accuracy of differentially private histograms through consistency.  Proc. VLDB Endow., 3(1–2):1021–1032.
  * He et al., (2016) He, X., Raval, N., and Machanavajjhala, A. (2016).  A demonstration of visdpt.  Proceedings of the VLDB Endowment, 9:1489–1492.
  * Hongde et al., (2014) Hongde, R., Shuo, W., and Hui, L. (2014).  Differential privacy data aggregation optimizing method and application to data visualization.  In 2014 IEEE Workshop on Electronics, Computer and Applications, pages 54–58.
  * Ji et al., (2014) Ji, Z., Lipton, Z. C., and Elkan, C. (2014).  Differential privacy and machine learning: a survey and review.
  * Kaggle, (2019) Kaggle (2019).  New york city airbnb open data.
  * Kaplan et al., (2022) Kaplan, H., Schnapp, S., and Stemmer, U. (2022).  Differentially private approximate quantiles.  In Chaudhuri, K., Jegelka, S., Song, L., Szepesvari, C., Niu, G., and Sabato, S., editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 10751–10761. PMLR.
  * Kellaris and Papadopoulos, (2013) Kellaris, G. and Papadopoulos, S. (2013).  Practical differential privacy via grouping and smoothing.  Proceedings of the VLDB Endowment, 6(5):301–312.
  * (21) Lalanne, C., Garivier, A., and Gribonval, R. (2023a).  Private statistical estimation of many quantiles.  In Krause, A., Brunskill, E., Cho, K., Engelhardt, B., Sabato, S., and Scarlett, J., editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 18399–18418. PMLR.
  * (22) Lalanne, C. S., Gastaud, C., Grislain, N., Garivier, A., and Gribonval, R. (2023b).  Private quantiles estimation in the presence of atoms.
  * Li et al., (2010) Li, C., Hay, M., Rastogi, V., Miklau, G., and McGregor, A. (2010).  Optimizing linear counting queries under differential privacy.  In Proceedings of the twenty-ninth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database systems, pages 123–134.
  * Liu et al., (2023) Liu, W., Zhang, Y., Yang, H., and Meng, Q. (2023).  A survey on differential privacy for medical data analysis.  Annals of Data Science, 11(2):733–747.
  * Nanayakkara et al., (2022) Nanayakkara, P., Bater, J., He, X., Hullman, J., and Rogers, J. (2022).  Visualizing privacy-utility trade-offs in differentially private data releases.  Proceedings on Privacy Enhancing Technologies, 2022:601–618.
  * Panavas et al., (2024) Panavas, L., Crnovrsanin, T., Adams, J. L., Ullman, J., Sargavad, A., Tory, M., and Dunne, C. (2024).  Investigating the visual utility of differentially private scatterplots.  IEEE Transactions on Visualization and Computer Graphics, pages 1–16.
  * Ramsay et al., (2024) Ramsay, K., Jagannath, A., and Chenouri, S. (2024).  Differentially private multivariate medians.
  * Smith, (2011) Smith, A. (2011).  Privacy-preserving statistical estimation with optimal convergence rates.  In Proceedings of the forty-third annual ACM symposium on Theory of computing, STOC’11. ACM.
  * Tukey et al., (1977) Tukey, J. W. et al. (1977).  Exploratory data analysis, volume 2.  Springer.
  * Tzamos and Vlatakis-Gkaragkounis, (2020) Tzamos, C. and Vlatakis-Gkaragkounis (2020).  Optimal private median estimation under minimal distributional assumptions.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 3301–3311.
  * Xu et al., (2013) Xu, J., Zhang, Z., Xiao, X., Yang, Y., Yu, G., and Winslett, M. (2013).  Differentially private histogram publication.  The VLDB journal, 22:797–822.
  * Zhang et al., (2016) Zhang, D., Hay, M., Miklau, G., and O’connor, B. (2016).  Challenges of visualizing differentially private data.
  * Zhang et al., (2014) Zhang, X., Chen, R., Xu, J., Meng, X., and Xie, Y. (2014).  Towards accurate histogram publication under differential privacy.  pages 587–595. Society for Industrial and Applied Mathematics.
  * Zhou et al., (2022) Zhou, J., Wang, X., Wong, J. K., Wang, H., Wang, Z., Yang, X., Yan, X., Feng, H., Qu, H., Ying, H., and Chen, W. (2022).  Dpviscreator: Incorporating pattern constraints to privacy-preserving visualizations via differential privacy.  IEEE Transactions on Visualization and Computer Graphics, pages 1–11.

##  Appendix A Technical proofs

Before presenting the proof of Theorem 2, we first prove a general sample
complexity bound for quantiles. For $t\geq 0$, $q=(q_{1},\ldots,q_{m})$ such
that $0<q_{1}<\ldots<q_{m}<1$, $-\infty<a<b<\infty$, and
$\nu\in\mathcal{M}_{1}(\mathbb{R})$, define:

$\displaystyle\alpha(t)=\inf_{a<x_{1}<\ldots<x_{m}<b,\ \left\lVert
x-\xi_{q,\nu}\right\rVert\geq
t}\sum_{i=1}^{m+1}|q_{j}-q_{j-1}-F_{\nu}(x_{j})+F_{\nu}(x_{j-1})|.$

One may wish to recall that $\inf\emptyset=\infty$, and so $\alpha(t)$ is
still defined if $\\{a<x_{1}<\ldots<x_{m}<b,\ \left\lVert
x-\xi_{q,\nu}\right\rVert\geq t\\}=\emptyset$.

######  Lemma 6.

For $-\infty<a<b<\infty$, if $\nu\in\mathcal{M}_{1}(\mathbb{R})$ is absolutely
continuous such that $\sup_{x\in[a,b]}f_{\nu}(x)\leq K<\infty$, then there
exists universal constants $C,c>0$ such that for all $q=(q_{1},\ldots,q_{m})$
such that $0<q_{1}<\ldots<q_{m}<1$ and
$a\leq\xi_{q_{1},\nu}\leq\ldots\xi_{q_{m},\nu}\leq b$, all $t\geq 0$ and
$0<\gamma<1$, it holds that $\lVert\xi_{q,\tilde{\nu}}-\xi_{q,\nu}\rVert\leq
t$, with probability $1-\gamma$, provided that

$n\geq Cm^{2}\frac{\log(1/\gamma)\vee m\log
m/c\alpha(t)}{\alpha^{2}(t)}\bigvee\
m\frac{\log\left(1/\gamma\right)\vee\log\left(\frac{K(b-a)}{q_{1}\wedge(1-q_{m})|\wedge\alpha(t)}\right)}{\alpha(t)\epsilon}.$
(3)

###### Proof of Lemma 6.

Given that $\tilde{\xi}_{q,\hat{\nu}}$ is a draw from the exponential
mechanism, the result follows from an application of Corollary 7 of Ramsay et
al., (2024). Corollary 7 of Ramsay et al., (2024) gives an upper bound on the
sample complexity of a draw from the exponential mechanism, provided the
utility function meets certain criteria. In order to apply Corollary 7 of
Ramsay et al., (2024), we must show that the following function

$\phi^{\prime}(x,\nu)=-\sum_{j=1}^{m+1}|F_{\nu}(x_{j})-F_{\nu}(x_{j-1})-(q_{j}-q_{j-1})|\mathbbm{1}\left\\{[x_{1},x_{m}]\in[a,b]\right\\}\\\
-\infty\cdot(1-\mathbbm{1}\left\\{[x_{1},x_{m}]\in[a,b]\right\\}),$

and $\nu$ satisfy three conditions. First, we must show that
$\phi^{\prime}(x,\nu)$ has a maximum, which is easily seen by definition,
$\phi(x,\nu)$ is maximized at $\xi_{q,\nu}$. The second requirement is that
$\phi^{\prime}(x,\nu)$ is $K$-Lipschitz, which follows by assumption. Lastly,
Corollary 7 of Ramsay et al., (2024) requires that $\phi^{\prime}(x,\nu)$ is a
$(L,\mathcal{F})$ regular function for some $L>0$ with
$\operatorname{VC}(\mathcal{F})<\infty$, see Definition 2. It is easy to see
that $\phi^{\prime}(x,\nu)$ is $(m,\mathcal{F})$ regular function with
$\operatorname{VC}(\mathcal{F})=1$. In order to apply the bound, we must
compute the discrepancy function of $\phi^{\prime},\nu$, which is given by:
$\alpha^{\prime}(t)=\phi^{\prime}(\xi_{q,\nu},\nu)-\sup_{\left\lVert
x-\xi_{q,\nu}\right\rVert\geq t}\phi^{\prime}(x,\nu).$ A straightforward
calculation yields that $\alpha^{\prime}(t)=\alpha(t)$. Applying Corollary 7
of Ramsay et al., (2024), yields that there exists a universal constant $C>0$
such that $\lVert\xi_{q,\tilde{\nu}}-\xi_{q,\nu}\rVert\leq t$ with probability
at least $1-\gamma$, if

$n\geq C\left(m^{2}\frac{\log(1/\gamma)\vee m\log
m/c\alpha(t)}{\alpha^{2}(t)}\bigvee\
m\frac{\log\left(1/\gamma\right)\vee\log\left(\frac{(b-a)}{|\xi_{q_{1},\nu}-a|\wedge|b-\xi_{q_{m},\nu}|\wedge\alpha(t)/K}\right)}{\alpha(t)\epsilon}\right).$

Define $a\gtrsim b\ (a\lesssim b)$ if there exists a universal constant $c>0$
such that $a\geq cb$ ($a\leq cb$). Next, using the Lipschitz assumption and
the mean value theorem, we have that
$|\xi_{q_{1},\nu}-a|\wedge|b-\xi_{q_{m},\nu}|\gtrsim[q_{1}\wedge(1-q_{m})]/K,$
which yields the desired result. ∎

We can now prove Theorem 2.

###### Proof of Theorem 2.

The proof is based on an application of Lemma 6. All of the conditions of
Lemma 6 are satisfied by assumption, and so it remains to lower bound
$\alpha$. To this end, let $k_{x}=\operatorname*{{\rm
argmax}}_{i\in[m]}|x_{i}-\xi_{q_{i},\nu}|$. Using this notation, for $0\leq
t\leq r$, the definition of $\alpha$ in conjunction with the mean value
theorem yields that

$\displaystyle\alpha(t)$ $\displaystyle=\inf_{a<x_{1}<\ldots<x_{m}<b,\
\left\lVert x-\xi_{q,\nu}\right\rVert\geq
t}\sum_{j=1}^{m+1}|q_{j}-q_{j-1}-F_{\nu}(x_{j})+F_{\nu}(x_{j-1})|$
$\displaystyle\geq\inf_{a<x_{1}<\ldots<x_{m}<b,\ \left\lVert
x-\xi_{q,\nu}\right\rVert\geq
t}\sum_{j=1}^{k_{x}}|q_{j}-q_{j-1}-F_{\nu}(x_{j})+F_{\nu}(x_{j-1})|$
$\displaystyle\geq\inf_{a<x_{1}<\ldots<x_{m}<b,\ \left\lVert
x-\xi_{q,\nu}\right\rVert\geq
t}|\sum_{j=1}^{k_{x}}q_{j}-q_{j-1}-F_{\nu}(x_{j})+F_{\nu}(x_{j-1})|$
$\displaystyle\geq\inf_{a<x_{1}<\ldots<x_{m}<b,\ \left\lVert
x-\xi_{q,\nu}\right\rVert\geq t}|q_{k_{x}}-F_{\nu}(x_{k_{x}})|$
$\displaystyle\geq tL/\sqrt{m}.$

Applying this bound, in conjunction with Lemma 6, yields the desired result. ∎

Next, we prove Theorem 3.

###### Proof of Theorem 3.

We apply Corollary 4 of (Acharya et al.,, 2021), of which a simpler version is
restated below for clarity. For
$\mathcal{P}\subset\mathcal{M}_{1}(\mathbb{R})$, let
$n^{*}=\inf\\{n\colon\inf_{T_{\epsilon}\sim\tilde{\mathcal{M}}_{1,\epsilon,n}(\mathbb{R})}\sup_{\nu\in\mathcal{P}}\int|T_{\epsilon}(\hat{\nu})-T_{\epsilon}(\nu)|d\nu^{n}\leq\tau\\}.$

######  Corollary 7 (Acharya et al., (2021)).

For all $\epsilon,\tau>0$ and any
$\mathcal{P}\subset\mathcal{M}_{1}(\mathbb{R})$, let
$\mathcal{Q}=\\{Q_{1},\ldots,Q_{m}\\}\subset\mathcal{P}$. If for all $i\neq
j$, it holds that $|T_{\epsilon}(Q_{i})-T_{\epsilon}(Q_{j})|\geq 3\tau$,
$\operatorname{KL}(Q_{i},Q_{j})\leq\beta$, and
$\operatorname{TV}(Q_{i},Q_{j})\leq\gamma$, then then $n^{*}=\Omega\left(\log
m(\beta^{-1}\vee(\gamma\epsilon)^{-1})\right)$.

In order to apply Corollary 7, we first define a class of Gaussians which lie
in $\mathcal{G}_{L,r,q}$. Recall that we consider $\nu$ with densities which
satisfy: $\inf_{x\in\xi_{q,\nu}\pm r}f_{\nu}(x)\geq L$. Let
$P_{\mu,\sigma}=\mathcal{N}(\mu,\sigma^{2})$ and let $h_{\mu,\sigma}$ denote
the associated quantile function. We have that for $q\leq 1/2$, it holds that

$\displaystyle\inf_{x\in h_{\mu,\sigma}(q)\pm r}f_{P_{\mu,\sigma}}(x)$
$\displaystyle\geq(2\pi\sigma^{2})^{-1/2}\exp\left(-[\text{erf}^{-1}(2q-1)-r/\sqrt{2}\sigma]^{2}\right)$
$\displaystyle\geq(2\pi\sigma^{2})^{-1/2}\exp\left(-r^{2}/2\sigma^{2}-[\text{erf}^{-1}(2q-1)]^{2}\right).$

Taking $\sigma=C_{q}/\sqrt{2\pi}L$ gives that

$\displaystyle\inf_{x\in h_{\mu,\sigma}(q)\pm r}f_{P_{\mu,\sigma}}(x)$
$\displaystyle\geq L/C_{q}\exp\left(-[\text{erf}^{-1}(2q-1)]^{2}-\pi
r^{2}L^{2}/C_{q}^{2}\right).$

Now, note that

$\displaystyle\text{erf}^{-1}(2q-1)$
$\displaystyle=\Phi^{-1}\left(q\right)/\sqrt{2},$

which in turn implies that

$\inf_{x\in h_{\mu,\sigma}(q)\pm r}f(x)\geq
L\exp\left(-[\Phi^{-1}\left(q\right)]^{2}/2-r^{2}L^{2}\right)/\sqrt{\pi}.$

Next, we have that

$\displaystyle L/C_{q}\exp\left(-[\Phi^{-1}\left(q\right)]^{2}/2-\pi
r^{2}L^{2}/C_{q}^{2}\right)\geq L$
$\displaystyle\iff\exp\left(-[\Phi^{-1}\left(q\right)]^{2}/2-\pi
r^{2}L^{2}/C_{q}^{2}\right)\geq C_{q}$
$\displaystyle\iff[\Phi^{-1}\left(q\right)]^{2}/2+\pi
r^{2}L^{2}/C_{q}^{2}\leq-\log C_{q}$ $\displaystyle\iff rL\leq
C_{q}\sqrt{(-\log C_{q}-[\Phi^{-1}\left(q\right)]^{2}/2)/\pi}.$

We must then have that $rL\leq C_{q}\sqrt{(-\log
C_{q}-[\Phi^{-1}(q)]^{2}/2)/\pi}$, which holds by assumption. Define
$\mathcal{F}_{L,r}=\left\\{\mathcal{N}(\mu,C_{q}^{2}/2\pi
L^{2})\colon\mu\in\mathbb{R}\right\\}$. We find the values of $\beta,\gamma$
and $\tau$ for $\mathcal{F}_{L,r}$. Consider the values $\mu_{0},\mu_{1}$, and
for $i$=1 or $i=0$, denote $Q_{i}=P_{\mu_{i},C_{q}^{2}/2\pi L^{2}}$. Note that
for any quantile $q$, we have that
$|\xi_{q,Q_{0}}-\xi_{q,Q_{1}}|=|\mu_{0}-\mu_{1}|$. That is, the distance
between the quantiles is just the difference between the means. For some
$\tau>0$, take any $\mu_{0},\mu_{1}$ such that $|\mu_{0}-\mu_{1}|\geq 3\tau$.
It follows that
$\operatorname{KL}(Q_{0},Q_{1})\lesssim\tau^{2}L^{2}/C_{q}^{2}$ and
$\operatorname{TV}(Q_{0},Q_{1})=2\Phi(c\tau L/C_{q})-1\lesssim L\tau/C_{q}$.
Therefore, applying Corollary 7 gives that

$n^{*}=\Omega\left(\frac{C_{q}^{2}}{L^{2}\tau^{2}}\vee\frac{C_{q}}{L\tau\epsilon}\right).\qed$

Before proceeding with the proof of Theorem 5, we first prove that the
unbounded algorithm estimates extreme quantiles consistently.

######  Lemma 8.

For all $\beta>1$, $q_{n}\rightarrow 0$ such that $nq_{n}\rightarrow\infty$,
it holds that

  1. i

if $\min(\nu)<b$, then
$\tilde{\psi}_{q_{n}}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}\min(\nu)$.

  2. ii

if $\max(\nu)>a$, then
$\tilde{\psi}_{1-q_{n}}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}\max(\nu)$.

###### Proof.

First, using the fact that $V_{0}$ is an exponential random variable, we have
that $\Pr\left(V_{0}>q_{n}\right)=e^{-nq_{n}\epsilon}\to 0$ as $n\to\infty$.
In addition, letting $c_{n}=1/\log n$, the Dvoretzky–Kiefer–Wolfowitz
inequality yields that
$\Pr\left(\sup_{x\in\mathbb{R}}|F_{\hat{\nu}}(x)-F_{\nu}(x)|>c_{n}\right)\leq
2e^{-2n/(\log n)^{2}}$. It suffices to show that
$\Pr\left(|\tilde{\xi}_{1-q_{n}}-\max(\nu)|>t\right)\rightarrow 0$ as
$n\to\infty$. (A symmetric argument then implies then that also,
$\Pr\left(|\tilde{\xi}_{q_{n}}-\min(\nu)|>t\right)\rightarrow 0$ as
$n\to\infty$.) We have that

$\displaystyle\Pr\left(|\tilde{\xi}_{1-q_{n}}-\max(\nu)|>t\right)$
$\displaystyle=\Pr\left(\max(\nu)-\beta^{k}-1\geq
t\right)+\Pr\left(\beta^{k}-1-\max(\nu)\geq t\right),$

for which we can write

$\displaystyle\Pr\left(\max(\nu)-\beta^{k}-1\geq t\right)$
$\displaystyle=\Pr\left(k\leq\frac{\log\left(\max(\nu)-t-1\right)}{\log\beta}\right)\coloneqq\Pr\left(k\leq
R_{t,1}\right)$ $\displaystyle\Pr\left(\beta^{k}-1-\max(\nu)\geq t\right)$
$\displaystyle=\Pr\left(k\geq\frac{\log\left(\max(\nu)+t+1\right)}{\log\beta}\right)\coloneqq\Pr\left(k\geq
R_{t,2}\right).$

It suffices to show that $\Pr\left(k\geq R_{t,1}\right),\Pr\left(k\geq
R_{t,2}\right)\to 0$ as $n\to\infty$. To this end, we have that

$\displaystyle\Pr\left(k\leq R_{t,1}\right)$
$\displaystyle=1-\prod_{i=1}^{R_{t,1}}\Pr\left(F_{\hat{\nu}}(\beta^{i}-1)+V_{i}<1-q_{n}+V_{0}\right)$
$\displaystyle\leq
1-\prod_{i=1}^{R_{t,1}}\Pr\left(F_{\nu}(\beta^{i}-1)+V_{i}<1-q_{n}-c_{n}\right)+2e^{-2n/(\log
n)^{2}}$ $\displaystyle\leq
1-\prod_{i=1}^{R_{t,1}}\left(1-\exp\left(-n\epsilon(1-q_{n}-c_{n}-F_{\nu}(\beta^{i}-1)\right)\right)+2e^{-2n/(\log
n)^{2}}$ $\displaystyle\coloneqq 1-I+2e^{-2n/(\log n)^{2}}.$

Now, for large $n$, we have that

$\displaystyle I$
$\displaystyle=\prod_{i=1}^{R_{t,1}}\left(1-\exp\left(-n\epsilon(1-q_{n}-c_{n}-F_{\nu}(\beta^{i}-1)\right)\right)$
$\displaystyle\geq\left(1-\exp\left(-n\epsilon(1-q_{n}-c_{n}-F_{\nu}(\beta^{R_{t,1}}-1)\right)\right)^{R_{t,1}}$
$\displaystyle\geq
1-R_{t,1}\exp\left(-n\epsilon(1-q_{n}-c_{n}-F_{\nu}(\beta^{R_{t,1}}-1)\right).$

where the last inequality uses the fact that $(1+x/n)^{n}\geq 1+x$ for all
$n\geq 1$ and $|x|\leq n$. Applying this inequality yields that

$\displaystyle\Pr\left(k\leq R_{t,1}\right)$ $\displaystyle\leq
R_{t,1}\exp\left(-n\epsilon(1-q_{n}-c_{n}-F_{\nu}(\beta^{R_{t,1}}-1)\right)+2e^{-2n/(\log
n)^{2}}\to 0,$

as $n\to\infty$. On the other hand, for the term $\Pr\left(k\geq
R_{t,2}\right)$, we have that

$\displaystyle\Pr\left(k\geq R_{t,2}\right)$
$\displaystyle\leq\prod_{i=1}^{R_{t,2}}\Pr\left(F_{\nu}(\beta^{i}-1)+V_{i}<1-q_{n}/2+c_{n}\right)+2e^{-2n/(\log
n)^{2}}+e^{-nq_{n}}$
$\displaystyle\leq\Pr\left(F_{\nu}(\beta^{R_{t,2}}-1)+V_{R_{t,2}}<1-q_{n}/2+c_{n}\right)+2e^{-2n/(\log
n)^{2}}+e^{-nq_{n}}$
$\displaystyle\leq\mathbbm{1}\left\\{1-q_{n}/2+c_{n}-F_{\nu}(\beta^{R_{t,2}}-1)\geq
0\right\\}+2e^{-2n/(\log n)^{2}}+e^{-nq_{n}}\to 0,$

as $n\to\infty$. ∎

We can now prove Theorem 5.

###### Proof of Theorem 5.

We first prove that the whiskers are consistent. Theorem 2 gives that
$\tilde{\ell}_{1}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}\ell_{1,\nu}$
and Lemma 8 gives that
$\tilde{\ell}_{2}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}\min(\nu)$.
Continuous mapping theorem and the fact that $\lambda_{n}\to 0$ as
$n\to\infty$ yields that
$\tilde{\ell}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}\ell_{\nu}$. The
same argument applies to the upper whisker.

For the outlyingness number, first, we have that
$|o_{\ell,\nu}/n-\tilde{o}_{\ell}|\leq|o_{\ell,\hat{\nu}}-\tilde{o}_{\ell}/n|+|o_{\ell,\nu}-o_{\ell,\hat{\nu}}|\coloneqq
I+II$. Next, the properties of the Laplace distribution give that
$\Pr\left(|o_{\ell,\hat{\nu}}-\tilde{o}_{\ell}/n|\geq t\right)\lesssim
e^{-n\epsilon t}$. Therefore,
$o_{\ell,\hat{\nu}}\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}\tilde{o}_{\ell}/n$,
or, $I\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}0$. Next, using the fact
that $\sup_{x\in\mathbb{R}}f_{\nu}(x)\leq K$, we have that $F_{\nu}$ is
$K$-Lipschitz. It follows that $|o_{\ell,\nu}/n-o_{\ell,\hat{\nu}}/n|\leq
K|\ell_{\hat{\nu}}-\ell_{\nu}|+\sup_{x\in\mathbb{R}}|F_{\nu}(x)-F_{\hat{\nu}}(x)|$.
Next, the Dvoretzky–Kiefer–Wolfowitz inequality yields that
$\sup_{x\in\mathbb{R}}|F_{\nu}(x)-F_{\hat{\nu}}(x)|\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}0$
and the fact that the whiskers are consistent implies that
$K|\ell_{\hat{\nu}}-\ell_{\nu}|\stackrel{{\scriptstyle\text{p}}}{{\rightarrow}}0$.
The same argument can be made for the upper outlyingness number. ∎

##  Appendix B Simulation details and additional results

###  B.1 Single boxplot estimation

This section presents more details on the simulation study described in 5. We
simulated data from five distributions, each with mean 0 and variance 1:
standard normal distribution (normal), skew normal distribution with scale
parameter 20 (skew), uniform distribution with interval $[-\sqrt{3},\sqrt{3}]$
(uniform), a normalized and mean-centered version of the beta distribution
with $\alpha=\beta=2$ (beta), and normalized and mean-centered empirical
distribution of data from 2019 NY airbnb listing prices sampled without
replacement. We considered sample sizes of
$n\in\\{1000,3500,10000,35000,100000\\}$ and privacy budgets of
$\epsilon\in\\{0.5,1,5,10\\}$. Each scenario was simulated 1000 times leading
to simulated data vectors $\mathbf{X}_{n}$. For each generated dataset, we
computed both the non-private boxplots $B(\mathbf{X}_{n})=B(\hat{\nu})$. We
also calculated differentially private boxplots
$\tilde{B}(\mathbf{X}_{n},\epsilon)$ using each of the methods. We then
quantified boxplot distances against location, scale, skewness, and tails
metrics between the differentially private boxplots and their corresponding
population boxplots. We also calculated the oracle distance between the sample
non-private boxplot and the population boxplot. Figures 4 and 5 give a more
detailed overview of the results from this simulation study.

Figure 4: Average error in each metric with 95% confidence intervals between
population boxplots and estimated differentially private boxplots (y-axis)
with $\epsilon\in\\{0.5,1,5,10\\}$ (line colors) and increasing $n$ (x-axis);
using DPBoxplot for private boxplot estimation; from all sampled distributions
(column-wise) comparing location, scale, skewness and tails (row-wise). Legend
oracle is the distance between the non-private boxplot and the boxplot
corresponding to the true underlying distribution. Figure 5: Average error in
each metric with 95% confidence intervals between population boxplots and
estimated differentially private boxplots (y-axis) with $\epsilon=1$ and
increasing $n$ (x-axis); using boxplot estimation methods JointExp,
ApproxQuantile, unbounded, PrivateQuantile, and DPBoxplot (line colors); from
all sampled distributions (column-wise) comparing location, scale, skewness
and tails (row-wise). Legend oracle is the distance between the non-private
boxplot and the boxplot corresponding to the true underlying distribution.

###  B.2 Parameter tuning

We performed simulations under similar settings as in the previous section. We
varied $\lambda_{n}\in\\{n^{-1/2},n^{-1/4},1\\}$ and $c\in\\{1,20,40\\}$,
fixing $\epsilon=1$ and estimated differentially private boxplots with
DPBoxplot. . Figure 6 summarizes the results for simulated distributions
(column-wise), boxplot metrics (row-wise). Variations in $\lambda_{n}$ and in
$c$ are represented with different line colors and line styles, respectively.
We observe the method is not highly sensitive to the chosen parameters except
for skewness and tails. Based on this observation and to account for better
results in small samples we concluded that $\lambda_{n}=n^{-1/4}$ and $c=20$
were reasonable choices.

Figure 6: Average error in each metric with 95% confidence intervals between
population boxplots and estimated differentially private boxplots (y-axis)
with $\epsilon=1$ and increasing $n$ (x-axis); using DPBoxplot, with varying
$\lambda_{n}$ (color) and $c$ (line style); from all sampled distributions
(column-wise) comparing location, scale, skewness and tails (row-wise). Legend
oracle is the distance between the non-private boxplot and the boxplot
corresponding to the true underlying distribution.

###  B.3 Multiple boxplot estimation

This section presents an extensive simulation study to assess the
differentially private boxplot’s ability to describe, and facilitate
comparison of multiple samples. Specifically, we analyze whether the metrics
(quantified by distances on location, scale, skewness and tails as in Section
5) between differentially private boxplots across distinct datasets mirror the
corresponding pairwise distances observed between their population
counterparts. The similitude between two distinct differentially private
boxplots is anticipated to align closely with the similitude observed between
their respective population counterparts, thereby facilitating analogous
visual comparisons and interpretation. For the purpose of quantifying this
phenomenon, we introduce the concept of relative similitude between two
differentially private boxplots as follows:

$\tilde{d}\left(\tilde{B}(\mathbf{X}_{n},\epsilon),\tilde{B}(\mathbf{Y}_{m},\epsilon)\right)=\left|1-\frac{d\left(\tilde{B}(\mathbf{X}_{n},\epsilon),\tilde{B}(\mathbf{Y}_{m},\epsilon)\right)+1}{d\left(B(\nu_{X}),B(\nu_{Y})\right)+1}\right|.$

Here, distance is related to the four metrics previously described. In order
to empirically validate the consistency of our proposed approach with this
behavior, we conducted Monte Carlo simulations by generating $t$ datasets
$(\mathbf{z}_{1},\mathbf{z}_{2},\cdots,\mathbf{z}_{t})$, and two vectors
$\mathbf{m}$ and $\mathbf{s}$ of size $t$ sampling from uniform distributions
with interval $[-1,1]$ and $[0.5,2]$, respectively. Each vector
$\mathbf{z}_{i}$ is size $n_{i}$ for $i\in\\{1,2,\cdots t\\}$ where $n_{i}$ is
chosen randomly such that $n=n_{1}+n_{2}+\cdots n_{t}$. Here, $t$ plays the
role of the number of treatments in an experiment. We replicated this
simulation 1000 times for each combination of $t\in{4,8}$ and
$n\in\\{500,5000,50000\\}$ leading to datasets
$(\mathbf{X}_{1},\mathbf{X}_{2},\cdots,\mathbf{X}_{t})$ such that
$\mathbf{X}_{i}=s_{i}\mathbf{z}_{i}+m_{i}$. We performed this simulation
generating the initial $t$ datasets using mixture of the five distributions
described in the previous section such that each distribution generates the
same amount of vectors on each replica. For each dataset we calculated non-
private boxplots
$B(\mathbf{X}_{1}),B(\mathbf{X}_{2}),\cdots,B(\mathbf{X}_{t})$, and
differentially private boxplots using each of the quantile estimation methods
and for each $\epsilon\in\\{0.5,1,5,10\\}$. We then calculated pair-wise
relative boxplot distances for every pair of differentially private boxplots.

Figure 7 and 8 offers a detailed overview of the results obtained from various
scenarios considered. Each boxplot within these plots represents the average
pairwise distances observed among all differentially private boxplots
generated in a specific scenario. The x-axis on each plot denotes the sample
size ($n$). Variations in hue color correspond to different methods and
different privacy budgets ($\epsilon$), respectively. The first and last two
plots in the first row correspond to relative location and relative scale,
respectively. The first and last two plots in the second row correspond to
relative skewness and relative outliers, respectively. For each pair of
relative distance plots, the first and second visualization correspond to
$t=5$ and $t=10$, respectively. Figure 8 uses $\epsilon=1$ and Figure 7 shows
results with DPBoxplot method.

Figure 7 shows that DPBoxplot is consistently good among different scenarions.
All other methods have poor performance in at least one scenario. In Figure 8
relative distances exhibit a diminishing trend with augmented sample sizes and
$\epsilon$, consistent with expectations illustrated by the oracle boxplot. A
marginal increase in the parameter $t$ leads to a slight augmentation in the
relative distances, attributable to the partitioning of data into smaller
subsets and consequent reduction in sample sizes; however, this impact appears
to be negligible. Hence, our analysis suggests that our proposed methodology
effectively maintains the visual coherence of relative similarities among
diverse boxplots during the execution of multiple comparisons.

Figure 7: Average pair-wise relative boxplot distances (y-axis) between
multiple sample boxplots in terms of relative location, scale, skewness and
outliers, with $\epsilon=1$ using different boxplot estimations (color),
$n\in\\{500,5000,50000\\}$ (x-axis) and $t\in\\{5,10\\}$ (columns). Legend
oracle is the relative similitude of the non-private boxplots. Figure 8:
Average pair-wise relative boxplot distances (y-axis) between multiple private
boxplots in terms of relative location, scale, skewness and outliers, under
varying conditions of $\epsilon\in\\{0.5,1,5,10\\}$ (color),
$n\in\\{500,5000,50000\\}$ (x-axis); using DPBoxplot quantile estimation and
$t\in\\{5,10\\}$ (columns). Legend oracle is the relative similitude between
the non-private boxplots

###  B.4 Computational requirements

All simulations were conducted on a single CPU and did not require significant
computational resources. The execution time for the simulations presented in
the paper did not exceed 24 hours.

##  Appendix C Auxiliary results

The next lemma says that JointExp is inconsistent for $q_{n}=O(1/n)$.

######  Lemma 9.

For $-\infty<a<b<\infty$, if there exists $M>a$ such that $\nu(\left\\{X\leq
M\right\\})=0$, then for any $0<t<M-a$, we have that

$\Pr\left(|\tilde{\xi}_{q_{n}}-\xi_{q_{n},\nu}|\geq t\right)\geq
e^{-\epsilon\frac{nq_{n}}{2}}\frac{M-t-a}{b-a}.$

###### Proof.

We have that

$\Pr\left(|\tilde{\xi}_{q_{n}}-\xi_{q_{n},\nu}|\geq
t\right)\geq\Pr\left(\tilde{\xi}_{q_{n}}\leq
M-t\right)\geq\frac{\int_{a}^{M-t}e^{-\epsilon\frac{n|F_{\hat{\nu}}(x)-q_{n}|}{2}}dx}{\int_{a}^{b}e^{-\frac{n|F_{\hat{\nu}}(x)-q_{n}|}{2\epsilon}}dx}\geq
e^{-\epsilon\frac{nq_{n}}{2}}\frac{M-t-a}{b-a}.\qed$

In this case, if $q_{n}=n^{-k}$, then the sample complexity is then at least
$O((\log(1/\gamma\cdot\frac{M-t-a}{b-a})/\epsilon)^{1/(1-k)})$.

######  Lemma 10.

For all $n\geq 1$, $a,c\in\mathbb{R}$, we have that $n\geq a\log n+c$ if
$n\geq 2(a\log a+c)\vee a$.

###### Proof.

First, note that $h(n)=n-a\log n-c$ is increasing for $n\geq a$. Now, for
$n\geq a$ we can use the fact that $x-\log x\geq x/2$ for $x\geq 1$ to get
$h(n)\geq 0$ whenever $n\geq 2(a\log a+c)$. ∎

Let $\mathscr{B}$ denote the space of Borel functions from $\mathbb{R}^{d}$ to
$[0,1]$. For a family of functions, $\mathscr{F}\subseteq\mathscr{B}$, define
a pseudometric on $\mathcal{M}_{1}(\mathbb{R}^{d})$,
$d_{\mathscr{F}}(\mu,\nu)=\sup_{g\in\mathscr{F}}|\int gd(\mu-\nu)|$, where
$\mu,\nu\in\mathcal{M}_{1}(\mathbb{R}^{d})$.

######  Definition 2 (Ramsay et al., (2024)).

We say that $\phi(x,\nu)$ is _$(K,\mathscr{F})$ -regular_ if there exists a
class of functions $\mathscr{F}\subset\mathscr{B}$ such that $\phi(x,\cdot)$
is $K$-Lipschitz with respect to the $\mathscr{F}$-pseudometric uniformly in
$x$, i.e., for all $\mu,\nu\in\mathcal{M}_{1}(\mathbb{R}^{d})$

$\sup_{x}|\phi(x,\mu)-\phi(x,\nu)|\leq Kd_{\mathscr{F}}(\mu,\nu).$

##  Appendix D Broader impacts

Our proposed methodology inherits the extensive implications, both beneficial
and detrimental, of differential privacy. Differential privacy significantly
enhances data confidentiality by introducing noise into datasets, which
complicates the ability to link specific data points to individual identities.
This greatly enhances the privacy awarded to individuals whose information is
contained in such datasets. However, this method is not without its drawbacks.
The primary challenge lies in the trade-off between privacy protection and
data accuracy. The introduction of noise, while safeguarding privacy, may
distort the data, potentially yielding inaccurate insights or conclusions.
Such inaccuracies are especially problematic in contexts requiring high
precision, such as policy development or clinical research, where exact data
is crucial.