Differentially Private Bootstrap: New Privacy Analysis and Inference Strategies

Differential private (DP) mechanisms protect individual-level information by introducing randomness into the statistical analysis procedure. While there are now many DP tools for various statistical problems, there is still a lack of general techniques to understand the sampling distribution of a DP estimator, which is crucial for uncertainty quantification in statistical inference. We analyze a DP bootstrap procedure that releases multiple private bootstrap estimates to infer the sampling distribution and construct confidence intervals. Our privacy analysis includes new results on the privacy cost of a single DP bootstrap estimate applicable to incorporate arbitrary DP mechanisms and identifies some misuses of the bootstrap in the existing literature. We show that the release of $B$ DP bootstrap estimates from mechanisms satisfying $(\mu/\sqrt{(2-2/\mathrm{e})B})$-Gaussian DP asymptotically satisfies $\mu$-Gaussian DP as $B$ goes to infinity. We also develop a statistical procedure based on the DP bootstrap estimates to correctly infer the sampling distribution using techniques related to the deconvolution of probability measures, an approach which is novel in analyzing DP procedures. From our density estimate, we construct confidence intervals and compare them to existing methods through simulations and real-world experiments using the 2016 Canada Census Public Use Microdata. The coverage of our private confidence intervals achieves the nominal confidence level, while other methods fail to meet this guarantee.


Introduction
In the big data era, individual privacy protection becomes more critical than ever because personal information is collected and used in many different ways; while the intention of the data collection is usually to improve the user experience or, more generally, to benefit society, there have also been rising concerns about malicious applications of these data. To protect individuals against arbitrary attacks on their data, Dwork et al. (2006) proposed differential privacy (DP) which has become the state-of-the-art framework in privacy protection. DP protects the sensitive information in the input of an algorithm by introducing additional randomness into the procedure. For statistical analysis, a great deal of prior work focused on producing private point estimates of a parameter, e.g., the sample mean, sample median, and the maximum of the data points. In contrast, while some prior work aims to quantify the uncertainty of a DP procedure, their techniques and results are usually restricted to specific settings, and there is still a lack of general-purpose methods (see related work).
One of the most widely used methods to approximate a sampling distribution is the bootstrap method (Efron, 1979). The bootstrap can be used to quantify uncertainty in many statistical ways, such as by producing a non-parametric confidence interval (CI).
Although the bootstrap has been studied very well in Statistics, it is still an open question of how to build and analyze a DP bootstrap for private statistical inference. Brawner and Honaker (2018) were the first to propose and analyze a DP bootstrap procedure and used it to produce a CI based on the private bootstrap outputs. However, a key step in their privacy proof is incorrect 1 , and we show in Section 3 that their stated privacy guarantee does not hold. Furthermore, the performance of their CI is not satisfactory in terms of both coverage and width.  developed the state-of-the-art analysis of resampling for DP, both with and without replacement. As a particular case, their results can be used to analyze the privacy guarantees of mechanisms with their input being bootstrap samples.
However,  did not consider the cumulative privacy cost of multiple samples of the resampling methods, which restricts the usage of their results on the bootstrap. This is non-trivial, mainly because their results are presented in (ε, δ)-DP which, unlike f -DP, does not have tight composition techniques.  also did not consider how to perform statistical inference on the DP bootstrap samples, which also limits the usability of their results.
Our contributions In this paper, we obtain a tight privacy analysis of a DP bootstrap method and develop inference strategies on the sampling distribution. Specifically, we derive the privacy guarantee of the DP bootstrap, which generalizes the result by  from (ε, δ)-DP to f -DP, and an aspect of our proof strategy applies to any mixture of DP mechanisms where the bootstrap is a special case. Our result also identifies misuses of resampling with replacement in the literature. We derive the asymptotic cumulative privacy cost of multiple DP bootstrap estimates via the central limit theorem (CLT) for f -DP (Dong et al., 2022). We demonstrate how to use deconvolution on the DP bootstrap results to obtain a private estimate of the non-private sampling distribution and develop CIs. We provide simulations to show the advantage of our method in terms of the coverage of the resulting CIs. We also conduct real-world experiments on the 2016 Canada Census Public Use Microdata, which reveals the income difference between two provinces and a positive correlation between income and shelter cost under the DP guarantee.
Related work For DP confidence intervals, Du et al. (2020) provided a discussion and comparison between different methods (including NoisyVar, which we discuss in our experimental sections) providing DP CIs for the mean of normally distributed data. Sheffet (2017) presented a DP estimator for the t-values in ordinary least squares and derived the CI based on the t-values. Karwa and Vadhan (2018) gave DP CIs for the population mean of a normal distribution along with finite sample coverage and lower bounds on the size of the DP CI. D' Orazio et al. (2015) presented algorithms for releasing DP estimates of causal effects, particularly the difference between means estimators along with their standard errors and CIs. Brawner and Honaker (2018) gave the first attempt to use the DP bootstrap to obtain a private CI. Covington et al. (2021) used the bag of little bootstraps (BLB) and the CoinPress algorithm (Biswas et al., 2020) to privately estimate the parameters' sampling distribution through its mean and covariance. Drechsler et al. (2021) proposed several strategies to compute non-parametric DP CIs for the median.
It is common to design DP mechanisms using either 1) resampling: replace the input by a dataset resampled from the original dataset (also known as subsampling if each entry of the dataset is resampled for at most once), and/or 2) composition: run multiple mechanisms and release all their outputs. For (ε, δ)-DP, the classical composition results are given by Dwork et al. (2010); Ullman (2017) derived a tight bound for subsampling, and  gave the state-of-the-art results for resampling with and without replacement. For privacy loss distributions (PLD), Koskela et al. (2020) derived the PLDs for resampling with and without replacement for the Gaussian mechanism. For Rényi-DP, the state-of-the-art in composition results are given by Mironov (2017), andWang et al. (2019) derived a tight bound for subsampling. For concentrated DP (CDP) and zero concentrated DP (zCDP), there are tight composition results (Dwork and Rothblum, 2016;Bun and Steinke, 2016), but no results for resampling. Bun et al. (2018) proposed truncated CDP (tCDP) as an extension of zCDP, which has results for both composition and resampling without replacement. For f -DP, Dong et al. (2022) gave results for both composition and subsampling, along with a central limit theorem under composition. Zhu et al. (2022) derived 'dominating pairs' for composition and subsampling (subset subsampling and Poisson subsampling), which could be used to calculate any of the above DP guarantees.
Organization The remainder of this paper is organized as follows. In Section 2, we review the definition of f -DP and results used in our DP bootstrap analysis. In Section 3, we provide our privacy guarantee for the DP bootstrap along with a central limit theorem result for the cumulative privacy cost of many DP bootstrap outputs. In Section 4, we propose multiple methods for performing statistical inference on the DP bootstrap, such as using estimates of the population parameter and its variance to construct standard confidence intervals, and using deconvolution to estimate the sampling distribution and construct percentile intervals with privacy guarantees. In Section 5, we demonstrate through a simulation study that the CIs derived from the DP bootstrap have better coverage than an existing method, NoisyVar, mentioned in related work. In Section 6, we apply the DP bootstrap to answer two questions on a Canadian census dataset under privacy requirements: 1) whether there is difference between the mean income in Ontario and Quebec; 2) whether the income is positively correlated with the shelter cost in each of the two provinces. In Section 7, we discuss the implications of our work and highlight some directions for future research. Proofs and technical details are deferred to the Supplement.

Background in differential privacy
In this section, we provide existing results that are closely related to our DP bootstrap method. First, we discuss the definition of f -DP. Then we introduce the Gaussian mechanism to guarantee DP, and prior results on subsampling, composition, and group privacy.
For a set X , a dataset D ∈ X n with cardinality |D| = n is a finite collection of elements from X . We define D 1 k D 2 if |D 1 | = |D 2 | and they differ in k entries. We call D 1 and D 2 as neighboring datasets if D 1 1 D 2 , and we also write it as D 1 D 2 . A mechanism M : X n → O is a randomized algorithm taking a dataset as input and outputting a value in some space O. Differential privacy measures how much the mechanism changes when it is run on two neighboring datasets.
In this paper, we use the f -DP framework which measures the difference between two distributions by hypothesis testing. With an observation M(D in ) = o, we consider a hypothesis test between H 0 : D in = D and H a : D in = D . Intuitively, the harder this test is, the stronger the privacy guarantee that the mechanism M provides. Wasserman and Zhou (2010) first used this hypothesis testing framework in DP, and Dong et al. (2022) generalized it to the f -DP framework using the tradeoff function which characterizes the difficulty of this test (also known as the receiver operating characteristic (ROC) curve.) an important subclass is G µ (α) = T N (0,1),N (µ,1) (α). G µ -DP is also called µ-Gaussian DP (GDP), which is shown to be the limit of many DP procedures under composition by Dong et al. (2022) (details are in Theorem B.9). Another important subclass is f ε,δ (α) = max{0, 1 − δ − e ε α, e −ε (1 − δ − α)} since f ε,δ -DP is equivalent to (ε, δ)-DP. We include the definition of (ε, δ)-DP in the Supplement Section A for easier reference.
The Gaussian mechanism Let the 2 -sensitivity of a function m : For any m : X n → R d , the Gaussian mechanism on m adds Gaussian noises to the output of m: Dong et al., 2022).
Subsampling Let D ∈ X n be a dataset and D := Sample m,n (D) ∈ X m be chosen uniformly at random among all the subsets of D with size m ≤ n (sampling without replacement). For a mechanism M : Composition The composition property of DP quantifies the cumulative privacy cost of several DP outputs. If f = T P,Q and g = T P ,Q , their tensor product is defined as f ⊗ g := T P ×P ,Q×Q where P × P is the product measure of P and P . Let M i be f i -DP, Dong et al., 2022).
Group privacy While f -DP guarantees protection of the privacy of each individual, it can be generalized to give a privacy guarantee for groups of size k. We say a mechanism M If a mechanism is µ-GDP, then it is kµ-GDP for groups of size k (Dong et al., 2022).

Privacy analysis of bootstrap resampling
In this section, we first use techniques from Dong et al. (2022) to convert the state-ofthe-art (ε, δ)-DP guarantee for the bootstrap  to f -DP. However, the resulting formula is computationally intractable. We then give a new f -DP bound with a direct proof, which agrees with the result from , but is more transparent and computationally friendly. At the end of this section, we analyze the composition of several DP bootstrap releases.
3.1 f -DP guarantee with one bootstrap sample as input Bootstrap sampling is denoted by a randomized mapping boot n : is the probability that a given entry of D is included i times in boot n (D). We also write boot n and p i,n as boot and p i respectively when n is known from the context.
Although this representation of f M•boot is seemingly simple, it is hard to compute or visualize this tradeoff function because evaluating f M•boot (α) requires solving over n optimization problems for each α. It is also hard to derive composition results from this f M•boot (α), which is crucial for using the bootstrap for statistical inference since multiple bootstrap samples will be used.
Due to the intractability of Proposition 3.1, we prove a new f -DP result for M • boot using the mixture of tradeoff functions. The fundamental idea of our proof is to view the bootstrap as a mixture distribution and decompose the overall tradeoff function into a mixture of tradeoff functions where each one is easy to obtain.
Definition 3.2 (Mixture of tradeoff functions). For i = 1, 2, . . . , k, let f i be a tradeoff function and p i ∈ (0, 1] satisfying as follows: For each α ∈ (0, 1), we find C such that α ∈ A(C) and the corresponding Intuitively, as illustrated in Figure 2, by matching the slopes of each tradeoff function f i , we minimize the overall type II error given a fixed type I error. In Theorem 3.3, we always gives a lower bound on the privacy cost of an arbitrary mixture mechanism. Note that this general result applies to any mixture of DP mechanisms.
Theorem 3.3. For i = 1, 2, . . . , k, let f i be a tradeoff function and M i : be a mechanism satisfying f i -DP. Let M : X n → Y be a mixture mechanism which randomly selects one mechanism from k mechanisms, {M i } k i=1 , with corresponding probabilities , and the output of M will be the output of M i if M i is selected. Then shown as the curves in the figures on the bottom row. Each circle dot on the mixture curve is the average of the circle dots on f 1 , f 2 , f 3 weighted by p 1 , p 2 , p 3 . Note that mix( , nor is it a member of the GDP family.
The proof of Theorem 3.3 is based on using different rejection rules when different M i is selected to freely allocate the overall type I error to each rejection rule to minimize the overall type II error. We can combine such rejection rules as one if Y i are disjoint; therefore, Theorem 3.3 is not improvable.
For the DP bootstrap setting, since each M i maps to the same Y i , we can leverage this fact to strengthen the privacy guarantee in Theorem 3.4. We separately consider the case that a given entry of D is not included in boot(D) -in this case, we use a similar idea to the subsampling result of f -DP to improve our bound.
Theorem 3.4. Let f i be a symmetric tradeoff function for i = 1, 2, . . . , n. For a mechanism M satisfying group privacy f i -DP with group size i, M • boot satisfies f M•boot -DP where We can use this f k in Theoerem 3.4 if no tighter bound is known.
Remark 3.6. The (ε, δ)-DP result by  can be derived from our f -DP result in Theorem 3.4; see Proposition B.8 in the Supplement. While  provided specific settings attaining their bound which also prove the optimality of our result, it is unknown for general settings how to construct rejection rules achieving each pair of the type I error and type II error given by our tradeoff function. Therefore, it may be possible to further improve the privacy analysis for the DP bootstrap. Nevertheless, we show in Figure 3 that our bound is fairly tight for some extreme settings and should suffice for most purposes.
Our result in Theorem 3.4 is easier to compute when ∂f i are all known, since for a subdifferential value C, we can immediately obtain its corresponding α and β. For a given α or β, we can use the bisection method to search for C. From the following example and the composition result in Section 3.2, we see that our bound can be easily evaluated for the Gaussian mechanism.
Example: DP bootstrap using the Gaussian mechanism We analyze the DP bootstrap using the Gaussian mechanism to illustrate the results in Theorem 3.4. In Figure   3, we assume that the Gaussian mechanism M satisfies 1-GDP, then the DP bootstrap The solid curve is tight as a lower bound of the transparent curves. The dashed and dotted dashed lines are misused as lower bounds by Brawner and Honaker (2018) and Koskela et al. (2020).
by Brawner and Honaker (2018); 2) the DP bootstrap cannot be correctly analyzed using the privacy loss distribution (PLD) method as in (Koskela et al., 2020). Details of this example are in the Supplement Section C.

Composition of bootstrap Gaussian mechanisms
In this section, we derive a composition result for the DP bootstrap because we need many bootstrap samples to obtain the sampling distribution when we use the bootstrap method.
Although numerical algorithms can be used to evaluate the privacy guarantee of the composition, the calculation is often burdensome (Zheng et al., 2020). We prove an asymptotic result to understand the behavior as the number of compositions goes to infinity. Note that for simplicity, we assume the initial mechanism satisfies GDP, but it may be possible Then Note that the asymptotic privacy guarantee of the composition result above is the same as running a ( 2 − 2/e)µ B -GDP mechanism on the original dataset (not on the bootstrap sample) for B times (Dong et al., 2022). Therefore, the factor ( 2 − 2/e) = 1.12438 . . . < 1.125 is the price we pay for using the DP bootstrap (see Figure 4). While there is a small increase in the privacy cost, the bootstrap samples now contain the randomness from sampling as well as from the privacy mechanism. In the next section, we will see how we can use this fact to perform statistical inference with the DP bootstrap estimates.

Private statistical inference with DP bootstrap
In this section, we review two classic methods for statistical inference with bootstrap and modify them to obtain valid inference with DP bootstrap. First, we demonstrate how to use DP bootstrap to estimate the mean and variance of the sampling distribution and construct standard confidence intervals (CIs). However, the coverage of this type of CI is not satisfactory (Brawner and Honaker, 2018). We propose another approach using deconvolution to recover the correct sampling distribution from the DP bootstrap estimates. The deconvolved distribution satisfies the same f -DP guarantee because of the post-processing property (Dong et al., 2022, Proposition 4). At the end of this section, we provide simulation results on deconvolution with the DP bootstrap Gaussian mechanism.
Throughout the rest of this paper, to simplify the statements and make our results comparable to other private inference methods, we use the Gaussian mechanism in DP bootstrap and its corresponding asymptotic privacy guarantee, µ-GDP. We include remarks when other types of mechanisms can be used. Note that the asymptotic privacy guarantee can always be replaced by accurate composition results on our f -DP result in Theorem 3.4 if strict privacy guarantee is needed (see Supplement Section A for numerical techniques.)

Statistical inference with Efron's bootstrap
Before explaining how to conduct private statistical inference with the DP bootstrap, we briefly review two basic inference methods with the original Efron's bootstrap estimates (Efron and Tibshirani, 1994), and we will develop our private inference methods based on these classic methods. We denote the original dataset by D ∈ X n and the estimator for a population parameter θ by g(D). The sampling distribution of g(D) can be estimated from the bootstrap estimates, {g(D j )} B j=1 where D j is the jth bootstrap sample of D.
Using the idea of normal approximation to the sampling distribution, we can obtain the standard confidence interval : we first estimate the standard error of g(D) byŝ 2 , then approximate the sampling distribution of g(D) by N (θ,ŝ 2 g,B ), and finally build CIs for θ using g(D) andŝ 2 g,B .
If we approximate the sampling distribution directly using the empirical distribution of {g(D j )} B j=1 , we can build the percentile interval for θ using the percentiles of {g(D j )} B j=1 .

DP point estimates and standard confidence intervals
The DP bootstrap estimates using the Gaussian mechanism are {g( ; 2) g(D) is no longer accessible due to the privacy constraint and cannot be used as the center of the confidence interval. Therefore, we estimate g(D) usingm g,B and replaceŝ 2 g,B with (s 2 g,B ) =s 2 g,B − (∆g) 2 µ 2
For better illustration, we look at a specific problem, the population mean estimation, and our results are in Table 1. The derivation is in the Supplement Section D where we assume the support of the population distribution is [0, 1], simplifying the sensitivity calculation. Note that our asymptotic privacy analysis requires a constant µ with respect to n and B. In Table 1, the variance of the point estimate from DP bootstrap, Var(θ 4 ), is at the same rate as the variance of the other three estimates. When B ∈ O(n), the Sample mean (Gaussian mechanism): Bootstrap (non-private): Unbiased estimate of Var(θ i ) Variance of the estimate Sample variance (non-private): Var(Var (θ 4 )) ∈ O( 1 n 2 B + 1 n 3 + B n 4 ) variance of the variance estimates from DP bootstrap, Var(Var (θ 4 )), is at the same rate as the non-private bootstrap, Var(Var (θ 3 )). Furthermore, when B ∈ Θ(n), all Var(Var (θ i )), i = 1, 2, 3, 4, are at the same rate.
Although the estimatesm g,B and (s 2 g,B ) have nice convergence rates, they have severe problems: 1 Bs 2 g,B + (s 2 g,B ) could be negative, and the coverage of the CIs is lower than their nominal confidence level as shown in the Figure 7b in (Brawner and Honaker, 2018). In the following subsection, we develop a new inference method to solve this issue based on deconvolution.

Deconvolution for estimating sampling distribution
Instead of using the standard confidence interval, which uses normal approximation and has potentially negative variance issues, we propose to obtain the percentile intervals using deconvolution on DP bootstrap estimates. Note that although we analyze this method on DP bootstrap with Gaussian mechanism, it can be applied to DP bootstrap with any additive noise mechanisms, e.g., Laplace, Truncated-Uniform-Laplace (Tulap) (Awan and Slavković, 2020), canonical noise distributions (Awan and Vadhan, 2021).
Deconvolution is the major tool to solve the contaminated measurement problem, i.e., we want to estimate the distribution of and e i is the measurement error. This is exactly the relationship between ). Since DP allows for all details of the privacy mechanism to be publicly revealed (except the private dataset), the distribution of the added noises can be incorporated into our post-processing without raising any privacy concerns.
Deconvolution has been well-studied in Statistics, and it is known that the difficulty of deconvolution is determined by the smoothness of the noise distribution: the smoother, the harder (Fan, 1991). For the model Y = X + e, if e follows a normal distribution N (0, σ 2 ), which is super-smooth, the optimal rates for recovering where k is the smoothness of f X , the unknown density of X, and N is the sample size of Y (Fan, 1991). Although this slow rate indicates the difficulty of deconvolution with normal noises in general, if the scale of the normal noise is not too high, nonparametric deconvolution can still be practical (Fan, 1992): the optimal rate is O(N −k/(2k+1) ) when σ = cN −1/(2k+1) for some constant c. This is the best possible rate since it is the same as the optimal rate for estimating f X directly from the sample of X. This optimal rate can be achieved by a kernel density estimatorf N, is the Fourier transform of a specified kernel function K(·) with φ K (0) = 1, φ is the characteristic function of , h N = dN −1/(2k+1) , and d > 0 is a constant.
In our DP bootstrap setting, the sample size N for deconvolution is the number of bootstrap estimates B, andf B,h B (x) can be seen as a generalization of the point estimates used in the standard confidence interval (Hall and Lahiri, 2008): if we fix B and definẽ The deconvolution method is better since it does not require normal approximation to build confidence intervals compared to the point estimates in Section 4.2. We also circumvent the problem caused by the possibly negative point estimation of the variance of the sampling distribution. However, it is difficult to analyze the convergence rate of deconvolution in DP bootstrap with respect to B and n because the distribution of the added noises e flattens when B increases, and the distribution to be recovered is from bootstrap estimates. At the end of this section, we propose a signal-noise-ratio (SNR) measure as a rule-of-thumb for the choice of B.
We compared multiple numerical deconvolution methods through preliminary simulations with different settings on n and B. Among different deconvolution methods, we choose to use deconvolveR (Efron, 2016) since it performs the best in our settings without tuning its hyper-parameters. We briefly summarize this method as follows. For the model Y = X + e, deconvolveR assumes that Y and X are distributed discretely with the sizes of their supports |Y| = k and |X | = m. Then it models the distribution of X by f (α) = e Qα /c(α) where Q is an m × p structure matrix with values from the natural spline basis with order p, ns(X , p), and α is the unknown p-dimensional parameter vector; c(α) is the divisor necessary to make f sum to 1. The estimation of the distribution of X is obtained through the estimation of α: We estimate α by maximizing a penalized log-likelihood m(α) = l(Y ; α) − s(α) with respect to α where s(α) is the penalty term, and l(Y ; α) is the log-likelihood function of Y derived from the distribution of X and e (the distribution of X is determined by f (α), and the distribution of e is known).
After deconvolution, we construct the CI with confidence level α using the α 2 and 1 − α 2 percentiles of the recovered sampling distribution. Figure 5 is an example of deconvolution on the DP bootstrap estimates with the following settings: n = 10000, iid ∼ Unif(0, 1) for 6 different replicates. Since the 'private bootstrap' density curve is much flatter than the 'non-private bootstrap' curve, conducting statistical inference directly with the distribution of the estimates from the DP bootstrap with Gaussian mechanisms gives results inconsistent with the original bootstrap. We see that in Figure   5, deconvolution corrects this problem.
Remark 4.1. Although the deconvolution method is specifically designed for additive noise mechanisms, e.g., Gaussian mechanism, we can potentially approximate more general mechanisms as additive when the sample size n is very large, e.g., exponential mechanism (Awan et al., 2019;Reimherr and Awan, 2019).

Remark 4.2.
As it is difficult to analyze how to choose B theoretically, we provide a rule of thumb where we define the signal-noise ratio (SNR) as the ratio between the variance of the sampling distribution and the variance of the noise added for DP: 1) If we have enough privacy budget, we choose the largest B satisfying SNR ≥ 1; 2) If we only have a very limited privacy budget, e.g., SNR ≤ 1 for any B ≥ 10, we choose the smallest B such that B ≥ 2/α for 1 − α confidence intervals. Note that for the case of population mean inference we have discussed in Section 4.2, the variance of non-private bootstrap estimates , and for DP bootstrap, the added noise µ 2 n 2 ); therefore, our rule of thumb for building percentile intervals using deconvolution suggests B ∈ O(n) which is consistent with our rate analysis in building standard confidence intervals.
Remark 4.3. To the best of our knowledge, we are the first to use deconvolution to approximate the non-private sampling distribution from the noisy DP estimates. Farokhi (2020) demonstrated a related application of deconvolution to estimate the data distribution from locally differentially private data.

Simulations
In this section, we use the confidence interval construction as a showcase for our DP bootstrap algorithm. First, we compare the non-private CI from the original bootstrap to the private CI from the recovered sampling distribution by the deconvolution on our DP bootstrap estimates. Then we compare our DP bootstrap with the NoisyVar (Du et al., 2020) and discuss the inherent difference between the two methods.

Private CI compared to non-private CI
Consider the setting n = 3000, The results for 90% CIs are in Table 2. The coverage and the width of the DP bootstrap are similar to the non-private bootstrap when µ = 1, while the DP bootstrap is more conservative when µ < 1. We also examine the results with other confidence levels in Figure   6: Since the CI is built with the quantiles of the recovered distribution and the coverage is determined by whether the true parameter value is in between the two quantiles, we can directly check the recovered CDF, F * , at the true parameter value θ, e.g., 0.05 ≤ F * (θ) ≤ 0.95 is equivalent to the 90% CI covering θ. Therefore, the coverage at different confidence levels can be calculated by E[I(p lower ≤ F * (θ) ≤ p upper )], and we want it to be close to p upper − p lower which is the nominal confidence level. This is achieved if the CDF of F * (θ) is close to the line F (x) = x, ∀x ∈ [0, 1]. In Figure 6, when µ = 1, we can see that our DP bootstrap result aligns with the F (x) = x, similar to the non-private bootstrap. For µ = 0.1, the curve for the DP bootstrap is lower than F (x) = x for x < 0.5 and higher than F (x) = x for x > 0.5, which means our CI is more conservative since the actual coverage is higher than confidence level.

Comparison between DP bootstrap and NoisyVar
In the non-private setting, the parametric CI for the population mean θ can be built with the t-statistic t = (x − θ)/ s 2 X /n wherex and s 2 X are the sample mean and variance.
Similarly, to construct a DP CI, one can obtain a DP t-statistic by replacingx and s 2 X with their corresponding DP statistics, but this DP t-statistic may not follow the t-distribution because of the additional noise for privacy. Du et al. (2020) adopted the idea of parametric bootstrap to construct the confidence interval for the population mean: They used simulation with normally distributed data generated by plugging in the differentially private releases of sample mean and sample variance. We compare our DP bootstrap with their algorithm, NoisyVar. We include NoisyVar in the Supplement (Algorithm 4) for easier reference.
In Table 2, the performance of NoisyVar is not as satisfactory as in (Du et al., 2020) in terms of the coverage, which is because the distribution of our original data, Uniform distribution, is not in the parametric distribution family used by NoisyVar. Our comparison highlights the importance of non-parametric inference.
Note that the NoisyVar is limited to building DP CIs for the population mean, while our DP bootstrap can be used on any DP statistic with additive noise mechanisms. As a non-parametric method, the bootstrap CI does not assume the family of the sampling distribution or data distribution. We use deconvolution to correct the distribution of DP bootstrap estimates to the true sample distribution. Therefore, our results do not suffer from the under-coverage issue. When the privacy guarantee is stronger, our results are more conservative, and the coverage is higher than the nominal coverage level in our simulation.

Real-world experiments
In this section, we conduct experiments with the 2016 Census Public Use Microdata Files (PUMF), which provide data on the characteristics of the Canadian population (Canada, 2019). Under specific privacy-protection requirements, we analyze the mean difference between the two provinces' average market income and the covariance between market income and shelter cost.

Experiment settings
The where the input datasets are sampled with size n = 200, 000 from the original dataset with replacement. The privacy guarantee is set to be 1-GDP, and we use B = 100 for bootstrap and DP bootstrap methods.

Mean difference in MRKINC between Quebec and Ontario
We answer the question of whether the mean of MRKINC is different between Quebec and Ontario by constructing the confidence intervals of the mean of MRKINC and checking the overlap of the confidence intervals. Similar to our simulations, we use the Gaussian mechanism in DP bootstrap and compare our results with NoisyVar. The results are shown in Table 3. For Ontario, the true mean is 42865.35. For Quebec, the true mean is 37074.50.
Although the distribution of MRKINC in this real-world dataset is different from the uniform distribution in our simulation, as shown in Figure 7, the results are consistent with what we have observed in the simulations: DP bootstrap is as good as or better than bootstrap in terms of CI coverage. When the confidence level is high (99%), the width of CI from DP bootstrap is much wider than that from bootstrap or NoisyVar due to the

Covariance between MRKINC and SHELCO
We answer whether the covariance between MRKINC and SHELCO is positive for people in Quebec and Ontario by constructing the CIs of covariance between MRKINC and SHELCO and checking whether those CIs only contain positive values. For Ontario, the true covariance is 3352465. For Quebec, the true covariance is 2405268. Similar to our simulations, we use Gaussian mechanism in DP bootstrap: If X ∈ [0, 1] and Y ∈ [0, 1], the 2 -sensitivity is 1 n for Cov(X, Y ) (the same as for Var(X) and Var(Y )). We propose NoisyCov by following the idea of NoisyVar. In Table 4, we compare our method with both the non-private bootstrap and NoisyCov.
From the results in Table 4, we can see the coverage of the CIs from DP bootstrap is the best and above the nominal confidence level, although their corresponding width is also the largest. The width of CIs from NoisyCov is similar to those from bootstrap, but the coverage of NoisyCov is unsatisfactory. As the true covariance values are much larger than all the widths in Table 4, we are confident that there is a positive correlation between MRKINC and SHELCO in both Quebec and Ontario.
Remark 6.1. DP bootstrap is easy to use since it does not require additional effort to construct the CIs for either mean or covariance. The NoisyVar method is designed only for mean estimation, and we need a much more detailed analysis to adopt the same idea in covariance estimation.

Conclusion
The analysis of the DP bootstrap provides a new perspective on resampling in DP by considering the output distribution as a mixture distribution which gives a tractable lower bound for the DP bootstrap in f -DP. Furthermore, the composition result for the bootstrap Gaussian mechanism gives a simple asymptotic privacy guarantee and highlights the minimal cost of the DP bootstrap. We use deconvolution to recover the non-private sampling distribution, compute CIs, and perform statistical inference. Our simulations and experiments show that the confidence intervals generated by the deconvolved distributional estimate achieve the nominal coverage.
One direction of future work is on the privacy analysis of the DP bootstrap: The lower bound could be tightened by considering all α i jointly rather than individually. Furthermore, our composition result is asymptotic, and it remains to develop a computable finite-sample composition result.
For the statistical inference using the DP bootstrap, the choice of the number of boot-strap samples B can be further optimized for a given sample size n and privacy parameter µ. One may pursue finite sample guarantees for the CIs or the distributional estimates obtained by the DP bootstrap deconvolution procedure. Furthermore, to handle non-additive noise mechanisms, new inference techniques are needed. We also notice in Figure 6 that under a strong privacy guarantee, our deconvolution method may still be much wider than the true sampling distribution, indicating that the deconvolution procedure can be further optimized to get a tighter estimate and improve the width of the CIs while maintaining the nominal coverage. Apart from the standard confidence interval and the percentile interval, there are also other inference methods based on the bootstrap estimates such as BC a (bias-corrected and accelerated) (Efron, 1987) or ABC (approximate bootstrap confidence) (Diciccio and Efron, 1992), we leave it as future work to implement the private versions of these methods.
We can also use the DP bootstrap framework for high-dimensional inference. Our privacy analysis applies to any mechanism, including those multivariate mechanisms, e.g., the Gaussian mechanism on each dimension of the estimate. Future work of this direction is to apply existing methods of multivariate deconvolution (Youndjé and Wells, 2008;Hazelton and Turlach, 2009;Sarkar et al., 2018) and bootstrap (Hall, 1987;Aelst and Willems, 2005) to our DP bootstrap framework.

SUPPLEMENTARY MATERIAL
This supplementary material is organized in the following way. In Section A, we provide a numerical computation way for evaluating the non-asymptotic privacy guarantee of our lower bound under composition. In Section B, we provide proofs for theorems and propositions in Section 3. In Section C, we first give an analysis of our lower bound on bootstrap Gaussian mechanisms; then we give the details to show why the privacy analysis by Brawner and Honaker (2018) is incorrect, finally, we give the details to show why the privacy loss distribution given by Koskela et al. (2020) is incorrect. In Section D, we analyze the inference of population mean with unknown variance and improve the normal confidence interval (Brawner and Honaker, 2018) using a t-confidence interval with adjusted degrees of freedom. In Section E, we represent our results in algorithmic format for easier practical usage.

A Exact evaluation of f -DP under composition
In this section, we first restate the definition of (ε, δ)-DP with some discussions. Then we provide a numerical computation method for f -DP under composition, which can be used to evaluate the non-asymptotic privacy guarantee of our lower bound under composition. For any M and ε ≥ 0, there exists δ ≤ 1 such that M is (ε, δ)-DP. The privacy profile , δ M (ε), is the smallest δ such that M is (ε, δ)-DP.
, also known as the convex conjugate of f .
There are also other composition results for (ε, δ)-DP.

B Proofs for Section 3
In this section, we provide the proofs for the theorems and propositions in Section 3.

B.1 Proofs for Section 3.1
We first restate the (ε, δ)-DP result in . Then we provide some useful results for the proof of Proposition 3.1.
Theorem B.1 (Theorem 10 in ). In the (ε, δ)-DP form, we use δ = δ M,i (ε) to denote the group privacy of M with group size i. Given ε ≥ 0 and ε = log(1 + Definition B.2. Let f be a tradeoff function,x = inf{x ∈ [0, 1] : −1 ∈ ∂f (x)}. The symmetrization operator which maps a possibly asymmetric tradeoff function to a symmetric tradeoff function is defined as We can use f -DP to restate the upper bound of δ M•boot (ε ) in Theorem B.1: Let From the Supplement to (Dong et al., 2022) (Page 42), we know the following two equations, ε = log(1 − p + pe ε ), δ = p 1 + f * (−e ε ) can be re-parameterized into δ = Id. Since f p may not be symmetric, we use Proposition et al., 2022), we have C p (f ) = Symm(f p ). Therefore, we have an f M•boot -DP result from (ε, δ M•boot (ε ))-DP result: where f M,i is the group privacy with group size i.

From the fact that
Since f = f −1 and f is strictly decreasing when f > 0, we know that f (y) > 0 and f (z) = f −1 (f (y)) = y for y ∈ [0, f (0)). We also know that f (z) = y when y = f (0) since z = f (y) = f (f (0)) = 0. From the fact that C is a sub-differential value of f at α i , we have through the analysis below.
Proof of Lemma B.4. For part 1, first we show that for every α ∈ (0, 1), there exists C ∈ (−∞, 0] such that α ∈ A(C). Since each f i is convex and non-increasing, its sub-differential ∂f i (α i ) is in (−∞, 0] and non-decreasing with respect to α i . Therefore, for any −∞ < C 1 < C 2 ≤ 0, we have a 1 ≤ a 2 ∀a 1 ∈ A i (C 1 ), a 2 ∈ A i (C 2 ), and for any 0 < a 1 < a 2 < 1, we have C 1 ≤ C 2 ∀C 1 ∈ ∂f i (a 1 ), C 2 ∈ ∂f i (a 2 ). We name these two properties as the monotonicity of the sub-differential mapping.
Since each f i is convex and continuous, its sub-differential ∂f i (α i ) is also continuous in the sense that for any ε > 0, there exists δ > 0 such that (Borwein and Vanderwerff, 2010)).
From the continuity and monotonicity of the sub-differential mapping, we know that for any C ∈ (−∞, 0], A i (C) is a closed interval, e.g., [a, b]. Note that A i (C) is always if C ∈ C i,range , there must be an α ∈ (0, 1) such that C ∈ ∂f i (α). By the same reasoning, we also have (0, 1) ⊂ C∈(−∞,0] A i (C) for any i. As A(C) = { k i=1 p i α i |α i ∈ A i (C)}, we have (0, 1) ⊂ C∈(−∞,0] A(C), and we also have the monotonicity of A(C) with respect to C.

Next, we show that mix({(p
) is a well-defined function, i.e., for a given α, although there could be multiple choices of {α i } k i=1 such that k i=1 p i α i = α, we will obtain the same value of k i=1 p i f i (α i ) for all those choices. We consider two choices, {α i } k i=1 and , that both correspond to the same α. We let C and C correspond to {α i } k i=1 and p i α i and p i > 0 for i = 1, 2, . . . , k), from the monotonicity of the sub-differential mapping, we know that C = C since C ≤ C and C ≤ C.
We also know f is non-increasing because of the monotonicity of A(C), the monotonicity of the sub-differential mapping, and f i being non-increasing.
Now we prove that f is convex. By the definition of convexity, we only need to show that for any α, α , t ∈ [0, 1], we have tf (α)+(1−t)f (α ) ≥ f (tα+(1−t)α ). From the construction corresponding to the α, α , andα = tα + (1 − t)α respectively, with their matched sub-differential being C, C , andC. Then we have Therefore, f is convex, and we have proved that f is a tradeoff function.
From Proposition B.6, we know that 1 is strictly decreasing because g is a tradeoff function. Therefore, for any g(α) > 0, there is a one-to-one mapping between α and g(α), and g −1 (g(α)) = α. Now we can view Since the mixture operation is well-defined, we know that g = g −1 .
Theorem 3.3. For i = 1, 2, . . . , k, let f i be a tradeoff function and M i : be a mechanism satisfying f i -DP. Let M : X n → Y be a mixture mechanism which randomly selects one mechanism from k mechanisms, {M i } k i=1 , with corresponding probabilities where k i=1 p i = 1, and the output of M will be the output of M i if M i is selected. Then Proof of Theorem 3.3. Consider the neighboring datasets D 1 , D 2 , and a rejection rule ψ giving the type I error α = E M(D 1 ) ψ, and type II error β = E M(D 2 ) (1 − ψ). We write In order to obtain the lower bound of β given α, which we denote as f min (α), we not only need the tradeoff between α i and β i , which is f i , but also consider the tradeoff between This is a convex optimization problem since the objective function is linear and the constraints f i are all convex (Proposition B.5). Therefore, by Karush-Kuhn-Tucker theorem (Boyd and Vandenberghe, 2004), let the Lagrangian function be if and only if the following condi-tions are satisfied Stationarity: p i − µ i = 0, i = 1, 2, . . . , k; (4) Primal feasibility: Dual feasibility: Complementary slackness: Therefore, we have µ i = p i = 0 from the stationarity condition, and f i (α i ) = β i from the complementary slackness condition.
For a given λ, • if ∂f i (α i ) > λ for all α i ∈ (0, 1), we define α λ,lower i = α λ,upper i = 0; Intuitively, we do not want to have any of the type I error of α on the f i part since the corresponding type II error would be larger otherwise; • if ∂f i (α i ) < λ for all α i ∈ (0, 1), we define α λ,lower i = α λ,upper i = 1; Intuitively, we want to have as much of the type I error of α on the f i part as possible since the corresponding type II error would be larger otherwise. Therefore, for any α ∈ [0, 1], we can find λ such that α ∈ [α λ,lower , α λ,upper ], and we can determine α i by the λ. 4 From the procedure above, we know , and M satisfies f min -DP.
We will use Lemma B.7 in the proof of Theorem 3.4.
Lemma B.7 (Equation (13) in (Dong et al., 2022)). For a symmetric tradeoff function f , Theorem 3.4. Let f i be a symmetric tradeoff function for i = 1, 2, . . . , n. For a mechanism M satisfying group privacy f i -DP with group size i, M • boot satisfies f M•boot -DP where Proof of Theorem 3.4. We are going to find a lower bound of T M(boot(D 1 )),M(boot(D 2 )) uniformly for any neighboring datasets D 1 = (x 1 , x 2 , · · · , x n ) and D 2 = (x 1 , x 2 , · · · , x n ) (without loss of generality, we let x 1 be the different data point in D 1 and D 2 ). We use boot i to denote the conditional bootstrap subsampling where x 1 or x 1 are drawn for exactly i times from D 1 or D 2 . Note that both boot i (D 1 ) and boot i (D 2 ) are random variables for any i = 0, 1, 2, · · · , n. Furthermore, we define boot > to denote the conditional bootstrap subsampling where x 1 or x 1 is drawn for at least once from D 1 or D 2 .
From Theorem 3.3, we know mix({(q i , g)} i∈I ) = g for any g, I and q i . Therefore, we have In order to obtain a better lower bound of T M(boot(D 1 )),M(boot(D 2 )) , we find the mixture of f > and f 0 with considering the bootstrap resampling context because M • boot is a mixture of M • boot > and M • boot 0 .
Consider a rejection rule ψ. Let We prove that there are two additional constraints when we calculate the mixture of f > and f 0 : Intuitively, this is because the hypothesis testing between M(boot 0 (D 1 )) and M(boot > (D 1 )) is similar to the hypothesis testing between M(boot > (D 2 )) and M(boot > (D 1 )).
We consider the hypothesis testing between M(boot 0 (D 1 )) and M(boot > (D 1 )). Similar to the idea in , we consider replacing each x 1 in M(boot > (D 1 )) with a data point uniformly randomly drawn from (x 2 , x 3 , . . . , x n ). We denote the distribution of M(boot > (D 1 )) as ω 1 , the distribution of M(boot 0 (D 1 )) as ω 0 , the distribution of M(boot(D 1 )) as ω 0&1 , and the replacement procedure as replace. Since Balle et al.
For one element in ω 1 , let its histogram be h > = (h 1 , h 2 , · · · , h n ) where h i is the number of occurrences of x i in this element, h 1 ≥ 1, h i ≥ 0 for i = 2, 3, . . . , n, n i=1 h i = n. For one element in ω 0 , we let its histogram be h 0 = (h 1 = 0, h 2 , · · · , h n ). Then Now we replace the h 1 elements (x 1 ) in a point h > of ω 1 by randomly drawing one of (x 2 , x 3 , · · · , x n ), where the total histogram of the replacement is h = (h 2 , h 3 , · · · , h n ), then the new distribution is Then the final distribution of such replacement is It is not easy to show that P replace(ω 1 ) (H = h 0 ) = P ω 0 (H = h 0 ) because of the re-quirement n i=2 h i = h 1 on h . So, we use another way: since for any element in ω 0 , the replacement does not change it, i.e., replace(ω 0 ) = ω 0 , we can perform the replacement on the distribution of ω 0&1 . Therefore, if we show replace(ω 0&1 ) = ω 0 , then Therefore, the replacement procedure indeed transforms ω 1 into ω 0 . Now we are ready to prove β 0 ≥ f > (α > ) and β > ≥ f > (α 0 ).
From the result replace(ω 1 ) = ω 0 , if we consider the tradeoff function between the two distributions M(boot 0 (D 1 )) and M(boot > (D 1 )), we can break each of the two mixture distributions into parts following the replacement procedure so that there is a one-to-one mapping between the parts. When the number of occurrence of x 1 is h 1 in the outcome of boot > (D 1 ), two parts in such a mapping pair have distance h 1 (between boot h 1 (D 1 ) and replace(boot h 1 (D 1 ))), so the tradeoff function between M(boot h 1 (D 1 )) and M(replace(boot h 1 (D 1 ))) is f h 1 . We use Theorem 3.3 to obtain the mixture of those f h 1 . Since f h 1 only depends on h 1 , its corresponding probability in boot > (D 1 ) is p h 1 1−p 0 . Therefore, the tradeoff function between M(boot > (D 1 )) and M(boot 0 (D 1 )) is From this result, and also recall that Now we have established the additional constraints for the mixture of f 0 and f > . We are ready to derive the final mixture tradeoff function.
Complementary slackness: Because our constrained optimization problem has both the convex objective function and convex constraints, we only need to find a solution of α 0 , β > satisfying the KKT conditions, and this solution will be a minimizer of the problem given fixed α.
-The stationarity conditions also hold: We plug in the value C for ∂f > (α), then the right-hand sides of both conditions are 0.
-The stationarity conditions also hold: We plug in the value −1 for ∂f > (ᾱ), then the right-hand sides of both conditions are 0.
-The primal feasibility is satisfied.
-The complementary slackness is also satisfied because -The dual feasibility conditions hold because p 0 ∈ (0, 1).
Therefore, we have the tradeoff function f min from f > as follows Now we verify the f min above is the same as C 1−p 0 (f > ). In Lemma B.7, we replace p and f with 1 − p 0 and f > respectively: is non-increasing, we have that (f > ) 1−p 0 (t) is strictly decreasing because p 0 ∈ (0, 1).

B.3 Proofs for Section 3.2
We first restate the results in (Dong et al., 2022) for comparing the results of f > and C 1−p 0 (f > ) which are used in our proof. Let f 1 , . . . , f n be symmetric tradeoff functions such that κ 3 (f i ) < ∞ for all 1 ≤ i ≤ n.
Corollary B.11. Let M : X n → O be a mechanism satisfying µ-GDP. Then M • boot Now we provide the exact representation of f > in Corollary B.11.
Then f > is the tradeoff function between two Gaussian mixtures, n Proof of Theorem B.12. First, we know that is the cumulative distribution function of the standard normal distribution. We first obtain . We let the type I error and type II error in f i be α i and β i , and we have This setting of (α i , β i ) can also be achieved by using the rejection rule φ(x) = I x≥log (−C) for the hypothesis testing between H 0 : x ∼ N − i 2 µ 2 2 , i 2 µ 2 and H 1 : x ∼ N i 2 µ 2 2 , i 2 µ 2 .
Then α and β as type I error and type II error can be achieved by using the rejection rule φ(x) = I {x≥log(−C)} for the hypothesis testing between We have log h 1 (x) Therefore, f > is exactly the tradeoff function between n i=1 p i 1−p 0 N − i 2 µ 2 2 , i 2 µ 2 and Theorem 3.7. Let µ ∈ (0, ∞) be a given constant, and {µ B } ∞ B=1 be a sequence such Then Proof of Theorem 3.7. Let n be the sample size for the bootstrap resampling. From the result of Corollary B.11 and Lemma B.12, the tradeoff function is f Bi,boot = C 1−p 0 (f > ) where p 0 = (1 − 1/n) n and f > is the tradeoff function between two Gaussian mixtures, From the proof of Lemma B.12, we can parameterize the tradeoff function f > (α C ) = β C we can transform the result in Lemma B.10 to ).

C Privacy analysis on bootstrap Gaussian mechanisms
In this section, we first give an analysis of our lower bound on bootstrap Gaussian mechanisms; then we give the details to show why the privacy analysis by Brawner and Honaker (2018) is incorrect, finally, we give the details to show why the PLDs by Koskela et al.

C.1 Our lower bound
In this section, we restate our lower bound on bootstrap Gaussian mechanisms and also explain how the other curves in Figure 3 were derived.
Corollary B.11. Let M : X n → O be a mechanism satisfying µ-GDP. Then M • boot Then f > is the tradeoff function between two Gaussian mixtures, n Based on this, we can calculate our lower bound C 1−p 0 (f > ) which is the final result of Theorem 3.4 in the case of the Gaussian mechanism, and we visualize our result in Figure   3. From the proof of Lemma B.12, our lower bound is easy to evaluate as the tradeoff and we can use Lemma B.7 to transform the curve of f > to C 1−p 0 (f > ).
In Figure 3, we also visualize the tradeoff functions for specific neighboring datasets.
Note that the tradeoff function between M G (D 1 ) and M G (D 2 ) is G µ . By the number of occurrences of a and 1 − a, we have where we are referring to the distribution of the output of M G applied to one bootstrap sample which includes the randomness of M G as well as the randomness of boot. Since the tradeoff function is a lower bound for the curve between type I error and type II error from any rejection rule, we consider the hypothesis tests H 0 : D = D 1 , H a : D = D 2 and the rejection rule {M G • boot(D) < C} (which need not be the optimal rejection rule). The type I and type II errors are α = n i=0 Φ(Cnµ − iaµ) and β = n i=0 Φ(i(a − 1)µ − Cnµ).
We can see that the curve for a = 0 is not lower bounded by 1-GDP which shows that the bootstrap cannot be used for free with the same f -DP guarantee.

C.2 Counterexample of the privacy analysis by Brawner and
Honaker (2018) In this section, we show an example disproving the statement 'bootstrap for free' by Brawner and Honaker (2018). Their result is under zCDP, a variant of differential privacy, based on Rényi divergences: , and dP dQ is the Radon-Nikodym derivative of P with respect to Q.

2020)
In this section, we show that the privacy loss distribution (PLD) of the bootstrap Gaussian mechanism cannot be the one shown in Koskela et al. (2020).
The privacy loss function is defined as below.
Definition C.2 (Definition 3 (Koskela et al., 2020)). Let M : X N → R be a randomised mechanism and let X Y . Let f X (t) denote the density function of M(X) and f Y (t) the density function of M(Y ). Assume f X (t) > 0 and f Y (t) > 0 for all t ∈ R. We define the privacy loss function of f X over f Y as The privacy loss distribution is defined as below.
Definition C.3 (Definition 4 in (Koskela et al., 2020)). Using the notation of Def. C.2, suppose that L X/Y : R → D, D ⊂ R is a continuously differentiable bijective function. The privacy loss distribution (PLD) of M(X) over M(Y ) is defined to be a random variable which has the density function Although Koskela et al. (2020) did not mention in this definition of the privacy loss distribution, their density function requires L X/Y (s) to be a monotonically increasing function with respect to s. When L X/Y (s) monotonically decreases, we need to replace . Then they state the following result on privacy profile.
Proposition C.5 (Section 6.3 in (Koskela et al., 2020)). Consider the sampling with replacement and the -neighbouring relation. The number of contributions of each member of the dataset is not limited. The size of the new sample is m. Then , the number of times the differing sample x is in the batch, is binomially distributed, i.e., ∼ Binomial(1/n, m), where n is the size of the original sample. Then the subsampled Gaussian mechanism, where B is a uniformly randomly drawn subset (with replacement) of D = {x 1 , . . . , x n } and f (x) 2 ≤ 1 for all x ∈ X, satisfies (ε, δ(ε))-DP After restating the results by Koskela et al. (2020), we prove that the privacy profile δ(ε) from the above theorem is not a valid privacy guarantee for some neighboring datasets.
For f D 1 ,1 vs f D 2 ,1 , we know that f D 2 ,1 f D 1 ,1 is monotonically decreasing since the ratio between each of the three parts in f D 2 ,1 and f D 1 ,1 is decreasing. Therefore, log( and we let it be L X/Y with f X = f D 2 ,1 and f Y = f D 1 ,1 . To calculate the δ, we find t 1 such that t 1 = L −1 X/Y (1) and t 2 = L −1 X/Y (∞), i.e., upper bound of δ(ε) is not true.

D Inference of the population mean with unknown variance
In this section, first we analyze the inference of population mean with unknown variance to obtain the results in Table 1, then we improve the normal confidence interval used in (Brawner and Honaker, 2018) by a t-confidence interval with an adjusted degrees of freedom.

D.1 Rate analysis of the unbiased estimates
We denote the original sample by D = X = (x 1 , x 2 , . . . , x n ). Assume that x i ∈ [0, 1], E(x i ) = θ and Var(x i ) = σ 2 x for i = 1, 2, . . . , n. We let g(D) =X = 1 n n i=1 x i , and its 2 -sensitivity is 1 n . Note that for the analysis in this subsection, we do not require any normal assumption or approximation.
We estimate θ byθ 1 = g(D) when there is no privacy constraint. If µ-GDP is required, we use Gaussian mechanism on g(D) and releaseθ 2 = g(D) + ξ, ξ ∼ N (0, 1 µ 2 n 2 ) as an estimate of θ. For the original bootstrap method which is non-private, we still estimate θ In our analysis below, the mean of B bootstrap estimates isX = 1 The sample variance of B DP bootstrap estimates We have For i = 1, 2, 3, 4, we need Var(θ i ) to construct confidence intervals of θ based on θ i . Since we do not have the population variance σ 2 x , we replace σ 2 x with its unbiased estimates.
Var(s 2 B ) = E(Var(s 2 B |X)) + Var(E(s 2 B |X)) Therefore, Var(Var (θ 3 )) = Var n n − 1s If we set B ∈ O(n), the rate of the width of confidence intervals from the DP bootstrap max(0, Var (θ 4 )) will be the same as those from the non-private bootstrap Var (θ 3 ). If we have B ∈ Θ(n), the rate of the width is O( 1 n 3 ) for all four methods.

D.2 From normal CI to t-CI
Following the discussion above, we useσ 2 x = n 2 n−1s 2 B − (2−2/e)B µ 2 (n−1) to construct confidence interval. Sinceσ 2 x does not exactly follow a χ 2 -distribution, we answer two questions and build our DP bootstrap-t confidence interval: 1) how to choose B so thatσ 2 x is non-negative with high probability; 2) how to modify the degrees of freedom (df) of the t distribution to obtain a better approximation.
We replace n − 1 and B − 1 by n and B respectively, and we have (1 + (2 − 2/e)B nµ 2 σ 2 To solve this with a real number B, we need nµ 2 σ 2 x ≥ 144(1 − 1/e) ≈ 91.025, and we obtain a range of B as Since σ x is unknown, the above inequality of µ and B is only an oracle guidance of choosing the privacy guarantee and the number of bootstrap estimates when constructing the standard confidence interval using the unbiased estimate.
As we have assumed that (B − 1)s 2 B /( n−1 n 2 σ 2 x + (2−2/e)B n 2 µ 2 ) follows a χ 2 -distribution, we want to build an approximated t-confidence interval which contains the uncertainty from bothθ 4 and Var (θ 4 ). For simpler illustration, we consider the case that Var (θ 4 ) > 0. We aim to analyze the distribution ofθ 4 −θ Var (θ 4 ) . Its distribution depends on n, µ, B, σ 2 x which is not a t-distribution, but sinceθ 4 − θ is independent to Var (θ 4 ) when conditioning on X, we try to use a t-distribution as an approximation. Now we need to find a proper degrees of freedom (df) for this approximated t-distribution. Since t n = z √ χ 2 n /n , this df n is determined by the χ 2 n -distribution in the denominator. If C ∼ χ 2 n , we know that Var(C/E(C)) = 2 n . Therefore, we can usedf = 2/Var(C/E(C)) to estimate the df of an approximated χ 2 -distribution.
Here we useσ 2 x ands 2 B in the above approximation since they are known to us. For the private bootstrap, we letdf := (σ 2 x /(ns 2 B )) 2 (B − 1) and use the fact thatθ 4 −θ Var (θ 4 ) follows td f -distribution to build confidence interval for θ.
We observe that this t-confidence interval with adjusteddf improves the normal confidence interval by Brawner and Honaker (2018) when the choice of µ and B satisfies the conditions in Eq. (14) and (15). Our analysis also reveals the reason for the under-coverage issue in the simulation by Brawner and Honaker (2018) when B is too large. However, when µ is too small, i.e., we have a strong privacy constraint, the under-coverage issue still exists since we do not always have Var (θ 4 ) > 0. Therefore, we develop the deconvolution method for inference.

E DP Bootstrap algorithm for the Gaussian Mechanism
In this section, we show our DP bootstrap method in Algorithm 1 for the Gaussian Mechanism, and based on it, we can obtain the DP sampling distribution estimate (Algorithm 2) and DP point estimate with corresponding standard error (Algorithm 3). We also include the NoisyVar (Algorithm 4) and NoisyCov (Algorithm 5) used in our simulations and experiments.
Note that our privacy analysis in this paper (on Algorithm 1) is valid for g : X n → R d for any d ∈ {1, 2, 3, . . .}, but the deconvolution (Algorithm 2) can only be applied when d = 1 because of deconvolveR. For d ∈ {2, 3, . . .}, one may try other deconvolution methods. This will not affect the DP guarantee because of the post-processing property. For Algorithm 3, we only show the case when d = 1, but one can generalize it to d ∈ {2, 3, . . .} intuitively. Similar ideas in Algorithm 3 are also used in (Brawner and Honaker, 2018). Du et al. (2020) proposed many algorithms in their paper for constructing DP CIs. We choose NoisyVar because other algorithms assume the data follows a normal distribution and use this property to estimate the standard deviation from quantile estimates which is restrictive. Du et al. (2020) used NoisyVar with a parametric bootstrap framework, SIM; Therefore, we put them together in Algorithm 4. We adopt their idea in Algorithm 5, named NoisyCov, to construct DP CIs for the covariance where we use a normal distribution to approximate the correlation structure of the clamped data, i.e., the original data restricted to a specific range for the sensitivity bound in the privacy guarantee. Note that in NoisyCov, we do not spend any privacy budget estimating the mean since it does not affect the covariance in the parametric bootstrap simulation.
Algorithm 2 DP Sampling Distribution (Differentially private sampling distribution estimation) 1: Input Dataset D ∈ X n , Statistic g : X n → R with 2 sensitivity ∆ g , µ > 0 for µ-GDP, number of bootstrap samples B.
Algorithm 3 DP Estimate with SE (Differentially private point estimate with corresponding standard error) 1: Input Dataset D ∈ X n , Statistic g : X n → R with 2 sensitivity ∆ g , µ > 0 for µ-GDP, number of bootstrap samples B.