{'fraction_non_alphanumeric': 0.03483613161032516, 'fraction_numerical': 0.009023234829686443, 'mean_word_length': 4.397379406307977, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-parametric Differentially Private Confidence Intervals for the Median

Differential privacy is a restriction on data processing algorithms that provides strong confidentiality guarantees for individual records in the data. However, research on proper statistical inference, that is, research on properly quantifying the uncertainty of the (noisy) sample estimate regarding the true value in the population, is currently still limited. This paper proposes and evaluates several strategies to compute valid differentially private confidence intervals for the median. Instead of computing a differentially private point estimate and deriving its uncertainty, we directly estimate the interval bounds and discuss why this approach is superior if ensuring privacy is important. We also illustrate that addressing both sources of uncertainty--the error from sampling and the error from protecting the output--simultaneously should be preferred over simpler approaches that incorporate the uncertainty in a sequential fashion. We evaluate the performance of the different algorithms under various parameter settings in extensive simulation studies and demonstrate how the findings could be applied in practical settings using data from the 1940 Decennial Census.


Introduction
Statistical agencies constantly need to find the right balance between the two competing goals of disseminating useful information from their collected data and ensuring the confidentiality of the units included in the database. Many methods have been developed in the past decades to address this trade-off. However, with the advent of modern computing and the massive amounts of data collected every day, many of the data protection strategies commonly used at statistical agencies are no longer adequate to sufficiently protect the data [Abowd, 2018, Garfinkel et al., 2019. The problem's difficulty is amplified by the continual appearance of new data sources that facilitate attacks.
One promising strategy to circumvent this dilemma is to rely on formal privacy guarantees such as those provided by differential privacy (DP) [Dwork et al., 2006b]. These guarantees hold no matter what background knowledge a potential attacker might possess, or how much computational power they have. However, methodology for differential private statistical inference has mostly been studied from a theoretical perspective under asymptotic regimes. Although many algorithms have been proposed to ensure formal privacy guarantees for various estimation tasks, evaluations of their relative performance on real data with limited sample sizes and complex distributional properties are still limited, and only a small fraction of that literature has focused on inference and associated measures of uncertainty. Section 2.3 surveys related work.
In this paper, we address these issues, focusing on one of the key measures of location: the median. We chose the median for two reasons. On one hand, it is a widely used summary statistic for skewed variables such as income (see, for example, the U.S. Census Bureau's tables of median incomes for various subgroups of the population [U.S. Census Bureau, 2020a]). On the other hand, medians provide an interesting technical challenge for differentially private computation. The accuracy of differentially private median computations depends on the exact data distribution; as a result, providing sound and narrow confidence intervals appears to require releasing strictly more information about the data than is required for point estimation.
The discussion of confidence intervals is an important contribution of our paper. None of the previously proposed algorithms for DP median estimation come equipped with a method for additionally releasing DP uncertainty estimates on the point estimator. In fact, the level of uncertainty in the point estimate is typically data dependent, and hence measuring it requires additional privacy budget. Thus, the optimal algorithm for differentially private point estimates can be different from the optimal algorithm for differentially private confidence intervals. Instead of deriving the variance of some differentially private point estimate, we suggest estimating DP confidence intervals directly. We show that our proposed methodology ensures proper confidence interval coverage in a frequentist sense and discuss why this strategy requires less privacy budget than starting from the protected point estimates.
When designing and analysing differentially private algorithms it is tempting to separate the error due to sampling from the error due to privacy and bound the two separately. A main finding in our work is the limitation of this approach. We find that one can obtain considerably tighter confidence intervals by analysing the relationship between the two sources of error. Unlike approaches which treat the analysis of the non-private algorithm as a black-box, this involves looking at the different ways that the sampling error can result in the confidence interval failing to capture the median, and considering how the error due to privacy affects each of these modalities.
We assume simple random sampling throughout the paper. This assumption is often violated in survey practice. However, understanding the implications of complex sampling designs on the privacy guarantees is an open research problem [Drechsler, 2021] and we are not aware of any DP applications that take complex sampling designs into account. We see our contribution as an important first step towards the goal of better serving the needs of statistical agencies, while acknowledging the limitations of the current findings. We will come back to this point in the conclusions.
We evaluate several algorithms for computing valid differentially private confidence intervals. We discuss algorithms that satisfy two versions of differential privacy: the strictest version [Dwork et al., 2006b], now known as pure differential privacy, as well as a slight relaxation, concentrated differential privacy Steinke, 2016, Dwork andRothblum, 2016]. The focus of our paper is on empirical evaluation, using a mix of simulated and real data. Nevertheless, we found that new methodology and theory was also needed to adapt existing algorithms for confidence interval computation. We include an application using data from the U.S. Census 1940 to illustrate how statistical agencies willing to adopt the methodology could decide which algorithm and parameter settings to pick for their data release.
The algorithms we developed are all sound in the nonparametric, frequentist sense: when run with nominal coverage 1 − α the probability that the true population median is contained in the computed confidence interval is at least 1 − α, where the probability is taken over the entire process of sampling from the population and computing the private confidence intervals based on the drawn sample. Since all algorithms rely on non-parametric strategies for computing the confidence intervals, the intervals are valid for every IID distribution on observations. We summarize our findings briefly: • In our comparison of several algorithms, the best choice across a range of settings was a variant of the exponential mechanism [McSherry and Talwar, 2007], a generic framework for DP algorithm design that we adapt for confidence interval estimation. This algorithm is tailored to the median, and releases only a single confidence interval.
• A different algorithm, based on a differentially private CDF estimate [Li et al., 2010], consistently produced confidence intervals that were slightly wider than those of the exponential mechanism. However, the algorithm's output can be used to produce a confidence interval for any quantile of the data set or even a confidence band for the entire CDF. In principle, the approximation to the entire CDF would also allow the incorporation of a parametric model or a Bayesian prior after the fact. Its flexibility makes it a better choice for settings where eventual users will be interested in more than a single median.
• For settings where only a very loose bound on the range of the data range is known a priori, a hybrid algorithm that uses binary search to narrow the range and then switches to the CDF-based estimator produces narrower intervals than other methods.
• The methods we tested exhibited noticeable bias that depends on the underlying distribution and appears to be hard to correct. The bias was low relative to the width of the confidence intervals, and so would not be an issue for one-shot applications. However, it might be a concern when aggregating estimates across many small areas. It is not clear whether bias is necessary for accurate nonparametric DP median approximations. 1 The remainder of the paper is organized as follows: In Section 2 we review some of the privacy definitions that are relevant for this paper and discuss confidence interval estimation for the median without privacy considerations. We extend these discussions to differentially private confidence intervals in Section 3. Section 4 contains a highlevel review of the algorithms we considered (detailed descriptions of the different algorithms can be found in the Appendix). In Section 5 we present the results from extensive simulation studies that evaluate the performance of the algorithms under various parameter settings. Section 6 illustrates how the methodology could be applied in practice by replicating one of the income tables published by the U.S. Census Bureau using publicly available data from the 1940 U.S. Census. The paper concludes with some final remarks.

Differential Privacy
The algorithms in this paper satisfy a version of differential privacy (DP) called concentrated differential privacy (CDP). This notion of privacy lies between the more common notions of pure differential privacy and approximate differential privacy. Since our algorithms often include hyperparameters, we state a definition of DP for algorithms that take as input not only the dataset, but also the desired privacy parameters and any required hyperparameters. Let X be a data universe (e.g., R for medians) and X n be the space of datasets of size n. Two datasets d, d ∈ X n are neighboring, denoted d ∼ d , if they differ on a single record. Let H be the space of hyperparameters and Y be an output space. In order to build some intuition, let us first define pure and approximate DP.
The key intuition for this definition is that the distribution of outputs on input dataset d is almost indistinguishable from the distribution on outputs on input dataset d . Therefore, given the output of a differentially private mechanism, it is impossible to confidently determine whether the input dataset was d or d . If δ = 0, then we refer to this as ε-pure differential privacy. If δ > 0, we refer to (ε, δ )-approximate differential privacy. For strong privacy guarantees, the privacy-loss parameter is typically taken to be a small constant less than 1 (note that e ε ≈ 1 + ε as ε → 0). However, in practice, larger values of ε are occasionally used to satisfy utility constraints while providing some level of non-trivial privacy guarantee.
Concentrated differential privacy has the same intuition; it bounds the divergence between the distributions M(d) and M(d ).
Figure 1: A graphical representation of the process of computing a confidence interval. When privacy is not a concern, no restrictions are placed on the function M. When computing a differentially private confidence interval, we require that M is differentially private. The probability is taken over all the randomness in the system, both the randomness due to sampling, and the randomness in M.
Definition 2.2 (ρ-Concentrated Differential Privacy [Bun and Steinke, 2016]). Given ρ ≥ 0, a randomized mechanism M : X n × H → Y satisfies ρ-concentrated differential privacy if for all datasets d ∼ d ∈ X n , hyperparams ∈ H , and α ∈ (1, ∞), where D α is the α-Rényi divergence and the probabilities are taken over the random coins of M.
In order to give some intuition for concentrated DP, let us elaborate more on its relationship with pure and approximate DP. Given data sets d ∼ d , and a randomised mechanism M, we can define a random variable, called the privacy loss random variable, denoted Z = Priv(M(d), M(d )), as follows. Let y ∼ M(d) (i.e. y is the output of the mechanism M on input d), then Z = ln Pr(M(d))=y Pr(M(d ))=y . Then M is ε-pure differentially private if and only if Pr(|Z| > ε) = 0, and M being (ε, δ )-approximately differentially private is (almost) captured by the requirement that Pr(|Z| > ε) ≤ δ . Now, ρ-concentrated differential privacy essentially translates to the requirement that Z is a subgaussian random variable with mean ρ and variance 2ρ. From this perspective, it is clear that concentrated differential privacy lies between pure and approximate DP.
We will focus in this paper on algorithms that satisfy concentrated differential privacy. While still satisfying a rigorous notion of privacy, this will allow our algorithms to be significantly more accurate than their corresponding purely differentially private counterparts. For most of our algorithms little accuracy is gained from transitioning to approximate differential privacy. Additionally, CDP has the desirable property of being a one-parameter property, which allows for simpler privacy accounting. The lemma below captures the fact that the class of ρ-CDP algorithms is closed under adaptive composition and post-processing.
If M is ρ-CDP and M is ρ -CDP then M is ρ + ρ -CDP.

Confidence Intervals for the Median
In many statistical applications, we assume that data are drawn i.i.d. from an underlying population distribution, and the statistic of interest is a property of the underlying population. However, one typically only has access to a sample from that population, so the statistic computed on the sample is used as an estimate of the true population statistic. We will refer to the median of the underlying population as the population median and the median of a given sample as the sample median. Since there is randomness in the sampling process, there is always uncertainty in how well the sample median matches the true population median. As this uncertainty can be large, sample statistics should be accompanied by a measure of the uncertainty. Providing a measure of uncertainty is even more important for differentially private statistics since randomness in the algorithm provides an additional source of uncertainty.
One method for capturing the uncertainty in an estimate is a confidence interval. We consider the standard set-up for statistical inference. Let P ⊂ ∆(R) be the set of possible population distributions over the data domain R. For any P ∈ P, a median of P is defined to be any value m such that m −∞ P(x)dx ≥ 1/2 and ∞ m P(x)dx ≥ 1/2.
For every distribution the set of medians is a non-empty, compact set. Since defining a convention here will be convenient, we will refer to the midpoint of the set of medians as the median, denoted med(P). Let I R be the set of intervals in R. Given n ≥ 0, let M : X n → I R be a randomised mechanism that takes as input a data set of size n and outputs an interval in R. Given a desired confidence level 1 − α, the goal of M is to, with probability 1 − α, output an interval that contains med(P).
Definition 2.3. For any α ∈ [0, 1] and n ∈ N, M : X n → I R is a (1 − α)-confidence interval for the median for P if for all P ∈ P, where the randomness is taken over both the randomness M and the randomness in the sample d ∼ P n .
A graphical representation of the framework for computing a confidence interval is given in Figure 1. We will refer to 1 − α as the coverage of the confidence interval.
When one is not concerned with privacy, a non-parametric confidence interval for the median can be computed using the order statistics of the sample. The rank of the median med(P) in a data set d ∈ P n is distributed as the binomial Bin(n, β ), for some β = Pr P (x < med(P)). We can exploit this to obtain a confidence interval for the median. For a data set d ∈ R n , let d (k) denote the k-th smallest value in d, referred to as the k-th order statistic. The median of d is the midpoint of d ( n/2 ) and d ( n/2 ) . Lemma 2.3 (Non-private (1−α)-confidence interval). Let C Bin be the CDF of the binomial random variable Bin(n, 1/2) and let N α For any data set d ∈ R n , let is a (1 − α)-confidence interval for the median for ∆(R).
We will often refer to the interval, [ci α L (d), ci α U (d)], output by the mechanism in Lemma 2.3 as the non-private (1 − α)-confidence interval for the median, but we note that it is the mechanism M that satisfies Definition 2.3, not the output.
In this vein, the goal of CDP confidence intervals is not to privately estimate the specific interval [ci α L (d), ci α U (d)], but to output valid confidence intervals. These confidence intervals may, or may not, contain [ci α L (d), ci α U (d)]. Referring to Figure 1, producing a CDP confidence interval involves the same procedure, with the additional requirement that M is CDP. That is, the realised sample d is only accessed through a CDP mechanism.
The definition of confidence intervals as stated in Definition 2.3 and Definition 2.4 only requires that the confidence interval is valid for distributions P ∈ P. Parametric estimation is when one defines P to be only distributions of a particular, often quite simple, form. For example, P might be the set of all log-normal distributions over R. In non-parametric estimation, we assume no knowledge of the underlying population and set P = ∆(R), the set of all distributions over R. If one has accurate knowledge of the underlying population, then parametric estimation can result in a tighter confidence interval. However, if there is model mismatch (for example if the underlying population is not exactly log-normal) then parametric estimation can result in invalid confidence intervals. This effect can be amplified by private algorithms which may rely on the modeling assumptions in non-trivial ways. In this paper our goal is to focus on non-parametric confidence intervals, which means our algorithms will always produce valid confidence intervals. As a minor caveat, we restrict ourselves to the set of distributions with continuous probability density functions on R. 2 Denote the set of all continuous distributions on R by ∆ C (R).
Note that a confidence interval does not directly output a point estimate for the median itself. In the absence of privacy constraints, one can simply additionally release the sample median med(d). However, under privacy constraints, it is typically desirable to compute as few statistics as possible, in order to allocate the maximum amount of privacy budget to each statistic. As such, rather than allocating some of the privacy budget to providing a point estimate of the median, it is often preferable to allocate the entire budget to estimating the confidence interval, then use the midpoint of that interval as a point estimate of the median.

Related Work
Computing confidence intervals for the median is one of the most fundamental statistical tasks. However, finding a differentially private estimator for this task that is accurate across a range of datasets and parameter regimes is surprisingly nuanced. There has been a significant amount of prior work on differentially private point estimators for the median [Nissim et al., 2007, Bun and Steinke, 2019, Asi and Duchi, 2020, Alabi et al., 2020, Tzamos et al., 2020 and other quantiles [Gillenwater et al., 2021]. To the best of our knowledge, none of these works addressed DP confidence intervals for the median. However, there has been significant work on DP confidence intervals for other estimation tasks like (Gaussian or sub-Gaussian) mean estimation [Karwa and Vadhan, 2018, Gaboardi et al., 2019, Du et al., 2020, Biswas et al., 2020, and linear regression [Barrientos et al., 2017, Evans and. There are also several works on designing more general DP confidence intervals using bootstrapping, or a technique called subsample-and-aggregate [Nissim et al., 2007], to account for the combined uncertainty from sampling and noise due to privacy [Barrientos et al., 2017, Ferrando et al., 2020, Brawner and Honaker, 2018, D'Orazio et al., 2015. These algorithms typically require a parametric model on the data or a normality assumption on the quantity being estimated; neither hold in our setting.
To the best of our knowledge, our work is unique in focusing on valid non-parametric differentially private confidence intervals for the median. This approach allows us to define algorithms that provide accurate and private confidence intervals without requiring distributional assumptions on the underlying population.
3 Designing DP Confidence Intervals

Roadblocks and first attempts
There are several roadblocks in designing CDP confidence intervals for the median. Firstly, CDP algorithms that estimate the median using data independent output perturbation methods (methods that involve simply adding noise to the non-private estimate) necessarily perform poorly since the median is very sensitive for worst-case data sets. Thus, in order to design algorithms that perform well on "typical" data sets, the noise addition must be data dependent. This creates difficulties when releasing information regarding the uncertainty in the private estimate since the uncertainty itself might reveal sensitive information. In order to explore these roadblocks in more detail, let us first consider the simpler task of designing a CDP point estimator for the median. A first attempt may be to consider the global sensitivity [Dwork et al., 2006b]: Definition 3.1 (Global Sensitivity). For a query f : X n → R, the global sensitivity is For any function f , one can create a differentially private mechanism by adding noise proportional to GS f / √ ρ. If one has no bound on the data, then GS med is infinite. Even if one knows that all the data lie in a bounded range [a, b], GS med = |b − a|, and adding noise proportional to |b − a|/ √ ρ essentially removes the signal for any reasonable value of ρ.
However, for the type of datasets that we typically see in practice, changing one data point, or even a few data points, does not result in a major change in the median. For such data sets, one might consider the local sensitivity [Nissim et al., 2007], which can be substantially smaller than the global sensitivity.
Definition 3.2 (Local Sensitivity [Nissim et al., 2007]). The local sensitivity of a query f : Unfortunately, since the local sensitivity itself is data dependent, adding noise proportional to the local sensitivity is not differentially private. Several approaches have been explored in the DP literature for adding noise that is close to the local sensitivity, or at least significantly less than the global sensitivity for typical data sets. In [Nissim et al., 2007], Nissim et al. define the smooth sensitivity, SS f , a smooth upper bound on the local sensitivity such that adding noise proportional to SS f / √ ρ is differentially private. They showed that for many statistics, including the median, the smooth sensitivity can be much smaller than the global sensitivity on typical data sets. Several other DP mechanisms for median estimation have been proposed that avoid the large GS med by calibrating the noise introduced to the specific data set [Nissim et al., 2007, Dwork andLei, 2009]. While these estimators can perform well as point estimators for the median, our goal is not just to provide an estimate of the median, but also to quantify the uncertainty in our estimate. In algorithms like the smooth sensitivity mechanism that tailor the noise to the specific data set, the uncertainty itself is data dependent and thus sensitive. The task of differentially privately releasing an estimate of this uncertainty is nontrivial. Even if one could release a DP estimate of the amount of noise added to the non-private median, the uncertainty in the non-private median is still unaccounted for. For this reason, we focus on DP algorithms that attempt to directly estimate the confidence interval.

Accounting for all sources of randomness
Accurate and tight coverage analysis is a crucial component of designing good algorithms since overly conservative coverage estimates can result in confidence intervals that are wider than necessary. Valid differentially private confidence intervals need to account for two sources of error; sampling error and error due to privacy. Sampling error, also present in the non-private context, captures how well the realised sample d represents the underlying population P. The error due to privacy takes into account the additional randomness in M as a result of the privacy guarantee. Our experimental results highlight that it is important to carefully exploit the dependence between the two sources of randomness.
As a primer, let us first consider the coverage analysis of the non-private algorithm described in Lemma 2.3. This coverage analysis relies on the fact that if P is continuous then for all m ∈ n, Pr(rank d (med(P)) = m) = Pr(Bin(n, 1/2) = m).
There are two ways that the interval [ci α L (d), ci α U (d)] can fail to capture med(P); med(P) < ci α L (d) or med(P) > ci α U (d). Let us focus on the probability of the first type of failure, med(P) < ci α L (d). For every P ∈ ∆ C (R), where C Bin is the CDF of the binomial random variable Bin(n, 1/2). The probability of failure at the upper end of the confidence interval is analogous. . A naive way to analyse the coverage error of M is to attempt to find β 1 and β 2 such that assuming β 1 denotes the failure probability for the non-private confidence interval, the M(d) contains the non-private interval [ci β 1 L , ci β 1 U ] with probability 1 − β 2 . Then M has coverage at least 1 − (β 1 + β 2 ). Even if β 1 and β 2 are chosen carefully, this analysis can be overly conservative. In particular, it assumes that the only way that M(d) can succeed in containing med(P) is if both med(P) ∈ [ci Figure 2, this corresponds to ensuring that M(d) contains the red dotted interval with high probability. It's clear from this figure that neither of these events are necessary.
A more careful analysis of the relationship between the sampling error and the error due to privacy results in a tighter coverage analysis. As in the non-private setting, there are two ways that M(d) can fail to contain med(P), and we will focus on analysing the probability that med(P) < M(d) L Now, we have reduced the problem to analysing the failure probability conditioned on the empirical rank of med(P) in the data set d. This is a helpful reduction since, as we will see in the following section, most of our algorithms will come with accuracy guarantees on the rank of M(d) L . Accuracy guarantees of this form can then be exploited, via Equation (1), to obtain a coverage analysis of M.
Our experiments show a stark difference between the performance of algorithms designed using the naive analysis, and those using the tighter, more careful analysis. In Figure 7 we directly compare the confidence intervals that arise from the different analyses. This highlights the importance of understanding the relationship between the two sources of error.

Algorithms
In this section we will introduce the four algorithms for releasing CDP confidence intervals for the median that will be the focus of this paper; ExpMech, CDFPostProcess, NoisyBinSearch, and BinSearch + CDF. The first algorithm, which we call ExpMech, is based on the exponential mechanism. This mechanism is efficient, satisfies the stronger privacy guarantee of pure differential privacy and outputs the tightest, or close to the tightest confidence intervals in a majority of parameter regimes we studied. The remaining three algorithms partially address a common frustration with differentially private data analysis; that exploratory data analysis to visualise the data set and verify findings typically requires additional privacy budget. For many tasks, this means allocating privacy budget away from the primary task resulting in a noisier algorithm. A key feature of the three algorithms CDFPostProcess, NoisyBinSearch and BinSearch + CDF is that they release additional information about the data set without consuming additional budget. In particular, CDFPostProcess releases a full CDP estimate to the empirical CDF. It is then notable, and perhaps surprising, that in many settings these algorithms perform almost as well as ExpMech, which releases no side information.
In this section we will give a high level description of each algorithm. Further information for all algorithms, including pseudo-code and proofs of the privacy and validity guarantees, can be found in the online supplement accompanying this paper. Real code is available in our GitHub repository. 3 We note that we also experimented with several other algorithms that are not discussed in this section. Brief descriptions of these additional algorithms can be found in the online supplement F, but we do not focus on them here since they are outperformed by other algorithms in every parameter regime we studied.

Confidence intervals based on exponential mechanism, ExpMech
Our first private mechanism is an instantiation of the exponential mechanism [McSherry and Talwar, 2007], a differentially private algorithm designed for general optimization problems. The exponential mechanism has been used in prior work to give DP point estimates for the median [Dwork and Lei, 2009, Thakurta and Smith, 2013, Johnson and Shmatikov, 2013, Alabi et al., 2020, Asi and Duchi, 2020. Our extension to providing confidence intervals for the median, while using similar ideas to prior work, requires a careful coverage analysis that is new to this work.
The exponential mechanism is defined with respect to a utility function u, which maps (data set, output) pairs to real values. For a data set d, the mechanism aims to output a value r that maximizes u(d, r) by sampling a value r with probability proportional to e εu(d, Recall that the rank of a value r ∈ R in a data set d, denoted rank d (r), is the number of data points in d that are less than or equal to r. For any k ∈ [n], one way to instantiate the exponential mechanism to compute the k-th order statistic is by using the following utility function. Let In practice we will use a slight variant of the exponential mechanism described in online supplement B, which uses a granularity parameter θ . This version of the exponential mechanism has the guarantee that with high probability it will output a value within θ of some r, such that rank d (r) is close to k. Let us denote the exponential mechanism with privacy parameter ε for the k-th order statistic by A ε k . While it is tempting to use the exponential mechanism to release DP estimates of ci α , this will not create a DP α-confidence interval. This is because the private algorithm may under or over estimate For each x ∈ R discrete , the red dot is atC(x), the yellow dot is atĈ(x), and the pink dot is at the upper and lower thresholds a L x and a U x . For each x the blue vertical line represents a point wise 99% confidence interval onĈ(x) based on the measurementC(x). The purple dashed vertical lines represent the DP interval CDFPostProcess(d).
these quantities resulting in an invalid confidence interval. Instead we need to choose k L and k U carefully so that In order to analyse the coverage, we return to Equation (1), In this computation, we use the inequality Pr(med( k L outputs something close to d (k L ) (which is the goal of this algorithm) with high probability. Some additional performance could be obtained by tighter analysis of these terms, but we expect the improvement to be small. Now the first term is only due to sampling error and easily bounded. The second term depends on both the sampling error and the error rate of the private algorithm A ε/2 k L . In Appendix B, we show how to obtain a distribution independent upper bound on the error rate Pr( . This allows us to give an upper bound for Pr(med(P) ≤ A ε/2 k L (d) − θ ) which holds for any distribution P ∈ ∆ C (R). Given this coverage analysis, we can search for k L and k U that maximise Pr(med(P) / ∈ M k L ,k U (d)) subject to the constraint that Pr(med(P) / ∈ M k L ,k U (d)) ≤ α. We will refer to this algorithm as ExpMech. Pseudo-code and details of the analysis can be found in online supplement B. A graph representing the distribution of ExpMech on a single data set can be found in Figure 3.

Confidence intervals based on CDF estimator, CDFPostProcess
Our second CDP confidence interval CDFPostProcess is obtained from post-processing the output of a CDP cumulative distribution function (CDF) estimator. Unlike ExpMech, which only released the confidence interval, CDFPostProcess can additionally release a CDP estimate of the CDF without consuming additional privacy budget. There has been considerable work in the DP literature on DP cumulative distribution function (CDF) estimators, both for the parametric and non-parametric models [Diakonikolas et al., 2015, Brunel andAvella-Medina, 2020]. We will focus on a particular CDF estimator based on the tree-based mechanism introduced in Li et al. [2010], Dwork et al. [2010], and Chan et al. [2011]. This mechanism was further refined in Honaker [2015], whose algorithm we base our mechanism on. Given a range for the data [r , r u ], and a discretization of this range R discrete = [r , r + θ , r + 2θ , · · · , r u ], the algorithm in Honaker [2015] outputs a DP estimate,C, of the empirical CDF,Ĉ, of a data set restricted to R discrete . The relevant feature of this estimator for our analysis is that we understand the marginal distribution ofĈ(x) for each x ∈ R discrete . Specifically, for all x ∈ R discrete , there exists σ x > 0 such that the value we release is equal to the empirical CDF with normally distributed noise, that is: The noise values are not independent across different values of x, nor are they perfectly correlated. The estimatesC(x) don't even increase monotonically with x (see Figure 4). Thus, we seek a procedure that uses only our knowledge of the marginal distribution. The main observation is that we can rewrite the distribution ofC(x) in terms of the true distribution CDF value, C(x): We will use this, for each value x in R discrete , to test the hypothesis that x is less than the median. Observe that which we can easily numerically evaluate. Thus, for each x, there exists a value a x (that just depends on σ x ) such that if x < med(P) then Pr(C(x) > a x ) ≤ α/2, where 1 − α is our desired coverage. Equivalently, ifC(x) > a x then we can be confident that x ≥ med(P). The upper end of our confidence interval is then defined as the grid point just to the right of the largest value x for which the test accepts, that is, To see why this leads to a valid confidence interval, let x * be the largest value in R discrete that is less than med(P). With probability at least 1 − α 2 , the test at x * accepts. In that case, ci α U (d) will be a grid point greater than x * , and thus we will have ci α U (d) ≥ med(P). The process and reasoning for the left end of the interval are symmetric. Pseudo-code and a proof that this process gives valid confidence intervals can be found in Appendix D. We will refer to this algorithm as CDFPostProcess. A graph representing a single run of this algorithm can be found in Figure 4. In particular, note that the horizontal pink dotted lines represent the values of a x , which are closer to N α L and N α U when σ x is small. 4

Confidence intervals based on noisy binary search, NoisyBinSearch
One draw-back of the CDF-based estimator CDFPostProcess is that it spends its privacy budget roughly evenly across the entire range R. This can result in substantially reduced performance if the data is concentrated in a small subset of the range. We can see this effect in Figure 4 where the values of theC below x = 0 are very unlikely to impact the confidence interval. NoisyBinSearch attempts to locate the data within the region R using a small amount of the privacy budget then delves more deeply into the region actually containing the data. This defines a CDP confidence interval in its own right but we will see the real power of it in the design of the next algorithm.
NoisyBinSearch uses noisy queriesĈ(x)+N (0, σ 2 ) to the empirical CDF to search for the relevant quantiles using binary search. Unlike CDFPostProcess that obtains a CDP estimate toĈ over the entire R discrete , NoisyBinSearch minimises the number of values we need to evaluateĈ on. Full details and pseudo-code are given in online supplement C. As in the design of ExpMech, this algorithm starts with two target quantiles k L and k U which are defined to ensure that the resulting confidence interval has the desired coverage. We will discuss this choice in online supplement C. LetĈ be the empirical CDF. In the non-private setting, we can find an approximation to d (k L ) by iteratively asking queries of the form "isĈ(x) ≤ k L ?" and adjusting our search accordingly. This search method is called binary search and minimises the number of values we need to evaluateĈ on. We can perform a CDP version of binary search using noisy queries of the form "isĈ(x) + N (0, σ 2 ) ≤ k L ?" In setting the variance σ we need to balance the desire for the noisy query to answer correctly (so the binary search moves in the right direction) and adding enough noise to guarantee privacy. Notice that if |Ĉ(x) − k L | is large, for example if x is inside R but far from the concentration of the data, then σ can be chosen to be quite large while still ensuring the response is correct with high probability. A single query of this form is 1 nσ -CDP [Bun and Steinke, 2016] so if |Ĉ(x) − k L | is large, then we only need to consume a small amount of privacy budget in order to ensure the binary search moves in the right direction. The privacy guarantee for the entire algorithm is sum of the privacy guarantees for each iterate. Since a priori we don't know |Ĉ(x) − k L |, at each iteration we'll start with a quite noisy query, then keep decreasing the noise until we are confident which direction to move in. We repeat until the privacy budget is consumed. We search for both the upper and lower end points of the confidence interval using noisy binary search then a final post-processing step on the noisy measurements ensures we release a valid confidence interval. Full details and pseudo-code are given in online supplement C. At the points x i at which NoisyBinSearch obtains noisy measurements, the noisy measurementĈ(x i ) + N (0, σ 2 ) can be released in addition to the confidence interval without consuming additional privacy budget. While this is not as informative as the full CDP CDF released using CDFPostProcess, it does provide useful additional information about the distribution. We note that unlike ExpMech and CDFPostProcess, our analysis of NoisyBinSearch separates the two sources of randomness coming from sampling and noise added for privacy. An interesting open problem is whether this analysis can be improved by considering the relationship between the two sources of randomness.

Range-robust estimator based on CDF estimator, BinSearch + CDF
The mechanism NoisyBinSearch was designed to improve the performance of CDFPostProcess when the data is concentrated in a small region within R. That is, when NoisyBinSearch can greatly reduce the search range for the median using a small amount of privacy budget. However, when the data is well spread out within the range, we do not expect NoisyBinSearch to perform as well as CDFPostProcess which is highly optimised to provide an accurate approximation to the CDF, C. Our final algorithm, BinSearch + CDF first uses NoisyBinSearch with a portion of its privacy budget to narrow down the search space R. It then clips the data to this range, and runs CDFPostProcess with the remaining privacy budget within this smaller range to obtain a confidence interval. Full details and validity proofs are given in online supplement E.
The privacy budget and coverage α both need to be partitioned between the two stages of the algorithm. We expect the optimal split to be distribution dependent. In particular, it likely depends on how large a region the data occupies within the range R. We found experimentally, for the parameter regimes we studied, using ρ/4 for the first step, and 3ρ/4 for the second step seemed to be a good choice. Similarly, we ensure that the region found in the first setp contains the median with probability 1 − α/4, and the second step finds a 1 − 3α/4-confidence interval within that region.A representation of a single run of BinSearch + CDF is shown in Figure 5. Notice that relative to CDFPostProcess in Figure 4, it uses only a few measurements to narrow into the range of interest, about [0.8, 3] then takes more measurements within this range.

A note on hyperparameters
All our algorithms require some hyperparameter tuning and domain knowledge. Throughout this manuscript we will attempt to give guidance on how one may set these parameters and how sensitive the performance is to these choices.
All our algorithms require as input a range space for the median. That is, an interval R ⊂ R that is promised to contain med(P). This should be chosen as small as is reasonable, but one can typically be quite conservative when defining R. We expect the dependence of all algorithms to be approximately log |R|, we see this explicitly for ExpMech in Lemma B.2 in the online supplement. We will empirically explore the dependence on |R| in Figure 6d below. Note that requiring a range on the median is different to requiring a range on the data points themselves. It is a preferable condition since the data points may occupy a considerably larger range than R, and, in practice, guaranteeing a bound on outliers in the data may be difficult. The fact that a bound on the median suffices comes from the fact that if R contains the median then projecting the data into R leaves the median unchanged.
All the algorithms we consider require an additional hyperparameter we refer to as the "granularity" parameter, θ . In order to gain some intuition, one can imagine each algorithm as discretizing the range R so that any two potential output points are θ apart. This intuition is not exact, please refer to the online supplement for more details. This granularity parameter affects each algorithm differently. ExpMech is not very sensitive to this parameter in general, and it can typically be taken to be very small. In fact, setting this parameter to 0 is a reasonable choice in a wide variety of parameter settings. CDFPostProcess and BinSearch + CDF can be sensitive to this parameter. In all algorithms, the width of the confidence interval will be at least θ .

Simulation Studies
In this section we present extensive simulation studies to evaluate the different algorithms under various parameter settings. 5 We focus on log-normal data, as this is the natural use-case for computing a non-private median and corresponding confidence interval. We expect many of our findings to extend to other types of skewed data. We evaluate the performance of the four algorithms described in Section 4, as well as the non-private confidence interval described in Lemma 2.3, in terms of width of confidence interval, coverage, and bias.

Data description
In order to visualize the distribution of the noisy confidence intervals, we run each private algorithm 5 times on 100 independently drawn datasets. Let x 1 , . . . , x 100 , each contain n = 1000 i.i.d. draws from the underlying log-normal distribution. The underlying normal random variable has mean µ = ln(1.5) and standard deviation of either σ = 1 or σ = 5. We run each DP algorithm on each x i for 5 trials. In our experiments, we show the relative performance of the algorithms as we vary n, ρ, σ , and |R|.

Utility measures
We consider two main utility measures in our experiments. The first measure is the relative width of the CDP confidence interval M(d) compared to the width of the non-private confidence interval, [ci α L , ci α U ]. For a data set d ∈ D n , α ∈ [0, 1], and interval I = [I L , I U ], the relative width is defined as We are concerned with the distribution of rel-width α (d, I) when I = M(d). We expect the private confidence intervals to be wider than the non-private confidence intervals, so if I = M(d) then we expect rel-width α (d) ≥ 1 with high probability. If rel-width α (d) ≤ 2, then, intuitively, the additional uncertainty due to privacy is less than the uncertainty due to sampling. We are interested in the distribution of the relative width over multiple trials, so for each algorithm we show box plots of this metric over 500 trials (100 datasets times 5 trials of the DP mechanism on each). The second utility measure is empirical coverage of the DP confidence interval. For intervals I 1 , · · · , I T and distribution P, let cov T (P, I 1 , · · · , Given n ∈ N and a 1 − α-confidence interval M, for all t ∈ [T ], let d t ∼ P n and I t = M(d t ) then cov T,n,P (M) = cov T (P, I 1 , · · · , I T ) gives an estimate of the empirical coverage of M on the distribution P. We estimate the coverage over 5,000 trials (1,000 drawn samples of size n = 1, 000 times 5 trials of the DP mechanisms on each). A key component of the confidence intervals presented in the previous section is that the coverage should be at least 1 − α. However, the empirical coverage may, and in many settings will, exceed 1 − α.

Results and Discussion
Comparison Among Algorithms Figure 6 demonstrates the performance of our four CDP confidence interval algorithms across a range of parameter regimes on log-normal data, in terms of the relative width metric. Notice that in a variety of regimes, including large n, large ρ and large σ d all of the CDP algorithms provide confidence intervals that are at most twice the width of the non-private confidence interval with high probability. Our results indicate that ExpMech provides the tightest, or close to the tightest, confidence intervals in most parameter regimes we studied. This algorithm is the most targeted of the CDP algorithms we discuss and is carefully calibrated to not waste privacy budget on estimating additional information about the underlying distribution. It is a good general choice when one is solely interested in confidence intervals for the median. There are a few regimes in which the other algorithms outperform ExpMech which we will discuss in this section. The CDFPostProcess algorithm is appealing in practice since it allows a CDP estimate of the CDF to be released without consuming additional privacy budget. This can be used not only to produce confidence intervals for the median, but to produce a wealth of other insights about the distribution P. Surprisingly, in a variety of parameter regimes, CDFPostProcess provides confidence intervals that are almost as tight as those obtained by ExpMech. In fact, when σ d is large, CDFPostProcess can result in tighter confidence intervals than ExpMech (Figure 6b). We explore this further in Appendix G. Conversely, when σ d is small, or |R| is large, CDFPostProcess is not a good  choice. As discussed earlier, these are regimes where CDFPostProcess spends a lot of its privacy budget estimating the CDF in regions that are far from the median. The remaining CDP confidence intervals BinSearch + CDF and NoisyBinSearch were designed to outperform CDFPostProcess when the data is concentrated in a small subset of R. This setting is shown in Figure 6b, when the data is very concentrated and Figure 6d when |R| is large. In both settings we see either BinSearch + CDF or NoisyBinSearch outperforming the other CDP algorithms. BinSearch + CDF involves a hyperparameter that decides how much of the privacy and coverage budget is spend on each step of the algorithm. In our experiments, we used one quarter of the privacy budget to perform the range-finding step (using NoisyBinSearch), and three quarters to find the confidence interval within that region (using CDFPostProcess). This choice was made experimentally by testing multiple choices for this partition. Since this parameter interpolates between CDFPostProcess and NoisyBinSearch, the optimal choice appears to be data dependent. In particular, we conjecture that the smaller the fraction of the range that is occupied by the data, the higher the fraction of the budget that should be allocated to the NoisyBinSearch step. Of course, one typically does not know a priori how concentrated the data is, so the (1/4, 3/4) seems to be a reasonably safe split that performs well in a variety of contexts.

Empirical coverage analysis.
A key component of the algorithmic design of each of the CDP confidence intervals was the coverage analysis. We discussed in Section 3.2 how a careful coverage analysis that leverages the relationship between the two sources of the randomness in the CDP confidence intervals potentially results in a much tighter coverage analysis than the naive analysis that separates the sources of randomness. Our experimental results presented in Figure 7 highlight two key findings regarding the coverage; that the careful analysis does result in substantially tighter intervals, and that the empirical coverage of the CDP confidence intervals is still notably above the target coverage. Figure 7a compares the relative width of the different confidence intervals. ExpMechUnion and CDFPostProcessUnion refer to the versions of ExpMech and CDFPostProcess resulting from the naive coverage analysis. While the relative width is only slightly reduced for ExpMech, the relative width of CDFPostProcess is almost halved if the improved approach is used to produce the confidence intervals. These findings are also reflected in Figure 7b, which compares the empirical coverage of the naive coverage analyses of ExpMechUnion and CDFPostProcessUnion and the more careful analyses described in Section 4. While theoretically we can show that the careful analysis will result in empirical coverage that is much closer to the target coverage, Figure 7b shows that this improvement is practically relevant. While the improved analysis only leads to modest reduction in the overcoverage for ExpMech, the changes for CDFPostProcess are more substantial. The naive approach results in coverage rates that are close to 1 irrespective of the selected α. The improved approach leads to coverage rates that are much closer to the nominal coverage rates. This highlights the importance of designing algorithms that consider the relationship between the two sources of randomness.
Despite the substantial improvement, Figure 8a shows that all the CDP algorithms exhibit empirical coverage higher than the target coverage for moderate values for ρ. As expected, Figure 8b reveals that the coverage improves with ρ as each algorithm trends towards outputting the non-private confidence interval [ci α L , ci α U ] when ρ = ∞. Overcoverage does not necessary correspond to substantially larger confidence intervals. We see in Figure 6 that in a wide range of parameter regimes our CDP algorithms still result in confidence intervals that are at most twice as wide as their non-private counterparts. However, it does suggest an opportunity for improvement. An important question for future work is to what degree this over-coverage is necessary? In particular, is there an inherent tension between privacy guarantee, and learning enough about the data set to accurately quantify the uncertainty?
In many estimation tasks defining a non-parametric confidence interval that gives close to nominal coverage rates is difficult. Without information regarding the underlying distribution, the confidence intervals need to be wide enough to ensure valid coverage rates for any possible distribution. This can result in the non-parametric confidence intervals having higher than expected empirical coverage when the data is drawn from a nice distribution, e.g a log-normal distribution. This effect is one possible explanation for the fact that the CDP confidence intervals have empirical coverage higher than 1 − α in Figure 8. In fact, we see evidence of this in the analysis of ExpMech. The error of the exponential mechanism is data dependent (and hence distribution dependent), but in our coverage analysis we are forced to use the worst case error of the exponential mechanism over all datasets. Figure 9: Bias of the algorithms (average of the difference between the value obtained and the true median) for lognormal data centered at 1 as the shape parameter of the distribution (the variance of the normal distribution that is exponentiated) increases from 2 to 10.
Bias. The goal of this work was to design algorithms that output valid confidence intervals for the median, not to estimate the median itself. An ad-hoc estimate of the median can be obtained from a confidence interval by taking the estimate to be the mid-point of the interval. This approach is preferable in the differential privacy context since it allocates its entire budget to the object of interest (the confidence interval) and we discussed in Section 3.1 some of the reasons why direct estimators of the median are difficult to generalise to CDP confidence interval algorithms. For all of our CDP algorithms, as well as the non-private confidence interval, this results in a biased estimator for the median, if the underlying distribution is skewed. In Figure 9, we explore the bias of the inherited median estimators. As expected, the bias increases with the skew of the data. The bias of most of the DP algorithms (except for NoisyBinSearch) is not substantially different from the bias for the non-private estimate. This implies that most of the bias can be attributed to the ad-hoc strategy of using the mid-point of the confidence interval as the point estimate for the median. As mentioned earlier, one benefit of CDFPostProcess, NoisyBinSearch and BinSearch + CDF is that they come with additional information about the distribution which could potentially be used to release a less biased estimate of the median. We leave this for future work. In this section, we illustrate how the findings from the previous sections could inform the implementation of a differentially private median release strategy in practice. We also demonstrate what level of accuracy one could reasonably expect for realistic applications. Our motivating example is the median income tables published by the U.S. Census Bureau for various subgroups of the population. Specifically, we aim to replicate a subset of statistics from Table A1 Since we want to assess the accuracy of the CDP estimates, we use income data from the 1940 Decennial Census [Ruggles et al., 2021], which enables us to compare the noisy estimates to the true values in the population. We restrict the population data to heads of households in the mountain division region (Arizona, Colorado, Idaho, Montana, Nevada, New Mexico, Utah, and Wyoming) and focus on the variables type of household (two categories), metropolitan area (two categories) and age (two categories). To mimic the illustrative application described above, we repeatedly sample from this population and treat the resulting data as the survey from which the (noisy) estimated medians will be computed. For simplicity, we draw 1% simple random samples without replacement. We acknowledge that the sampling design for the CPS ASEC is far more complex. However, as indicated in the introduction understanding the subtle effects of complex sampling designs on the privacy guarantees is currently an area of active research and is beyond the scope of this paper.

Real data application
The variable we use for our evaluations is 'INCWAGE,' which "reports each respondent's total pre-tax wage and salary income for the previous year." The amounts are displayed in "contemporary dollars," which means they are not adjusted for inflation. In the 1940's dataset, the variable is topcoded at 5, 001 dollars. We remove all N/A and missing values from the dataset, and only consider records corresponding to the head of each household. We note that we do not propose simply dropping all cases with missing values in practice as this will likely introduce bias. However, properly integrating any non-response adjustments into the DP algorithms is beyond the scope of the paper. Thus, we treat the fully observed data on household heads as our population of interest. Finally, for the purpose of error analysis we treat the empirical distribution of the entire population dataset (from which we sample 1%) as the true underlying distribution P. 6

Selecting the algorithm and (hyper)parameters
To generate the privatized confidence intervals, we use the algorithm identified as the winner in a wide range of regimes in the simulation studies: ExpMech. The hyperparameters are set to R = [0, 5001], and θ = 5. The lower and upper bounds are chosen based on the assumption that the threshold used for top coding is public knowledge and the reasonable assumption that the median income will not be less than zero. The granularity parameter θ is chosen based on Census Bureau data visualizations that report median incomes from 1967 to present, which are rounded to the nearest $100 [U.S. Census Bureau, 2020b] indicating that a granularity of $5 for 1940 median incomes is likely sufficient for data users. We split the overall privacy budget of ρ = 0.5 equally across three characteristics: type of household, metropolitan status, and age. Results based on the first simulation run are included in Table 1. Note that the CDP confidence intervals and median estimate are the result of a single run of ExpMech so this table is indicative of what we would expect in practice. The CDP point estimates for the median incomes are chosen as the midpoint of the corresponding CDP confidence intervals. Note that we could also leverage a prior assumption of the right-skewness of income data by choosing the CDP point estimator from the left half of the CDP confidence interval, rather than from its center, but we leave this type of parametric estimation to future work. We leverage the assumption that the incomes are non-negative, so we set the lower endpoint of the CDP confidence intervals at the maximum of the output of the algorithm and 0. However, we compute the point estimates before the truncation step to avoid introducing bias. The table also provides non-private and private 90% confidence intervals (the margin of error reported in the Census tables could be computed as the half-width of these intervals). The non-private median estimates are closer to the true values than the DP estimates for all sub-populations. However, for many statistics the difference between the point estimates is small relative to the width of the confidence interval indicating that the bias introduced by the ad-hoc approach of using the center of the confidence interval as the point estimator for the median is only minor. Except for the householders aged 65 and above, the relative increase in uncertainty also seems to be acceptable. The relative increase in the length of the confidence intervals ranges between 20.7% and 81.7%, that is, the uncertainty from data protection is always less than the uncertainty from sampling. The large relative increase of the confidence intervals for householders aged 65 and above can be explained if we note that the width of the CDP intervals is lower bounded by the granularity hyperparameter (which we chose to be θ = 5), which leads to a large relative uncertainty if the non-private interval has width close to 0. However, the absolute increase in uncertainty is still acceptable.

Results
In Figure 10 and Figure 11 we explore the performance of ExpMech and the non-private algorithm over 1,000 randomly sampled datasets. Figure 10, which contains boxplots showing the width of the private and non-private confidence intervals, confirms the findings based on one simulation run. While the private confidence intervals are typically wider than the non-private intervals, the increase in width is less than a multiplicative factor of two for all sub-populations except for head of households older than 65, and less than $100 in all sub-populations. The figure also reports the coverage rates for the non-private and private confidence intervals. The coverage rates are computed over 1,000 simulation runs and 20 trails of the CDP algorithm within each simulation run. While the non-private coverage rates are close to the nominal 90% coverage, the CDP confidence intervals overcover substantially. The coverage rates vary between 95.9% and 100%. The results are in line with the simulation studies. Figure 11 contains boxplots showing the variability of the private and non-private point estimates (where ExpMechPoint is the exponential mechanism point estimator). The whiskers of the boxplots indicate the 5th and 95th quantile of the empirical estimates to ensure consistency with the 90% confidence intervals reported in Figure 10. The true medians from the population are indicated by cyan lines for each of the sub-populations.
We find that from an inferential perspective the difference between the private and non-private estimates is small. For most of the estimates, the range of the boxplots overlap to a large extent and similar to Table 1 the bias is small relative to the variability in the estimates. The only estimate for which we find noticeable bias is the private estimate for non-family households. The bias arises because of the large fraction of zeros among this sub-population in the original data. Since 51% of the records in the original data report an income that is essentially zero (except for the small amount of noise that we introduce to make our data approximately continuous), the sample median will also be close to zero in many simulation runs. However, the CDP point estimate is based on the midpoint of the CDP confidence interval and the upper limit of this confidence interval will almost always be larger than the 51st quantile in the population, that is, it will almost always be larger than zero. Thus, the estimated CDP median will be biased. This highlights that the midpoint strategy can run into problems if the data is highly concentrated in certain areas of the distribution. However, such a scenario does not necessarily introduce substantial bias. This can be seen for the head of households that are 65 years or older. For this sub-population, the percentage of householders with zero reported income is about 81%. As a result, the upper limit of the CDP confidence interval is still close to zero in most simulation runs and the point estimate remains almost unbiased.

Conclusion
Measuring the uncertainty in differentially private estimates is a challenging task, especially if the input data is a sample from a larger population. In this case, both sources of randomness-the sampling error as well as the error from the CDP algorithm-need to be taken into account. If the mechanism is data dependent, the two sources are no longer independent making it difficult to quantify the uncertainty in the final output.
In this paper we addressed this challenge for the median, evaluating several strategies to obtain differentially private confidence intervals for this commonly used statistic. All the algorithms proposed produced valid non-parametric confidence intervals. We also demonstrated that directly accounting for both sources of uncertainty simultaneously allowed us to give tighter confidence bounds than relying on naive approaches that account for the two components sequentially. Our simulation results showed that an algorithm we called ExpMech produced reliable and consistent confidence intervals which were less than twice the width of the non-private confidence intervals in a wide variety of parameter regimes. A pair of algorithms called CDFPostProcess and BinSearch + CDF provide confidence intervals that are almost as tight, or slightly tighter, than ExpMech in a variety of regimes. These algorithms are practically appealing since they release a wealth of additional information about the distribution P without consuming additional privacy budget.
The private confidence intervals in the application based on the 1940 Decennial Census were not substantially wider than the confidence intervals released by the non-private algorithm, illustrating that the extra uncertainty due to data protection can be small in practice. It should be noted that the total privacy budget always needs to be divided among all the characteristics of interest (type of household, metropolitan status and age in our application) so the accuracy will necessarily decrease if more statistics are to be released under the same privacy budget.
We also found that the bias introduced by the ad-hoc strategy of using the midpoint of the confidence interval as an estimate for the median was limited for most estimates in our real data application. One strategy to further reduce this bias would be to use available information regarding the skewness of the data to come up with a better point estimate for the median. The fact that most of the algorithms provide additional information regarding the CDF of the data could be helpful for this endeavour as the information could be exploited in a post-processing step to model the distribution of the data. We leave this for future research.
We saw in our experiments on both simulated and real data that the empirical coverage rate of our private confidence intervals was often (sometimes substantially) higher than the nominal coverage rate. An interesting open question is whether this is inherent for non-parametric CDP confidence intervals for the median. Further, if this is unavoidable, then what distributional assumptions are required to narrow the gap between the empirical and nominal coverage rates?
Finally, perhaps the strongest limitation of our paper is the reliance on the assumption that the sample is drawn using simple random sampling with replacement. Such a sampling design will never be used in the survey context in practice. Thus, the important next step will be to extend the methodology to allow for more complex designs. However, this raises many challenging problems. First, all algorithms assume that the sample is iid, which is typically not true in practice and it is not obvious which adjustments would be necessary to account for this. Second, many sampling designs are informative, meaning that the sampling design is data dependent, which has consequences on the privacy guarantees that are difficult to quantify. Third, the sampling weights that would be used to compute the median (or any other quantile) would influence the sensitivity of the statistic. Fourth, further adjustments such as calibration or dealing with non-response have additional impacts on the privacy guarantees. Fifth, computing confidence intervals for the median is challenging for many sampling designs even for the non-private case. Addressing all these aspects is well beyond the scope of this paper. However, each of these aspects would be an interesting and important area of future research with impact well beyond the median application considered in this paper.

B Details: Confidence intervals based on exponential mechanism, ExpMech
Theorem B.1 (Exponential Mechanism McSherry and Talwar [2007]). Given a space of datasets X n and an arbitrary range, R, let u : X n × R → R be a utility function that maps dataset/output pairs to utility scores. For a fixed dataset d ∈ X n and privacy parameter ε ∈ R + , the exponential mechanism outputs x ∈ [r , r u ] with probability proportional to exp εu(d,x) 2∆u , where The exponential mechanism is ε-DP.
Let the range of possible outputs be R = [r , r u ]. The standard exponential mechanism to estimate the value d (k) , k ∈ [1, 2, . . . , n] (described in [Smith, 2011] for k = n/2) uses the following utility function.
This utility function captures how far x is in rank from d (k) . However, this standard mechanism estimates d (k) poorly if the datapoints in d are highly concentrated around this value. 7 Below is a variant on the standard exponential mechanism designed to perform well even in this situation.
The θ -widened utility function can be implemented in different ways; Algorithm 2 offers one method of doing so. To sample efficiently from the distribution defined by the utility function, we implement a two-step strategy as shown in prior work [Alabi et al., 2020, Cormode]: First, we sample an interval, using the fact that sampling from the exponential mechanism is equivalent to choosing the value with maximum utility score after i.i.d. Gumbel-distributed noise has been added to the utility scores [Abernethy et al., 2016]. Second, we sample an output uniformly at random from that interval.
Proof. Follows directly from Theorem B.1.
Definition B.2 ((t, θ , β )-good). Let A be a randomized mechanism that outputs a real-valued variable m. For a fixed dataset d ∈ X , θ ∈ R + , and β ∈ (0, 1), m is (t, θ , β )-good with respect to target rank k ∈ [1, 2, . . . , n] if there exists a datapoint a ∈ d such that with probability at least 1 − β , where the probability is over the randomness of A.
Proof. We will upper bound the probability density of outputs that are not (t, θ )-good with respect to target rank k, ie. outputs m for which there does not exist a datapoint a ∈ d such that |m − a| ≤ θ and |rank d (a) − k| ≤ t.
To do so, recall that the θ -widened exponential mechanism assigns utility scores to dataset/output pairs according to Definition B.1. For a given t, let us define good outputs as those having a utility score ≥ −t, which are assigned unnormalized probability density of at least 1, and bad outputs as those having a utility score < −t, which are assigned unnormalized probability density of at most exp(−tε). By definition of the θ -widened utility function, the good Algorithm 2: ExpMechPoint: θ -Widened Exponential Mechanism for Quantile Estimation Data: d = (d 1 , . . . , d n ) ∈ R n Privacy params: ε ∈ R + Hyperparams: k ∈ [n], R = [r , r u ] ⊂ R, θ ∈ R + Clip d to the range [r , r u ], setting values less than r or greater than r u to r and r u respectively. n = |d| Sort d in increasing order outputs must span an interval of size at least 2θ and the bad outputs span an interval of size at most |R| − 2θ . Therefore, we have that Setting this probability to be within β , we can solve for t as The resulting bound is tight by virtue of the worst-case example.
Next, we will consider two ways in which we can use ExpMechPoint to create a confidence interval for the median. The first (ExpMechUnion) consists of taking a union bound over the probability that the non-private interval fails to capture the true median, and the probability that the private interval fails to capture the non-private interval. The second (ExpMech) is a more nuanced approach that accounts for the noise due to sampling and noise due to privacy together.

B.1 Union bound confidence interval
The following pseudocode describes ExpMechUnion or ExpMech (depending on the boolean hyperparameter Union). The sub-algorithm ComputeExpMechTargets will be described in the next subsection, as it is only called when Union = 0.
First, we show that both ExpMechUnion and ExpMech are ε-DP.
Proof. The computations of t, β 1 , N α ε,L , and N α ε,U do not depend on the dataset d. Therefore, when analyzing the privacy loss, we simply need to consider the two calls the algorithm makes to ExpMechPoint, each with privacy parameter ε/2. By Lemma B.1, each of these algorithms is ε/2-DP, so by composition, ExpMech is ε-DP.
Then, we show that ExpMechUnion produces a valid confidence interval.
Proof. First, letting ci and ci , for any β 1 ∈ (0, 1) we have by Lemma 2.3 that Pr d med(P) < ci Then, for N α ε,L = N β 1 to rank N α ε,L . By Definition B.2, this means that Putting these together, we consider the lower endpoint of the interval [ ci α L , ci α U ]. We can upper bound the failure probability as follows.
where the second inequality follows from (3) and (4), and the final equality follows from the definition of β 1 in Algorithm 3. A similar inequality holds for the upper endpoint of the interval, so a union bound gives the desired result.

B.2 Tighter Confidence Interval
Next, we consider the more nuanced approach. Let P ∈ ∆ C (R) be a population distribution function, where med = med(P) is the population median. For a dataset d = (d 1 , . . . , d n ) where d i is sampled i.i.d. from distribution P, let rank d (a) denote the rank of real value a within dataset d. Let A ε k (d) be the output of the θ -widened exponential mechanism on dataset d that estimates the value at rank k. For a given k L , k U , we would like to control the probability that the interval [A(d, k L ) − θ , A(d, k U ) + θ ] fails to contain the true median med. In particular, for α ∈ (0, 1), we would like to find the target ranks k L and k U closest to n/2 such that In Algorithm 4 (ComputeExpMechTargets), we find these target ranks by first computing the probabilities above for all possible k L and k U 's, and then by numerically searching for the target ranks closest to n/2 such that the probabilities above are both within α/2. 8 The following lemma characterizes these probabilities.
Lemma B.5. Let P ∈ ∆ C (R) be a population distribution function, where med = med(P) is the population median. For a dataset d = (d 1 , . . . , d n ) where d i is sampled iid. from distribution P, let rank d (a) denote the rank of real value a within dataset d. Let k L , k U ∈ [1, 2, . . . , n] be target ranks, and let A ε k L (d) and A ε k U (d) be θ -widened exponential mechanisms on dataset d that estimate the value at rank k L and k U , respectively. Let C Bin and C Bin be the CDF and PDF of the binomial random variable Bin(n, 1/2). Then, Proof. For simplicity, we consider just the first statement pertaining to the lower endpoint of the interval. We can split up the probability into two cases: first, when rank d (med) < k L , and second, when rank d (med) ≥ k L . Pr For the first case, where rank d (med) < k L , we simply upper bound the first probability in the summation by 1 and note that the random variable 1 rank d (med)=m follows a binomial distribution.
In the second case, we first observe that the probability of any real value a being greater than med is monotonically increasing in a, which gives where the last inequality follows from Lemma B.2. Therefore, we have that Validity of the ExpMech confidence interval then follows directly from the selection of the target ranks.
8 This search can be implemented more efficiently by noting that k L is greater than or equal to N β 1 L − t as defined in Algorithm 3, and similarly k U is less than or equal to N β 1 U + t . Lemma B.6. Let dataset d be drawn i.i.d. from a distribution P ∈ ∆ C (R) with population median med(P). For (hyper)parameters ε > 0, R ∈ R, θ > 0, and α ∈ (0, 1), let [ ci α L , ci α U ] be the output of ExpMech(d, ε, (α, Union = 0, R, θ , ·)). If med(P) ∈ R, then with probability at least 1 − α, where the probability is over the randomness in both the dataset d and the mechanism ExpMech.
Proof. ComputeExpMechTargets (Algorithm 4, relying on Lemma B.5) returns target ranks N α ε,L and N α ε,U such that The result follows from a union bound.
Pseudo-code for the first step, which we will call GetNoisyCounts is given in Algorithm 6. Let us focus on finding the lower limit of the confidence interval. Given a target quantile q target , this algorithm iterates reduces the search domain by querying the rank of the mid-point x t of the range. If it is confident that the mid-point is to left of the target quantile then it cuts the domain in half and only keep the right half (similarly if it is confident that the mid-point is to the right, it keep the left half of the range). It continues this process until the entire privacy budget is consumed.
At each iteration we use a portion of the privacy budget ρ to release the noisy rank of the query point x t . The total privacy budget consumed by the algorithm is the sum of the privacy budget consumed by each step (Lemma 2.2). One Algorithm 6: GetNoisyCounts: ρ-CDP Algorithm Data: d = (d 1 , . . . , d n ) ∈ R n Privacy params: ρ ∈ R + Hyperparams: n ∈ N, β 2 ∈ (0, 1), q target , q L , q U ∈ (0, 1), prev-queries, R = [r , r u ] ⊂ R, ρ step , β step // prev-queries = {(x, r x , σ x )} is a collection of noisy measurements where x ∈ [r l , r u ] and r x = rank d (x) + N (0, σ 2 x ) lower = r , upper = r u // The initial search space is the entire range. ρ init = ρ step /10, β init = β step /10 // Budget for initial measurement at every query point; can be arbitrarily small. t = 0 // Counter for number of query points. ρ used = 0 // Counter for used privacy budget. while ρ used + ρ init ≤ ρ do x t = (lower + upper)/2 // Query point est-good-enough = False if there exists r x t and σ x t such that (x t , r x t , σ x t ) ∈prev-queries then avg-noisy-count t = r x t , avg-var t = σ x t est-good-enough = True numMeasurements = 0, ρ t = 0, β t = 0 while est-good-enough = False and ρ t + ρ init ≤ ρ step and ρ used n or avg-noisy-count t + K < q U · n) and (avg-noisy-count t − K > q L · n or avg-noisy-count t + K < q L · n) then est-good-enough = True if avg-noisy-count t < n · q target then lower= x t else upper = x t ρ used = ρ used + ρ t t = t + 1 return (x 1 , avg-noisy-count 1 , avg-var 1 ), (x 2 , avg-noisy-count 2 , avg-var 2 ), · · · Algorithm 7: PostProcessUnion Input: (x 1 , ns 1 , var 1 ), (x 2 , ns 2 , var 2 ), · · · , (x T , ns T , var T ), n ∈ N, q L , q U , β 2 ∈ (0, 1) l, u] option for allocating the privacy budget is to decide in advance the number of iterations and divide the privacy budget by the number of iterations to obtain a per step privacy budget. However, we can actually improve on this approach by noticing that if |rank d (x t ) − q target · n| is large then we can tolerate a lot of noise in our estimate of rank d (x t ) and still determine with high confidence whether rank d (x t ) > q target · n or rank d (x t ) < q target · n. Thus, we may be able to only allocate a very small amount of privacy budget to some steps. For a given query point x t we do not know a priori how large |rank d (x t ) − q target · n| is, and hence how much noise the query can handle. Thus, we start by adding a large amount of noise (using only a small amount of the privacy budget ρ init ) to rank d (x t ). If this noisy estimate is far enough from q target that we can confidently determine which direction to continue with the binary search, then we move and this step has only consumed ρ init privacy budget. Otherwise, we take another noisy measurement of rank d (x t ) and average the two together. This produces a less noisy estimate, and consumes 2ρ init . We continue in this way until either the variance of the estimate is low enough that we can confidently move, or this step has consumed the maximum amount of privacy budget per step ρ step , and we move in the more likely direction. While we search for the left and right hand limit of the confidence interval separately, in many settings the early query points of the binary search will be same for both. Thus, we can improve accuracy by not repeating these noisy queries. This is why we pass prev-queries into GetNoisyCounts.
The next step is processing the noisy counts to obtain a valid confidence interval. Pseudo-code is given in Algorithm 7 and a validity proof is given in Lemma C.2. This step does not consume additional privacy budget since it is simply post-processing on top of the ρ-CDP output of GetNoisyCounts.
Proof. The lemma follows immediately from Lemma 2.2 in the main text. The privacy budget ρ is divided between the two calls to GetNoisyCounts, which each use privacy budget ρ/2. By [Bun and Steinke, 2016, Proposition 1.6], each new noisy measurement noisy-count numMeasurements is ρ init -CDP, so ρ t and thus ρ used accurately capture the privacy budget consumed per step, and in total at any point during the algorithms run.
Proof. While GetNoisyCounts is designed to ensure the final output is close to the right quantile, the validity of the confidence interval really comes from the post-processing function PostProcessUnion. We have that where the subscript A denotes that the probability is over the randomness of the mechanism NoisyBinSearch. Now, let noisy-counts = (x 1 , ns 1 , var 1 ), (x 2 , ns 2 , var 2 ), · · · , (x T , ns T , var T ) be the concatenated outputs of the two runs of GetNoisyCounts in NoisyBinSearch; these are inputs to PostProcessUnion. We have the guarantee that for all t ∈ T , ns t = rank d (x t ) + N (0, var t ), and therefore with probability 1 − β 2 /2, for all implies that there exists x t such that x t < med(P) but ns t ≥ N β 1 L + R t . But this would imply that |ns t − rank d (x t )| ≥ R t . Therefore, Pr(med(P) < ci α L ) ≤ β 1 /2 + β 2 /2(1 − β 1 /2) = α/2. Similarly we can argue that Pr(med(P) > ci α U ) ≤ β 1 /2 + β 2 /2(1 − β 1 /2) = α/2, so we are done. levels of granularity ought to sum to the same values. Honaker [2015] proposes an optimal method to leverage this information, which we use here. Following Honaker's notation, label each node in the tree in binary, so the root node is 1, the left child of the root is node 10 and its right child is node 11, and so on as in Figure 12. Let the count at node i be t i , and let the leaf nodes of the tree in Figure 12 represent a histogram with granularity θ , so that t 100 is the number of data-points that lie between the left of the histogram's range, r , to r + θ , t 101 is the number of data-points between r + θ and r + 2θ , and so on. Let the histogram at any level have granularity twice that of its children. Note that if the counts were perfectly accurate, it then follows that the parent node's value should be equal to the sum of its children nodes' value, e.g. t 10 = t 100 + t 101 in Figure 12. Similarly, a child nodes' value in a perfectly accurate tree should be equal to its parent node, minus the adjacent child node's count, e.g. t 100 = t 10 − t 101 . Let iΛ1 be the neighboring child of node i (e.g. 11Φ1 = t 10 ), and let iΦ1 be the parent of node i. Honaker leverages these relationships between the counts at each of the nodes to generate an optimal tree in a recursive process. First, the counts at each node using the children node are incorporated into a weighted estimate of each of the counts, where the weight of the child counts at node i is denoted w − i and the optimized count "from below" is denoted t − i . These are then recursively combined with the counts "from above" (i.e. the count using the parent and the adjacent child) to generate t + i , which weights the count from above with w + i . Finally, the two are combined to get a fully efficient estimate of each of the nodes, with weight w and efficient count t * i . 9 In our setting, we only consider trees with an equal amount of noise on each of the nodes in the tree. This results in several simplifications of the equations in Honaker 2015.
Lemma D.1. If each noisy count in the tree has noise with variance s added to it, then the weight vectors w − and w + only need to be calculated once per level of the tree, the weight at node i for the summation from below is , the weight at node i for the summation from above is , and the optimal weight w is equal to w + .
Note that if each node has noise with variance s added to it, then ∀i, the variance at node i, . Then, the weight at node i, w i may be recursively defined as where the first line comes from Honaker's definition in his Equation 10. Similarly, from Honaker's Equation 11 it follows that when the noise at each node has the same variance, 9 See Honaker equations 10, 11, and 13 for the full statement of the values of these weights and counts.
which is equivalent to then the expression for the optimal weights in this setting. Once there is a tree of fully optimized counts, one can read off the CDF at an arbitrary point by traversing the tree in a root-to-leaf path, summing as few of the values together as possible to get the desired estimate, as shown in Algorithm 9.

D.2 Confidence Intervals for Quantiles Estimated from Tree-Based CDF
In order to generate a confidence interval for the desired quantiles, we would need to understand what the uncertainty of each of the counts in our estimated CDF was. Since each of these counts is a combination of all of the different counts in the tree, weighted in a way that is recursively defined, this is not trivial to do in a closed-form manner. We can compute the effect of each node on any other node by generating a tree with every node valued at 0 except the node we are interested in, then running the recursive weighting algorithm on that tree, as shown in Alg. 10.

Algorithm 10: ComputeNodeEffect
Input: σ 2 ∈ R ≥0 , i * ∈ N, m ∈ N Construct a binary tree T ∈ BinTree(N, m), where for 0 ≤ i < 2 m , Let T ∈ BinTree(N, m) be the output of the CDF post-processing algorithm from Honaker [2015] on differentially private tree T , where each noisy count in T has variance σ 2 . return T We now have a method to understand how much each node effects any other node. If we run this on every single node of the tree, we can then combine them to generate a tree for every node on the tree that describes how much its optimized count is affected by any other node. 10 When summing the counts to generate the CDF, we can then keep track of the total weight of each node in the final count, and from here generate the variance of the count.
interval generation process is summarized in Alg. 13.
Proof. Note that the only step in CDFPostProcess(Union) that touches the dataset d is the call to DPTree, which creates a tree of m differentially private histograms. Each histogram is ρ/m-CDP, and by composition ([Bun and Steinke, 2016, Proposition 1.6]), DPTree is a ρ-CDP algorithm. The rest of the computations in CDFPostProcess(Union) apply post-processing to the output of DPTree, so they do not affect the privacy guarantee. Lemma D.3. For any dataset d iid ∼ P, where P ∈ ∆ C (R), and any hyperparameters θ , R = [r , r u ], γ ∈ (0, 1), failure rate α and privacy parameter ρ, let CDFPostProcess(d, ρ, (α, Union = 0, R, θ , γ)) return an interval [ ci α L (d), ci α U (d)]. If med(P) ∈ R, then with probability at least 1 − α, where the probability is taken over the randomness of both the dataset d and the mechansim CDFPostProcess.
Proof. Let us consider the upper endpoint of the interval. First, given a set of DP measurementsC(x) =Ĉ(x) + Proof. BinSearch + CDF is a composition of two algorithms -NoisyBinSearch which by Lemma C.1 is γρ-CDP, and CDFPostProcess which by Lemma D.2 is (1 − γ)ρ − CDP. By Lemma 2.2, this means BinSearch + CDF satisfies ρ-CDP.
The coverage analysis of BinSearch + CDF follows immediately from Lemma D.3 and Lemma C.2, and a union bound.
where the probability is over the randomness of both the dataset d and the mechanism BSCDF.

F Details: Other Algorithms Explored
In this section we give a brief overview of additional CDP confidence intervals and CDP median estimators that we explored. These algorithms were not included in the main body of this paper since they are outperformed by other algorithms in every parameter regime we studied. The additional CDP confidence interval algorithms were: • CDF+BS CI computes a CDP estimate to the empirical CDF in the same way as CDFPostProcess. However, instead of using the post-processing algorithm described in Algorithm 12, it performs binary search using the noisy CDF measurements.
• BinSearch is the same as NoisyBinSearch except it uses the same privacy budget at every iteration. We expect this algorithm to perform strictly worse than NoisyBinSearch which uses its budget more carefully. Figure 13 shows the performance of CDF+BS CI and BinSearch, as well as the naive estimators ExpMechUnion and CDFPostProcessUnion and the four CDP estimators we presented in the main body. We can see that for all values of ρ, at least one of the four main CDP estimators outperforms each of the other algorithms.   Figure 13: Performance of various CDP confidence intervals for the median as we vary ρ. Performance is measured in terms of the relative width with α = 0.05. Box plots are computed using 100 random datasets of 1000 data points drawn i.i.d. from Lognormal(ln(1.5), 1). Each CDP algorithm is run 5 times on each dataset.
We also explored several CDP median estimators. For point estimators that directly correspond to analogues of our CDP confidence intervals, we use the same name to denote both. However, note that these point estimators for the median are different to the mid-point of the confidence interval estimators that we used in the main body. The algorithms presented in Figure 14 are all directly estimating the median, and do not additionally release a confidence interval for the median.
• ExpMech is the point estimator version of our confidence interval algorithm ExpMech. It uses the exponential mechanism with target quantile n/2 to estimate the median.
• SmoothSens releases the median using the smooth sensitivity framework Nissim et al. [2007], Bun and Steinke [2019]. This algorithm Gaussian noise to the empirical median where the standard deviation of the noise is data dependent, and carefully calibrated to ensure differential privacy.
• BinSearch is the point estimator version of BinSearch described above. It uses binary search with target quantile n/2.
• NoisyStartBinarySearch was a preliminary version of altering the privacy budget through the iterations of the algorithm, with little budget initially then increasing through the search process.
• NoisyBinSearchis the point estimator version of our confidence interval algorithm NoisyBinSearch. It uses noisy binary search to search for the quantile n/2.
• FancyBinarySearch is similar to BinSearch. However, instead of halving the range at each iterate, it makes more conservative steps when it is not confident whether the median is to the left or right.
• CDFPostProcess is the point estimator version of our confidence interval algorithm CDFPostProcess. It computes a CDP estimate the CDF in the same way, then computes the median based on the CDP CDF. An important note is that since both this algorithm and CDFPostProcess are post-processing on the CDP CDF estimate, they can be performed at the same time without additional privacy budget.
• GradDescent uses CDP gradient descent to solve the optimisation problem arg min ∑ n i=1 |m − d i |. We use the private stochastic gradient descent technique proposed by Bassily et al. [2014]. Figure 14 shows the performance of each of our point estimators on log-normal data. We can see that SmoothSens, the only unbiased estimator, has among the highest variability in all regimes, and particularly poor performance for small ρ. It has comparable performance to the other algorithms for large ρ, but extending this point estimator to a confidence interval algorithm remains an open problem. All the variants of binary search perform similarly as point estimators for the median. Even as a point estimator, ExpMech slightly outperforms the other algorithms, except GradDescent. Extending GradDescent to a confidence interval remains an open problem. Figure 14: Performance of various CDP point estimators for the median as we vary ρ. Box plots are computed using 100 random datasets of 1000 data points drawn i.i.d. from Lognormal(ln(1.5), 1). Each CDP algorithm is run 5 times on each dataset.

G Other Regimes
The relative ordering of algorithms can depend on the scale of the data (σ d ) relative to the range. In the figures below, we display the relative widths of the algorithms on data sampled from a Lognormal(ln(1.5, σ 2 d ) distribution, where σ d = 5.0, as we vary the size of the dataset n and the privacy loss parameter ρ. Note that although we have drastically increased the scale of the data, the range is left the same as in Figure 6: R = [−5, 15]. From these plots, we can see that when n, ρ, and σ d are large, CDFPostProcess performs slightly better than ExpMech. In Figure 6, we saw that when σ d is small, ExpMech remains the best performing algorithm in both the large n and large ρ regimes. Hence we conjecture that σ d needs to be large, and either ρ or n need to be large for CDFPostProcess begins outperforming ExpMech. This conjecture is supported by Figure 15, where we see CDFPostProcess only beginning to outperform ExpMech when either n and ρ are large.