Exact Inference with Approximate Computation for Differentially Private Data via Perturbations
Ruobin Gong


  
Exact Inference with Approximate Computation for DP Data]Exact Inference with Approximate Computation for Differentially Private Data via Perturbations
  

R. Gong]Ruobin Gong
Department of Statistics, Rutgers University, New Brunswick, NJ 08854
rg915@stat.rutgers.edu





This paper discusses how two classes of approximate computation algorithms can be adapted, in a modular fashion, to achieve exact statistical inference from differentially private data products. Considered are approximate Bayesian computation for Bayesian inference, and Monte Carlo Expectation-Maximization for likelihood inference. Up to Monte Carlo error, inference from these algorithms is exact with respect to the joint specification of both the analyst's original data model, and the curator's differential privacy mechanism. Highlighted is a duality between approximate computation on exact data, and exact computation on approximate data, which can be leveraged by a well-designed computational procedure for statistical inference.


[
    [
    December 30, 2023
=====================







	First version: September 2019     This version: September 2022





§ INTRODUCTION



Differential privacy <cit.> advances statistical disclosure limitation by putting forth a formal and practical framework. In addition to grounding the concept of privacy on a mathematical footing, differential privacy distinguishes itself from traditional approaches by offering transparent probabilistic mechanisms, whose specifications can be made public without sabotaging the privacy guarantee. Differential privacy has been adapted by major data curators in the industry, research and government sectors. 
As a prime example, the U.S. Census Bureau deploys differential privacy to protect the 2020 Decennial Census data products <cit.>. The P.L. 94-171 redistricting data files have been released on August 12, 2021 <cit.>.


In this work, we adopt the perspective of a data analyst operating under the dissemination mode of data access <cit.>. A data curator, such as the Census Bureau, collects potentially sensitive data and releases differentially private data products to the analyst. The analyst in turn conducts statistical inference for their quantities of interest based on the privatized data. The analyst's goal is to draw trustworthy inference from the statistical model they wish to fit, knowing that the data have undergone privacy protection. This may not be a trivial task. The curator instills differential privacy in the data product via a data processing mechanism. Naïvely treating processed data as if they are unprocessed may result in erroneous and misleading statistical inference. With the wide adoption of differential privacy for disclosure limitation, social scientists and policy researchers are faced with the challenge to revise their preferred statistical analyses to account for the privacy mechanism, however complex they may be. To keep up with advances in privacy protection, we  need theoretically sound and computationally efficient statistical methodologies to supplant their predecessors <cit.>. 







This paper discusses the adaptation of two classes of approximate computation algorithms, approximate Bayesian computation  () and Monte Carlo Expectation-Maximization (MCEM), to obtain exact Bayesian and likelihood statistical inferences based on differentially private data products. The word exact means that, up to Monte Carlo error, the resulting inference corresponds precisely to the joint statistical model that accounts for both the analyst's specifications and the differential privacy mechanism. This paper draws a concrete connection between the novel disclosure limitation mechanisms that obey differential privacy, and the vast reserve of computational strategies available for likelihood and Bayesian statistical inference. The hope is that users of traditional, non-differentially private data can smoothly transition their existing methodologies to suit novel, differentially private data products while maintaining statistical validity. The two methods discussed in this work are applicable to a wide range of existing models, dispensing the need to analytically recompute the new joint model to account for the privacy mechanism.  Both classes of algorithms discussed in this paper do not assume specific structures of the likelihood, prior, and privacy mechanism. Indeed, the likelihood approach only requires that the analyst's original model is suitable for EM, and the Bayesian approach only requires that the original likelihood can be simulated and that the prior is proper. Should specific and convenient model structures be available, the proposed mechanisms would be amenable to adaptation to enhance computational efficiency.


The remainder of this paper is organized as follows. Section <ref> lays out the mathematical formalism and notation for differential privacy and perturbation mechanisms. Section <ref> proposes a rejection  algorithm, and shows that with kernel and bandwidth chosen to correspond to the perturbation mechanism underlying the privatized data, it produces exact posterior inference in the form of independent and identically distributed samples from the true posterior distribution. Section <ref> discusses an importance sampling implementation of Monte Carlo EM for likelihood inference. The validity of both approximate computation methods owes to the fact that their tuning elements can be chosen in accordance with the differentially private perturbation mechanism that is used to generate the privatized data product. Section <ref> provides two numerical demonstrations, on the Bayesian and likelihood inference for privatized count, and a differentially private adaptation of the Lalonde dataset for inference on job training program efficacy. Section <ref> concludes with a discussion on the duality between approximate computation on exact data and exact computation on approximate data, and the various challenges to the efficiency of these proposals.

 



§ DIFFERENTIAL PRIVACY AND PERTURBATION MECHANISM


Differential privacy aims to protect the confidential information of individual respondents in a dataset ∈, without undue sacrifice of accuracy in learning about aggregate features of the underlying population as represented by . Here, an aggregate feature is a query : →^p, a deterministic function of , such as the sample average, variance, quantiles and so on. Queries are the means through which analysts learn from the dataset. Counting queries, including histograms and contingency tables which are ordered multivariate counts over a partition of , constitute a most useful class of queries. It is the main query type for the 2020 U.S. Census data products, tabulated across various geographic levels such as states, counties, and Census blocks.


Differential privacy is realized via a probabilistic mechanism based on the intended query. A differentially private query reflects as truthful as possible about the status of , meanwhile behaves similarly should it be calculated based on any neighboring dataset of . The notion of differentially privacy is defined in probabilistic terms.



A random function :→^p is (ϵ,δ)-differentially private if for all neighboring datasets (,') ∈^2 and all A ∈ℬ(^p), 

    Pr((') ∈ A) ≤ e^ϵ· Pr(() ∈ A) + δ.
 
 is ϵ-differentially private if it is (ϵ,0)-differentially private.



The pair (, ') ∈^2 constitutes neighboring datasets if they differ by precisely one entry, either by adding or dropping one respondent, or by taking a different value (as used in the definition of bounded differential privacy; ). When operating on neighboring datasets, the random function  induces pairs of probability measures, associated respectively with  () and ('), that are close to each other. The degree of closeness is controlled by the privacy loss budget ϵ and δ. In the extreme case that both are zero, the two measures must be equal on every Borel set A, which for general  can only happen if  does not depend on the data at all. In other words, differential privacy requires that the distribution of  to be stable within the small neighborhood around the observable dataset.
 

Differential privacy is a property pertinent to the random function . Many widely employed differential privacy mechanisms take the form of perturbation mechanisms.


For a dataset ∈ and a deterministic function : →^p, the random function  is a perturbation mechanism based on  if

    () |() ∼η_( · |() ),

for η_ a known conditional probability distribution. In particular,  is an additive perturbation mechanism based on  if

    () = () + h,

where the noise component  is a p-dimensional random variable with known distribution η, 𝔼() =  0, and h > 0 is a scale  (or bandwidth) parameter.	


The differentially private query  is a noisy version of its deterministic counterpart . The protection of privacy is achieved through randomly perturbing what would otherwise be a deterministic query calculated based on . The subscript “” in  emphasizes that it instantiates the privacy  mechanism η_, rather than the data generation mechanism of , as the analyst might posit. The perturbation mechanism embodied by  is said to be unbiased if it satisfies 𝔼(() |() )=(). Additive perturbation mechanisms, by construction of (<ref>), are unbiased. Furthermore, if the scale parameter h does not depend on the confidential dataset , the mechanism may be called a data-independent mechanism <cit.>. Note that the additive perturbation mechanism resembles the classical measurement error model <cit.>, where the noise h has a known distribution, and the noisy measurement  is observed precisely once.  Appendix <ref> gives three examples of widely used differentially private mechanisms, with additive perturbation using Gaussian and Laplace noises. Their definitions invoke three notions of functional sensitivity, (<ref>)-(<ref>), which we generally denote as Δ(), to capture the idea that certain  is more revealing of individual information in  than others. It is crucial that the scale parameter of the additive perturbation mechanism is chosen as a function of both the sensitivity of  and the privacy budget, i.e. h = h(ϵ,δ,Δ()). Additional examples of additive differentially private mechanisms include the generalized Cauchy <cit.>, double Geometric <cit.>, correlated multivariate Gaussian <cit.> and the k-norm <cit.> mechanisms. Examples of non-additive perturbation mechanisms include the randomized response mechanism <cit.>, exponential mechanism <cit.>, objective perturbation <cit.>, among others.







A primary strength of differential privacy over traditional disclosure limitation frameworks is its transparency, which means that the specification of the perturbation mechanism η_ may be fully revealed to the data analyst (and indeed the public) while keeping the privacy guarantee intact. For additive mechanisms, this specification consists of 's distribution η, scale parameter h, and the privacy loss budget ϵ and δ. Perturbation mechanisms can be correctly accounted for in the probabilistic modeling of privatized data. Despite the necessary sacrifice of statistical efficiency, likelihood and Bayesian models utilizing privatized data can still retain validity, in the sense that any inference drawn based on  can still be drawn based on _ correctly while accounting for its generative process. As Section <ref> will discuss, for Bayesian analysis, an  rejection algorithm guarantees the exactness of draws from the true posterior distribution, when properly tuned according to the parameters of the perturbation mechanism. The  nature of the privatized query makes  an appealing choice for posterior computation, even when the model is not as complex as to necessitate its use.













§ EXACT BAYESIAN INFERENCE WITH DIFFERENTIALLY PRIVATE DATA


In the absence of privacy protection, suppose a Bayesian model was posited based on the confidential query  as a function of . Let  ()|θ∼π(|θ) be the confidential data likelihood, and θ∼π_0(θ) the prior distribution for θ. The posterior distribution of θ given  is 

    π(θ|) ∝π_0(θ)π(|θ).

If the query  isn't privacy-protected, quantities calculated based on (<ref>), either analytically or via simulation, would conclude the Bayesian analysis. With the privacy protection mechanism in place, however, we no longer observe the confidential query , but rather the  privatized (perturbed) query _ as a single realization of the privacy mechanism (<ref>). The joint distribution of θ and _ is 

    π(θ,_)=∫π(θ,)η_(_|)d,

marginalized over the latent . This identity holds because the conditional distribution of _ given  and θ is free of θ, as it is precisely the known perturbation mechanism: π(_|, θ) = η_(_|). The posterior distribution of θ given _ is

    π(θ|_) = ∫π(,_,θ)/π(_)d = π_0(θ)∫η_(_|)π(|θ)d/∫π_0(θ)∫η_(_|)π(|θ)d dθ.

As (<ref>) is the true posterior distribution for θ given the observable information,  analytical or simulated computation based on (<ref>) would conclude the exact Bayesian analysis. However, computation of (<ref>) may not be trivial, as part of it involves the observed likelihood ∫η_(_|)π(|θ)d, which is an integral of the product between the confidential data likelihood and the privacy mechanism. The challenge is exacerbated by the fact that the confidential likelihood is specified by the data analyst, whereas the privacy mechanism is specified by the data curator. These choices are typically made independent of one another, and either of them may already be complex and computationally demanding on its own. 

Algorithm <ref> presents a recipe to generate independent and identically distributed samples from the exact posterior distribution (<ref>). It demands little of the tractability of the confidential likelihood. The only requirement is that for given values of θ, one can simulate from π(|θ), but otherwise it need not be available in closed form. Algorithm <ref> is a type of  algorithm, which was designed to supply practical solutions to large-scale models for which the likelihood may be implicit or intractable and posteriors without closed-form expressions.  brought computational feasibility to stochastic differential equation models for complex dynamic systems in population genetics <cit.>, systems biology <cit.> and ecology <cit.>, albeit  posteriors are typically only approximate relative to the true target posterior. However, as will be shown in Theorem <ref> and discussed in Section <ref>, the employment of  for differentially private data serendipitously eradicates the “approximate” nature of the resulting posterior samples, which otherwise would be the case if the data were noise-free.
 







 
Let π(|θ) be the likelihood for the unobserved confidential query , π_0(θ)  a proper prior distribution, and  η_(_|) a perturbation mechanism. Then, Algorithm <ref> samples independently and identically from the exact posterior distribution π(θ|_) defined in (<ref>).



Proof of Theorem <ref> can be found in Appendix <ref>. Key to the validity of  Theorem <ref> is that the differentially private perturbation mechanism is ignorable for θ <cit.>, or in other words, the unobserved confidential query  is sufficient with respect to the complete likelihood π(,_|θ).  
Traditional statistical disclosure limitation mechanisms may or may not enjoy ignorability, a matter further complicated by their non-transparency to impact the quality of downstream statistical analysis  <cit.>. By contrast, the ignorability property of differential privacy enables exact statistical inference and may substantially simplify the computational task.

An intuitive connection with traditional  can be drawn if we restrict attention to the case of additive perturbation. As defined in (<ref>), assume η_(_|) = η((_-)/h) where η(·) is the density of the additive noise  and h a scale parameter, both known precisely to the analyst. Algorithm <ref> adopts the kernel density η, properly scaled by a factor of c, with bandwidth h = h(ϵ,δ,Δ()) and center _ to be its acceptance probability at step 3, thus reduces to a classic rejection  algorithm with tuning parameters (i.e. kernel and bandwidth) set to match precisely the additive perturbation mechanism employed to generate _.



One way to understand Theorem <ref> is that the privacy mechanism plays the role of the “random summary statistic” in the noisy  algorithm of <cit.>. Noisy  is calibrated with respect to the joint Bayesian model, whereas  typically isn't. However, the kernel and bandwidth in noisy  are merely parameters to fine-tune the tradeoff between approximation error and the Monte Carlo error in the posterior, which in turn controls the efficiency of the sampler. In contrast, both the kernel and the bandwidth of Algorithm <ref> are dictated externally by the perturbation mechanism and the privacy loss budget. The computational tradeoff and the privacy tradeoff are “bundled” together: specifying the parameters of  also specifies those of the privacy mechanism, and vice versa.  




The overall acceptance probability of Algorithm <ref> is π(_)/maxη_, or the model evidence evaluated at _ divided by the modal density of η_  (see Appendix <ref>). This means that rejection can be frequent if model evidence is low, such as when the prior and the observed likelihood are in disagreement (termed prior-data conflict; ), or if the privacy bandwidth h is too small. 


To address the concern, Algorithm <ref> can be adapted to work with a variety of alternative  sampling techniques to produce consistent posterior estimates for functions of interest. As an example, we discuss an importance sampling variation to Algorithm <ref> as follows. At step 1 of each iteration, sample θ_i ∼ g(θ), a proposal distribution that is positive wherever the prior π_0(θ) is positive. At step 3, no rejection is performed, but instead θ_i is assigned a weight 

    ω_i = ω(_i,θ_i) = η_(_|_i)π_0(θ_i)/g(θ_i).

The algorithm returns weighted draws {θ_i, ω_i}_i=1^N. For a square-integrable function of interest a(θ), the weighted average estimator converges in probability to its posterior expectation given _ as N →∞ <cit.>:

    ∑_i=1^Nω_ia(θ_i)/∑_i=1^Nω_ip→𝔼_g(ω(θ,)a(θ))/𝔼_g(ω(θ,))=𝔼(a(θ)|_),
 
where 𝔼_g(·) is with respect to the joint distribution g(θ)π(|θ), and 𝔼(·|_) is with respect to the true posterior in (<ref>). The proposal distribution g(·) can be chosen to minimize the variance of the estimator in (<ref>), such as a density that is close in shape to a(θ)π_0(θ) <cit.>. Further adaptations of and beyond , such as hybrid importance-rejection sampling <cit.>, rejection control <cit.>, Markov chain Monte Carlo <cit.> and sequential Monte Carlo <cit.> can be developed likewise, while the consistency result of (<ref>) remains standing. 




§ EXACT LIKELIHOOD INFERENCE WITH DIFFERENTIALLY PRIVATE DATA


This section discusses a Monte Carlo Expectation-Maximization <cit.> implementation for likelihood inference with differentially private data. Under the classic setting, when a likelihood involves both observed and latent data, EM seeks the maximum likelihood estimate of the parameter by iteratively integrating the log likelihood over the conditional predictive distribution of the latent data given the observed data and a current parameter value (the E-step), and maximizing the parameter value over this integral (the M-step).

In the context of differential privacy, the complete data is (, _), in which the latent data is the confidential query , and the observed data is the privatized query _. In the special case of additive perturbation, _ =  + h is a convolution of  and the noise component . The complete likelihood is written as L(θ;,_) ∝π(,_|θ), as defined in Section <ref>. The EM algorithm for maximum likelihood inference for θ given the differentially private _ is schematically described in Algorithm <ref>. 


The E- and M-steps are iterated until convergence, that is when θ^(t) stabilizes so that its distance (somehow measured) from the previous iteration, dist(θ^(t),θ^(t-1)), is sufficiently small. It is worth noting that the constant term in (<ref>) is equal to 

    𝔼( logη_(_|)|_,θ^(t)),

which can be ignored within the EM algorithm. This is because, as discussed in Section <ref>, the privacy mechanism η_ is known and independent of θ, and so is the conditional predictive expectation of its log density. 











As alluded to in Section <ref>, for likelihood modeling of differentially private data, the confidential data likelihood and the privacy mechanism are typically specified by separate parties without coordination with one another. Thus in general, one cannot expect the observed data likelihood (which is an integral of their product) to come from an exponential family <cit.>, nor be able to perform both the E- and the M-steps analytically. Monte Carlo implementation of one or both steps may be needed, which amounts to implementing the E-step of Algorithm <ref> via an importance sampling scheme. We  describe this scheme in Algorithm <ref>. The set of weighted samples {_i, ω_i}_i=1^N produced by Algorithm <ref> may be used in two ways, depending on whether the confidential data likelihood is or isn't an exponential family. We discuss both cases below.






 §.§ Confidential data with exponential family likelihood.


In the simpler scenario that the confidential data likelihood π(|θ) as specified by the analyst belongs to the exponential family, it admits a sufficient statistic to the parameter θ which we denote as b().  The function Q(θ; θ^(t)) in  (<ref>) can be written as an explicit function of θ and 

    𝔼( b()|_,θ^(t)),
	
the conditional expectation of b() given _ and the current maximizer θ^(t). With this simplification, however, (<ref>) may still not be evaluable in closed form, in which case we utilize the set of weighted samples {_i, ω_i}_i=1^N produced by Algorithm <ref> to consistently estimate it at every iteration t. Indeed, as N →∞, the weighted estimator

    ∑_i=1^Nω_ib(_i)/∑_i=1^Nω_i

converges in probability to (<ref>). For the E-step of the (t+1)st iteration, θ^(t+1) can be found by maximizing Q(θ; θ^(t)), replacing (<ref>) therein with (<ref>). The
 effective sample size at the tth iteration is

    ess^(t)(N) =  Nπ^2(_|θ^(t))𝔼^-1_|θ^(t)( η^2_( _|) ),
 
where the subscript “|θ^(t)” signifies that the expectation is taken with respect to the current approximation to the confidential data likelihood, or equivalently, the proposal distribution of the E-step importance sampler. Derivation of (<ref>) may be found in Appendix <ref>.





In Algorithm <ref>, the _i's are simulated from the current approximation to the analyst's confidential data likelihood, and the weights ω_i's are separately determined by the curator's privacy mechanism. Similar in spirit to Algorithm <ref>, this separation allows the computation to easily accommodate independently derived choices of data likelihood and privacy mechanisms, and does not require the evaluation or integration of quantities that are nontrivial functions of both.
Whenever appropriate, however, one may modify Algorithm <ref> to sample from the conditional predictive distribution in more efficient ways. For example, with rejection or Markov chain-based samplers, _i follows a proposal distribution and ω_i = 1 if _i is accepted and 0 otherwise <cit.>. One may also perform importance sampling where _i ∼π( |_,θ^(t-1)), the approximation to the conditional predictive distribution at the previous iteration, and ω_i = π( |_,θ^(t))/π( |_,θ^(t-1)) the ratio between the current and previous approximations, thereby reweighting and recycling the multiply-imputed _i's to save computational effort <cit.>. One may also resample the simulated _i's according to their associated weights to obtain an unweighted rejection sample, as long as the goal is to construct as accurate as possible an estimate for (<ref>) as part of the E-step.




 §.§ Confidential data with general likelihood.

If the confidential data likelihood does not come from an exponential family, Q(θ; θ^(t)) of (<ref>) may not reduce to a straightforward expression involving θ and (<ref>). In this case, the E-step requires a full approximation to Q(θ; θ^(t)) as a mixture of augmented log likelihoods, constructed as follows. 

Let {_i,ω_i}_i=1^N be a weighted sample from the target distribution π( |_,θ^(t)), the tth approximation to the conditional predictive distribution. Specifically {_i,ω_i}_i=1^N can be the importance sample generated by Algorithm <ref>, or by one of its variations such as as described above. Then,

    Q̂(θ ;  θ^(t)) = m∑_i=1^Nω_ilogπ(_i|θ)

serves as a consistent approximation to Q(θ; θ^(t)). The constant m^-1 = ∑_i=1^Nω_i in (<ref>) is inconsequential to the maximizer in the ensuing M-step, as long as the ω_i's do not involve the unknown parameter θ. That is indeed the case since, again, the perturbation mechanism is ignorable for θ. Writing λ_θ() = ∇_θlogπ(|θ), the observed score function  ∇_θlogπ(_|θ^(t))  can be approximated at the tth iteration according to 

    𝔼( λ_θ() |_,θ^(t))·=  m∑_i=1^Nω_iλ_θ(_i).

The observed Fisher information can also be approximated according to 

    -∇_θ^2logπ(_|θ^(t))  ·=
    
    m∑_i=1^Nω_i{ -∇_θλ_θ(_i)-λ_θ(_i)λ_θ(_i)^⊤} + m^2∑_i=1^N∑_j=1^Nω_iω_jλ_θ(_i)λ_θ(_j)^⊤.

Derivations of the observed score function and observed Fisher information can be found in Appendix <ref>. Both (<ref>) and (<ref>) may be used for quantifying the inferential uncertainty under the normal approximation to the likelihood <cit.>, as well as accelerating and assessing convergence for Newton-type implementations of the M-step. The approximations given above rely only on that the first and second derivatives of the confidential likelihood be evaluable at the simulated s_i's.




For any EM algorithm (and not just Monte Carlo EM) to be applicable to likelihood inference from differentially private data, one must be able to evaluate the confidential data likelihood π(|θ), to the extent that maximization of the Q function can be done at least numerically. 
The vast literature on Monte Carlo EM has much to instruct on implementing both the E- and the M-steps with better convergence rates, sampling efficiency, or under computational capacity constraints, for adapting modeling scenarios to differentially private data. The additive perturbation mechanism of (<ref>) is a special instance of a linear mixed effects model, which is particularly suitable for Monte Carlo EM and has been studied extensively in the literature, e.g. <cit.>. 













§ NUMERICAL DEMONSTRATIONS




 §.§ Bayesian and likelihood inference from privatized count


In this simple example, we consider modeling the number of respondents from a sample  in possession of a certain feature. (·) is the univariate counting query, for which we posit the sampling model ()|θ∼ Pois(θ). θ is the population expectation parameter for which we wish to draw Bayesian and likelihood inference. 

First consider a Bayesian model for θ. We implement rejection  as described in Algorithm <ref> to draw from the exact Bayesian posterior based on the privatized count _. Suppose _  is produced by the ϵ-differentially private Laplace mechanism  (Example <ref> in Appendix <ref>), where the additive noise follows ∼ Lap_p(1) with bandwidth h = ϵ^-1. As with general  samplers, Algorithm <ref> can work with arbitrary choices of prior and likelihood that need not be conjugate, so long as the prior is proper. For the purpose of illustration, we consider the prior θ∼ Gamma(α,β), where α and β are fixed hyperparameters, so that an analytically tractable posterior can be obtained for visual comparison. 









Figure <ref> depicts both the correct and the naïve analyses, with hyperparameters α = 25, β = 1, privacy loss budget ϵ = 0.2, and _ = 37.4. The true analytical posterior (green solid density), normalized via numerical integration, coincides with the differentially private  posterior histogram tabulated from 10^4 draws from Algorithm <ref>. The correct analysis differs substantially from the incorrect naïve  posterior (red dotted density), which treats _ as if it were an observed confidential query. (The latter posterior amounts to the posterior from the standard  Gamma-Poisson conjugate model.) Compared to the correct posterior, the naïve posterior succumbs less to the shrinkage effect imposed by the prior. It assigns a heavier weight of evidence to the observed value of _, more so than it deserves. It is furthermore overly concentrated at the mode, underestimating the posterior uncertainty associated with θ. 
 
Appendix <ref> displays additional experiments that employ Gamma prior distributions with hyperparameters α = 2, 5, 50, 75 and β = 1.
Worth noting is that when the privatized observation _ appears highly unlikely under the chosen prior (or is in conflict with it, in other words), the correct posterior heavily discounts the contribution by _. For example, such is the case when α = 2 or 5, as seen in Figure <ref> (a) and (b):  the correct posterior is in close alignment with the prior and differs significantly from the naïve posterior. As alluded to in Section <ref>, prior-data conflict presents a challenge for  algorithms in general, because forward sampling tends to explore the area with higher prior predictive concentration. A realized observation far from that area would result in a low acceptance rate. To see this, Table <ref> reports the average acceptance rates and their standard errors over 20 direct repetitions of Algorithm <ref> under various choices of Gamma priors. In comparison with the observed query, these priors range from congruent to conflicting, as can be seen from the varied differences between _ and its prior predictive expectation: 𝔼_θ(_) = ∫_∫η_(_|)∫π(|θ)π_0 (θ)dθ d d_.
  

















Maximum likelihood estimation for θ is carried out as follows. The confidential data likelihood is the Poisson density. Importance sampling as described in Algorithm <ref> is used to construct estimates for (<ref>) at every iteration of the E-step, followed by an analytical M-step.  Appendix <ref> describes details of the implementation using three stages of successively more stringent tolerance levels. With θ^(1) = 1, the algorithm converges to the maximizer θ̂ = 37.237, with observed Fisher information estimated to be 1.582× 10^-2. If _ were erroneously treated as the confidential data, the MLE for θ would've been 37.4, and the observed Fisher information would've been 2.674 × 10^-2, or 69% larger than the correct estimate, again displaying an underestimation of inferential uncertainty. The reduction of Fisher information content reflects a loss of statistical efficiency induced by the privatization mechanism, and is expected in typical inference problems whenever confidential data are replaced with their privatized counterparts. Details of the above calculations can be found in Appendix <ref>.





 §.§ Lalonde dataset


The Lalonde dataset <cit.>  was curated from the the randomized trial of the National Supported Work (NSW) Demonstration and nonexperimental comparison data, for the purpose of studying the efficacy of the job training program on recipients' future earnings. The dataset, with a total of 185 treated and 260 control units, is well-studied in the causal inference and econometrics literatures using regression and propensity matching methods, see e.g. <cit.>. We employ the example here to illustrate a Bayesian analysis that compares the 1978 earnings of the treatment and control groups, if ϵ-differentially private versions of the key descriptive statistics were released instead.

Let z_i be the observed indicator for whether subject i received treatment (z_i = 1) or control (z_i = 0), and y_i their earning in 1978 (in $1k).  The full parameter of the model is θ=(τ,μ,σ_t^2,σ_c^2), in which τ is the difference in average earnings between the treatment and control groups, and is the primary parameter of interest.  We posit independent priors for elements of θ, as well as the sampling model

    y_i| z_i,θ∼ N(τ z_i+μ,σ_t^2z_i+σ_c^2(1-z_i)).

For simplicity's sake, we do not consider additional covariates that distinguish the treatment and control subjects. 


Among  the descriptive statistics that the publisher plans to release,  relevant to the inferential task at hand are the within-group sample means and sample variances: = (y̅_t, y̅_c, s_t^2,s_c^2). Together they make up the sufficient statistic for the full parameter θ.  The top row of Figure <ref> displays the posterior inference for θ by repeatedly fitting this model in RStan using the actual value of . Discrepancies among the ten boxplots within each figure, all of them minor, are due to Monte Carlo errors. According to the model, there is a discernible positive treatment effect since the posterior mass of τ is overwhelmingly positive.

Suppose that the data publisher releases ϵ-differentially private version of sample means and variances. Since the mean and the variance are real-valued functions, they do not have a finite global sensitivity Δ_GS as defined in (<ref>), hence the Laplace mechanism cannot apply directly to them. To circumvent this issue, the publisher may clamp the underlying query, that is to enforce its value to stay within a bounded range. For simplicity's sake, suppose that the clamping range on individual income is conservatively set, say to between zero and $100k, and the treatment and control groups are guaranteed to exceed 100 people. This effectively restrict the global sensitivity of y̅_t and y̅_c  to 1 and that of s_t^2 and s_c^2 to 100. For reference, the maximum observed individual income in the dataset is $60.3k, and the treatment and control group are respectively of sizes n_t = 185 and n_c = 260, ensuring that all confidential query values fall well within the clamping range.  The benefit of conservative clamping is that the privatized statistics would not require truncation correction, even though it amounts to an inefficient privacy budget allocation strategy.   Further suppose two separate privacy loss budgets of ϵ = 1/3 and 100/6  are respectively expended on the sample means and variances, through Laplace mechanisms employing independent zero-mean noise components with bandwidths h^-1 = 1/3 for each of the sample means  y̅_t, y̅_c, and h^-1 = 1/6 for each of the sample variances s_t^2, s_c^2.



 
  
  
  
The middle and bottom rows of Figure <ref> respectively display posterior inferences from naïvely fitting the original model (i.e. disregarding the privacy mechanism)  in RStan, and correctly fitting the exact posterior (i.e. accounting for the privacy mechanism) using rejection  of Algorithm <ref>. Both methods were fitted to the same ten independent realizations of _ from the Laplace mechanism. Discrepancies among the ten boxplots within each figure in these two rows are due to the random privacy mechanism and to Monte Carlo errors – the latter to a much lesser extent. We see that with the correct analysis, posterior uncertainty for all parameters are substantially inflated, in part due to the highly inefficient allocation of the privacy loss budget. As a result, we can no longer conclude that the treatment effect is significantly in either direction. However, the posterior quantiles overlap substantially with their counterparts from the original posterior on the top row, indicating that the cost of privacy manifests more as an estimation precision loss rather than bias. This stand in contrast against the naïve analysis which delivers tight, yet idiosyncratically displaced, posterior masses. Details of this analysis can be found in Appendix <ref>.











§ CONCLUSION AND DISCUSSION











Modern likelihood and Bayesian inferences face the challenge of model complexity. They appeal to Monte Carlo and approximate methods to carry out needed computation, even if the resulting inferences are only approximate with respect to the full statistical model. This paper discussed how approximate computation algorithms, specifically  and Monte Carlo EM, can be adapted to obtain exact Bayesian and likelihood statistical inferences based on differentially private data products. In both cases, the tuning elements of the approximate computation algorithms are chosen to accord to the specifications of the differentially private perturbation mechanism, which can be made transparent to the data analyst. Both methods are applicable to a wide range of modeling scenarios, and may help data users transition existing methodologies to apply to differentially private data products while maintaining statistical validity of their analysis.


When no privacy mechanism is involved,  algorithms exhibits a bias whenever it cannot enforce an exact match between the observed and the simulated data <cit.>, which is typically the case in practice. The justification of  relies on that in the limit as the bandwidth h → 0, the  posterior π_(θ|) approaches the true posterior π(θ|) <cit.>. However in practice,  h cannot be too small in order for the algorithm to generate an adequate number of samples, trading off a larger approximation error with a smaller Monte Carlo error. 

The statistical insight underscored by this paper is the duality

    approximate computation on exact data  ↔  exact computation on approximate data.

Differentially private data is approximate data. The perturbation mechanism with which the data were treated serves coincidentally as the attributable cause of the approximation error. When differentially private data are employed, the Monte Carlo error becomes the sole kind of error attributable to the  algorithm, and vanishes as N → +∞ as any other consistent method of simulation.  


The pursuit of differential privacy pits a direct tradeoff against statistical efficiency <cit.>. But the efficiency-privacy tradeoff as a statistical consideration is interweaved with the approximation-exactness tradeoff as a computational consideration, a sentiment that is shared by explorations of other simulation-based Bayesian computational algorithms with differentially private data, including stochastic gradient Monte Carlo <cit.> and Gibbs sampling <cit.>. For  algorithms, to insist on maximal statistical efficiency necessitates computational approximation. Whereas the act of data perturbation not only gains differential privacy, but also computational exactness for free. Both the  and Monte Carlo EM approaches  adapt to differentially private data using the same logic, by setting the tuning parameters governing their numerical performance based on the privacy parameters. Tailoring an algorithm according to the data generative specification exploits the alignment between the statistical and computational tradeoffs, hitting two birds with one stone, so to speak. 



There are several computational challenges to the practical implementation of the proposed frameworks. These challenges are of two types: those intrinsic to  and other forward sampling techniques, and those induced by the privacy mechanism. A weakness in either of these aspects may impact the computational efficiency of these proposals, or in the worst case, render them infeasible. We discuss the two types of challenges below.






A data analyst operating under the dissemination mode of data access is on the receiving end of data products which are designed and privatized by the data curator. As this paper discusses the migration of existing statistical methodology to accommodate privacy-protected data products, we assume that the analyst knows how to perform their preferred analysis on the data product were it not privatized, i.e. if the curator releases  rather than _. That is, π(θ|) in (<ref>) is taken to be the ultimate posterior the analyst targets. 
Depending on the model, however, the analyst may or may not prefer to use  or other forward sampling techniques to draw inference from π(θ|). The strength of  lies in its ability to handle intractable likelihoods, but it presents several limitations. In the construction of the current paper, the intended query function  (and hence the private query _) may be multi-dimensional,  where each dimension is generated in isolation, in conjunction, or sequentially. In particular, we do not preclude the identity function,  () =, in which case the privacy perturbation is performed element-wise on the full dataset for publication, such as may be encountered in the local differential privacy setting. 

Whenever the full data likelihood does not admit a low-dimensional sufficient summary to , the computational efficiency of both proposed algorithms will likely suffer. For classic , the synthetic data matching step (step 3 in Algorithm <ref>) will be computationally wasteful. The  literature explores the use of approximate summary statistics <cit.> to achieve dimension reduction and efficient matching. Unless carefully designed, however, general approximate summary reduction to  will complicate the expression of the privacy kernel η_, and will destroy the “exact” nature of the proposed algorithm. The lack of sufficient reduction challenges the feasibility of other modes of computational for privacy-aware Bayesian inference as well; see e.g. <cit.>.  The question remains with the data curator: in anticipation of a broad range of data analysis needs, how to choose the query  that provides better statistical utility and computational efficiency?


Another limitation of  methods is that their performance depends on the prior and the nature of the state space. As mentioned previously,  must work with proper prior distributions. This  minimal requirement speaks nothing about the algorithm's efficiency.
As the numerical experiment in Section <ref> demonstrates, the acceptance probability of Algorithm <ref> is low when the observed data is in conflict with the prior. The remedy is to devote more sampling resource to areas of the parameter space for which the data exhibit more support. This is a tautology of sorts, since the area we seek is precisely the area with high posterior density, which may be particularly difficult to locate when the parameter space is high dimensional, and when the prior distribution is diffuse (despite being proper).


There are also computational challenges brought forth by the privacy mechanism. Since both proposed algorithms require the transparency of η_(_|·), any act that deprives the analyst's ability to evaluate this quantity also hinders the proposed computational schemes. Two notable causes to diminished transparency of the privacy mechanism are clamping and post-processing. As discussed in Section <ref>, the curator performs clamping when the query has unbounded global sensitivity. While naïve and conservative clamping (such as presented in Section <ref>) requires little additional work from the analyst, carefully designed clamping procedures typically involve the underlying confidential dataset in a nontrivial fashion <cit.>. The resulting privacy mechanism may not be simply captured by an analytically tractable η_.  In addition, the post-processing of differentially private data products may also complicate an otherwise simple expression of η_. Such is the case if the post-processing operation depends nontrivially on aspects of the observed data. For example, the TopDown algorithm imposes invariants on the differentially private noisy measurements via optimization-based post-processing <cit.>. As a result, the output of the algorithm does not permit a straightforward probabilistic description, which threatens its congeniality as a building block in the data processing pipeline <cit.>. From the statistical point of view, a transparent  privacy mechanism is instrumental to the feasibility of conducting exact statistical inference from privacy-protected data <cit.>. To ensure transparency of the privacy mechanism is yet another challenging task that lies with the data curator. 



































 











 






§ ACKNOWLEDGMENT

The author wishes to thank Xiao-Li Meng for inspiring discussions, as well as John Abowd, Gary King, Zhiqiang Tan, and three anonymous reviewers for helpful comments. The author gratefully acknowledges research support by the National Science Foundation (DMS-1916002).




abbrvnat








§ EXAMPLES OF ADDITIVE PERTURBATION DP MECHANISMS



[ϵ-Laplace mechanism; <cit.>]
In (<ref>), let ∼Lap_p(1), the p-dimensional product of independent and identically distributed standard Laplace variables, and h = ϵ^-1Δ_GS(), where

    Δ_GS()=sup_,'{‖()-(')‖: d(,')=1},

is the global sensitivity of , with ‖·‖ denoting the ℓ_1 norm. Then,  is ϵ-differentially private.


[(ϵ, δ)-Laplace mechanism; <cit.>]
In (<ref>), let ∼Lap_p(1), h = ϵ^-1Δ_ξ(, ), and ξ=ϵ{ 4(p+log(2/δ))}^-1, where   

    Δ_ξ(,)=sup_'{ e^-ξ d(,')Δ_LS(,'):'∈}

is the ξ-smooth sensitivity (ξ>0) of  at , and

    Δ_LS(,)=sup_'{‖()-(')‖:d(,')=1}

is the local sensitivity of  at . Then,  is (ϵ, δ)-differentially private. 


[Gaussian mechanism; <cit.>]
In (<ref>), let ∼  N( 0,  I_p) the p-dimensional standard multivariate Normal variable, h = ϵ^-15√(2log(2/δ))Δ_ξ(, ), and ξ=ϵ{ 4(p+log(2/δ))}^-1. Then,  is (ϵ, δ)-differentially private.


The above examples invoke three notions of functional sensitivity (<ref>)-(<ref>), generally denoted as Δ(), to capture the idea that certain choices of  may be more revealing of individual information in  than others. The global sensitivity measures the extent to which  varies between all conceivable pairs of neighboring datasets, whether or not realized in the observed sample. For example, the global sensitivity of the counting query is 1. On the other hand, the local sensitivity of  measures its maximum variability among neighboring datasets to a given observed dataset . The smooth sensitivity strikes a balance between the two, by providing an upper bound on the local sensitivity at  in such a way that the bound does not vary too quickly as a function of . It is crucial that the scale parameter of the additive perturbation mechanism is chosen as a function of both the sensitivity of  as well as the privacy budget, that is, h = h(ϵ,δ,Δ()).




§ PROOF OF THEOREM <REF>




Let I be the indicator of the event that a draw of θ is accepted. The joint distribution of all quantities produced by the ith iteration is π̃(θ,, I) = π_0(θ) π(|θ) π̃(I |), where π̃(I |) is the Bernoulli mass function with proportion parameter cη_(_|). The marginal distribution of an accepted θ sample is

    π̃(θ| I = 1 ) 
    	= ∫π̃(θ,,I=1)/π̃(I=1)d
    	= ∫π_0(θ)π(|θ)cη_(_|)d/∫∫π_0(θ)π(|θ)cη_(_|)d dθ,

which is equal to π(θ|_) as defined in (<ref>). From here, one can see that the overall acceptance probability of Algorithm <ref> is 

    π̃(I = 1) = π(_)/maxη_(·).



Note that under the special case of additive perturbation, the proof of Theorem <ref> parallels Theorem 1 of <cit.>. However, there is an important conceptual difference. In <cit.>, the conditioning query is a query that was observed noiselessly, but construed as if subject to additive error. The -induced posterior of θ therein, while essentially identical to (<ref>), is not the true posterior of θ but that of a “best model input θ̂” given _. With _ being a privatized query, no pretense is necessary in treating it as observed with error, since it indeed was.
 



§ EFFECTIVE SAMPLE SIZE FOR MONTE CARLO EM


In reference to Algorithm <ref>, at the tth iteration, the normalized version of the importance
sampling weights is

    ω̃_i=c_(t)η_(_|_i) = c_(t)ω_i

where c_(t)=1/π(_|θ^(t)) is the reciprocal
of the current approximation to the observed likelihood and is free
of _i. The weighted estimator ∑_i=1^Nω̃_ib(_i) is a consistent estimator of (<ref>) because

    𝔼(b()|_,θ^(t))     =    ∫ b()π(|_,θ^(t))d
        =    ∫ b()π(|θ^(t))η_(_|)/∫π(|θ^(t))η_(_|)dd
        =    ∫ω̃() b()π(|θ^(t))d.

We have that at the expectation of weights for the tth iteration is

    𝔼_|θ^(t)(ω̃) 
    	   =   ∫ c_(t)η_( _|)π(|θ^(t))d
       =   	c_(t)π(_|θ^(t))=1,

where the subscript |θ^(t) signifies the expectation is evaluated with respect to the current approximation to the latent data likelihood, or equivalently, the proposal distribution of the importance sampler. Similarly,

    var_|θ^(t)(ω̃)     =    𝔼_|θ^(t)(ω̃^2)-𝔼_|θ^(t)^2(ω̃)
        =     c_(t)^2𝔼_|θ^(t)(η^2_(_|))-1.

This gives rise to the effective sample size

    ess^(t)(N)     =     N/(1+var_|θ^(t)(ω̃))
       =    Nπ^2(_|θ^(t))𝔼^-1_|θ^(t)(η^2_(_|)).

See also section 2.5.3 of <cit.>. 





§ OBSERVED SCORE AND FISHER INFORMATION FOR MONTE CARLO EM


We have that the observed data log likelihood

    logπ(_|θ)=log∫π(_|)π(|θ)d,

thus the observed score

    ∇_θlogπ(_|θ)     =    ∫π(_|)∇_θπ(|θ)d/∫π(_|)π(|θ)d
        =    ∫π(_|)∇_θπ(|θ)/π(_|)π(|θ)π(_|)π(|θ)d/∫π(_|)π(|θ)d
        =    ∫∇_θπ(|θ)/π(|θ)π(_,|θ)/∫π(_,|θ)dd
        =    ∫∇_θlogπ(|θ)π(|_,θ)d
        =    𝔼( ∇_θlogπ(|θ) |_,θ).

Writing λ_θ() = ∇_θlogπ(|θ), we have that 𝔼( λ_θ() |_,θ^(t)) serves as the tth approximation to the observed score ∇_θlogπ(_|θ^(t)), giving rise to the expression

    𝔼( λ_θ() |_,θ^(t))≈ 
    m∑_i=1^Nω_iλ_θ(_i).



Similarly the Hessian, or the negative of the observed Fisher information matrix, is

    ∇_θ^2logπ(_|θ)	
    	   =   ∫∇_θ^2π(|θ)/π(|θ)π(|_,θ^(t))d
          	- (∇_θlogπ(_|θ))(∇_θlogπ(_|θ))^⊤
       =   𝔼(∇_θ^2logπ(|θ)+∇_θlogπ(|θ)∇_θlogπ(|θ)^⊤|_,θ) 
          	- (∇_θlogπ(_|θ))(∇_θlogπ(_|θ))^⊤.

Substituting again λ_θ() and the expression for the observed score into the above equation, we have that the tth approximation to the observed Fisher information -∇_θ^2logπ(_|θ^(t))  takes the form

    𝔼(-∇_θλ_θ()-λ_θ()λ_θ()^⊤|_,θ^(t))+𝔼(λ_θ()|_,θ^(t))𝔼(λ_θ()|_,θ^(t))^⊤
          ≈  m∑_i=1^Nω_i{ -∇_θλ_θ(_i)-λ_θ(_i)λ_θ(_i)^⊤} + m^2∑_i=1^N∑_j=1^Nω_iω_jλ_θ(_i)λ_θ(_j)^⊤.

See also the appendix of <cit.>.






§ DETAILS OF SECTION <REF>: PRIVATIZED COUNT INFERENCE


For the Bayesian analysis, by the ϵ-Laplace perturbation mechanism, the conditional distribution of 
_ given  is Lap(,ϵ^-1), which
has density ϵ/2exp(-ϵ|_-|). By construction, _ is not an integer with probability one, hence

    π(_|θ)     =   ∫π(|θ)π(_|)d
       ∝     e^-θ{∑_s=0^⌊_⌋θ^/!e^-ϵ_+ϵ+∑_s=⌈_⌉^∞θ^/!e^ϵ_-ϵ} .

Adopting the notations θ_ϵ^+=θ e^ϵ
and θ_ϵ^-=θ e^-ϵ, the first sum within the brackets can be written as

    e^-ϵ_∑_s=0^⌊_⌋(θ_ϵ^+)^/!=e^θ_ϵ^+-ϵ_F_θ_ϵ^+(⌊_⌋)

where F_λ(a) stands for the Pois(λ)
CDF evaluated at a. Similarly, the second sum can be written as

    e^ϵ_∑_s=⌈_⌉^∞(θ_ϵ^-)^/!=e^θ_ϵ^-+ϵ_(1-F_θ_ϵ^-(⌊_⌋)).

Combining the above with the Gamma prior, π_0(θ)∝θ^α-1e^-βθ,
we have that the posterior π(θ|_) takes the form

    π(θ|_)∝θ^α-1e^-(β+1)θ[Γ(⌈_⌉ ,θ_ϵ^+)/Γ(⌈_⌉)e^θ_ϵ^+-ϵ_+γ(⌈_⌉ ,θ_ϵ^-)/Γ(⌈_⌉)e^θ_ϵ^-+ϵ_],

where θ_ϵ^+=θ e^ϵ,
θ_ϵ^-=θ e^-ϵ, ⌈·⌉ is the ceiling function, and Γ(s,x)=∫_x^∞r^s-1e^-rdr is the incomplete Gamma function with Γ(s)=Γ(s,0) and γ(s,x) = Γ(s) - Γ(s,x). 





Figure <ref> displays additional comparisons between the true posterior and the naïve posterior for the same privatized count  under other choices of prior distributions in the Gamma family. Notice that when the observed count appears highly unlikely under a chosen prior (such as  Gamma(2,1) or Gamma(5,1)), a situation known as prior-data conflict <cit.>, the correct posterior heavily discounts the contribution by the privatized observation. The discounting can be seen from the close alignment between the correct posterior (represented by either the solid green density or the black histogram) and the prior (blue dashed density), which in contrast differ drastically from the naïve posterior (red dotted density) in Figure <ref> (a) and (b). The acceptance rate of Algorithm <ref> in these situations are reported in Table <ref>.

For the implementation of the Monte Carlo EM, three stages of iterations were performed with successively more stringent tolerance levels (|θ^(t)-θ^(t-1)|< 10^-3, 10^-4, and 10^-5) and larger Monte Carlo sample size (N = 10^3, 10^5, and 10^7), using the stable maximizer from the last stage as the starting point. This is a crude rule to let N increase, hence the Monte Carlo error decrease, as θ^(t) approaches the true mle. Advanced adaptive techniques, such as the ascent-based modification of <cit.>, can be employed achieve better performance. 




§ DETAILS OF SECTION <REF>: LALONDE DATASET



Let z_i be the indicator variable for whether subject i received
treatment (z_i = 1) or control (z_i = 0), and y_i the earning in 1978 (in
$1k). The full parameter of the model is θ=(τ,μ,σ_t^2,σ_c^2), for which we posit independent priors

    θ∼π_0(τ)×π_0(μ)×π_0(σ_t^2)×π_0(σ_c^2),

where for concreteness, we use π_0(τ) ∼ N(0, 5), π_0(μ) ∼ N(4, 5), π_0(σ_t) ∼ Gamma(2, 0.2) and π_0(σ_c) ∼ Gamma(2, 0.2) for the analysis. The sampling model is

    y_i| z_i,θ∼ N(τ z_i+μ,σ_t^2z_i+σ_c^2(1-z_i)),

where τ is the difference in average earnings between the treatment and control groups. Equivalently stated, treatment group earnings have the distribution N(μ+τ,σ_t^2) and the control group earnings have distribution N(μ,σ_c^2).

The sufficient statistics for θ are the within-group mean and sample variances

    =(y̅_t,y̅_c,s_t^2,s_c^2) =(1/n_t∑_i:z_i=1y_i,1/n_c∑_i:z_i=0y_i,1/n_t-1∑_i:z_i=1(y_i-y̅_t)^2,1/n_c-1∑_i:z_i=0(y_i-y̅_c)^2).

Due to statistical independence of the sample mean and variance of normal random variables, the likelihood can be equivalently represented by the generative model

    y̅_t,y̅_c,s_t^2,s_c^2| z,θ∼ N(μ+τ,σ_t^2/n_t)× N(μ,σ_c^2/n_c)×σ_t^2/n_t-1χ_n_t-1^2×σ_c^2/n_c-1χ_n_c-1^2.


Through a conservative clamping treatment described in Section <ref>, the ϵ-differentially private statistic  _ is obtained via a Laplace mechanism with independent Laplace noise components with bandwidth h^-1 = (1/3, 1/3, 1/6, 1/6) corresponding to . Since the clamping range well exceeds the anticipated range of observable data, we do not perform inferential correction for truncation. Both the original analysis using  (top row of Figure <ref>) and the naïve analysis using _ (bottom row of Figure <ref>) are carried out in RStan, whereas the correct analysis  (middle row of Figure <ref>) is carried out using rejection  of Algorithm <ref>.



