



\subsection{General Reweighting Framework}
\label{sec:reweighting}

Assume for now that we have access to (an approximation of) the propensity score function \(e(d) = \Pr(Z=1\mid D=d)\). We describe a method for approximating \(e(d)\) in the next subsection.

With the knowledge of \(X^1\) and the size of \(X^0\), we can approximate
\begin{equation*}
    \Pr(D=d\mid Z=1)\approx \frac{|\{\mskip-4mu\{ d: d\in X^1\}\mskip-4mu\}|}{|X^1|}
\end{equation*}
and
\begin{equation*}
    \Pr(Z=0)\approx \frac{|X^0|}{|X^0|+|X^1|},\quad \Pr(Z=1)\approx \frac{|X^1|}{|X^0|+|X^1|}.
\end{equation*}
With these probabilities, we can compute \(\Pr(D=d, Z=1) = \Pr(D=d\mid Z=1)\Pr(Z=1)\). Hence,
\begin{align*}
    \Pr(Z=1\mid D=d) &= \frac{Pr(D=d, Z=1)}{\Pr(D=d)}\\
    &= \frac{\Pr(D=d\mid Z=1)\Pr(Z=1)}{\Pr(D=d)},
\end{align*}
and thus
\begin{equation*}
    \Pr(D=d) = \frac{\Pr(D=d\mid Z=1)\Pr(Z=1)}{\Pr(Z=1\mid D=d)}.
\end{equation*}

If instead of estimating \(\Pr(D=d)\) we want to estimate \(\Pr(D=d\mid Z=0)\), we can do this as follows:
\begin{align*}
    &\Pr(D=d\mid Z=0) = \frac{\Pr(D=d,Z=0)}{\Pr(Z=0)}\\
    &= \frac{\Pr(D=d) - \Pr(D=d,Z=1)}{\Pr(Z=0)}\\
    &= \frac{\Pr(D=d) - \Pr(D=d\mid Z=1)\Pr(Z=1)}{\Pr(Z=0)}.
\end{align*}

Note that we only need the propensity scores for the points in \(X^1\), because these are the only points that we have direct access to and thus the only points that we reweight.


\subsection{Propensity Score Computation}
\label{sec:prop_scores}

To apply the reweighting described in the previous subsection, we need to know the propensity scores \(e(d) = \Pr(Z=1\mid D=d)\). For computing them, we rely on the following assumption:
\begin{assumption}
\(D\) can be well approximated by a mixture model with \(K\) components or classes (e.g., clusters or LDA topics), and the participant and non\hyp participant data differ only in the mixture weights.
\end{assumption}
We denote the class membership random variable with values in \(\{1,\dots,K\}\) by \(C\). Each point \(d\) has an associated distribution \(C\mid D=d\) over the classes, given by a probability vector \(\rho_d\). The idea is to compute the distributions \(C\mid Z=0\) and \(C\mid Z=1\) of classes in \(X^0\) and \(X^1\), respectively, and use these distributions to compute the propensity scores for the points in each class. \textit{DiPPS} approximates the propensity scores in a locally differentially private way via a multi\hyp step procedure:
\begin{enumerate}
    \item Train a model on \(X^1\) that learns the \(K\) different classes and assigns to each point \(d\) that it is invoked on a probability distribution over the classes, that is, an approximation of \(\rho_d\). This model might for example be an LDA topic model, and \(\rho_d\) could be the topic vector for document \(d\) that is normalized to a probability distribution.
    \item Send the model to the non\hyp participating users, i.e., to those with data \(X^0\).
    \item Each of the non\hyp participating users uses the model to compute \(\rho_d\) for their data point \(d\). They then use \(\rho_d\) to sample and return one of the \(k\) classes, by invoking the exponential mechanism with \(\rho_d\) as a utility function (we explain the exponential mechanism below).\label{enum:exponential}
    \item Collect the noisy class counts and apply the postprocessing as described later in this subsection to estimate the distribution of \(C\mid Z=0\), i.e., which fraction of the non\hyp participants lies in each class.
    \item Compute a propensity score for each class (how likely is it that \(Z=1\) for a point in a given class) and use these class propensity scores to compute propensity scores for all \(d\in X^1\) based on \(\rho_d\).
\end{enumerate}

We now describe the different steps in detail. We will start with the model that learns the \(K\) classes, continue with the exponential mechanism and the postprocessing of its outputs, and end with the computation of the propensity scores.

\xhdr{Implementation of the class assignment}
For our experiments we train a clustering model to compute class membership probabilities, where the classes in our case are clusters. The model consists of a dimensionality reduction via PCA, followed by a Gaussian mixture model (GMM) with \(K\) components. This model could be replaced by any other model that can learn different classes in a dataset, such as other clustering algorithms.
Note that whenever the classes are learned from the data (as in our case, where we learn the mixture components) instead of being predetermined, we require that all classes that are represented in the participant data are also represented in the non\hyp participant data. I.e., there needs to be at least one non\hyp participant record from each class of the class assignment model.

\xhdr{The exponential mechanism}
The exponential mechanism \cite{mcsherry2007mechanism} used in step~\ref{enum:exponential} is an algorithm whose outputs fulfill (local) differential privacy (see the definition in the introduction). It gets an input dataset \(x\) and is parametrized by a utility function \(u\) that describes the utility of each potential output given an input. It returns a single value \(r\) out of some range \(\mathcal{R}\) with probability proportional to \(\exp(\frac{\varepsilon u(x, r)}{2\Delta u})\), where \(\Delta u\) is the sensitivity of \(u\), that is, how much the value of \(u\) can change at most due to a change by a single record in the input database. In our setting, for a given point \(d\), \(x = \rho_d\) (the dataset consists only of a single record), and \(\mathcal{R} = {1,\dots,K}\), i.e., the different possible classes. The utility function is given by \(u(\rho_d, r) = \rho_{dr}\), meaning the utility of each class is the probability of \(d\) lying in it. Because \(\rho_d\) is a probability vector and thus \(0\leq u(\rho_d, r) \leq 1\) for any \(\rho_d\) and any \(r\), the sensitivity of \(u\) is 1. We hence sample proportional to \(\exp(\frac{\varepsilon u(x, r)}{2})\), which is a weighted softmax. From the resulting histogram we would like to estimate the distribution of \(C\mid Z=0\).

\xhdr{Postprocessing of the class counts}
If the classes sent by the users were directly sampled from \(\rho_d\), we could simply count how often each of the values \(1,\dots,K\) has been sent to get an estimate of \(C\mid Z=0\). However, since we do sampling via the exponential mechanism, we need to revert the distortion introduced by the weighted softmax function. Let \(p\) be the probability vector describing the probability distribution of \(C\mid Z=0\) and let \(n_0=|X^0|\). Via the exponential mechanism in step~\ref{enum:exponential} we get \(n_0\) samples \(A_1,\dots,A_{n_0}\) with values in \(1,\dots,K\). By randomly shuffling them we can ensure that they are i.i.d.\ with probabilities
\begin{equation*}
    \Pr(A_i = k) = \frac{e^{\varepsilon p_k/2}}{\sum_{l=1}^K e^{\varepsilon p_l/2}}.
\end{equation*}
Let \(\tilde{U}_k = |\{i:A_i=k\}|\), i.e., a random variable counting how often class \(k\) has occurred. Define \(U_k\) for \(k,l=1,\dots,K\) via
\begin{equation}
    \label{eq:Uk}
    U_k - U_l = \frac{2}{\varepsilon}\log\left(\frac{\tilde{U}_k}{\tilde{U}_l}\right)
\end{equation}
and \(\sum_{k=1}^K U_k = 1\). We can thus compute \(U_1 = \frac{1}{K} (1 + \sum_{k=2}^K (U_1 - U_k))\) and from there \(U_2,\dots,U_K\).
Because \((\tilde{U}_1,\dots,\tilde{U}_K)\) follows a multinomial distribution, a zeroth\hyp order Taylor approximation yields \(\tilde{U}_k/\tilde{U}_l \approx \mathbb{E}(\tilde{U}_k) / \mathbb{E}(\tilde{U}_l) = \exp\left(\frac{\varepsilon(p_k-p_l)}{2}\right)\) \cite{duris2018mean} and hence \(U_k - U_l \approx p_k - p_l\). Therefore, \(U_k\) is an approximation of \(\Pr(C=k\mid Z=0)\).

\xhdr{Approximating the propensity scores}
We next define the \emph{cluster propensity score} \(\tilde{e}\), which can be interpreted as the propensity score of points that lie in exactly one cluster \(k\), i.e., whose cluster probability distribution is binary:
\begin{align*}
    \tilde{e}(k) &= \Pr(Z=1\mid C=k)\\
    &=\frac{\Pr(C=k\mid Z=1)\Pr(Z=1)}{\Pr(C=k)}\\
    &=\frac{\sum_d \Pr(D=d, C=k\mid Z=1)\Pr(Z=1)}{\sum_{d}\Pr(D=d, C=k)}\\
    &\approx\frac{\frac{\sum_{d\in X^1} \rho_d(k)}{|X^1|} \frac{|X^1|}{|X^0| + |X^1|}}{\frac{\sum_{d\in X^1} \rho_d(k) + U_k}{|X^0| + |X^1|}}\\
    &= \frac{\sum_{d\in X^1} \rho_d(k)}{\sum_{d\in X^1} \rho_d(k) + U_k}.
\end{align*}
With the cluster propensity scores, we are finally able to compute an approximation of the propensity score \(e(d)\) as
\begin{align*}
    e(d) &= \Pr(Z=1\mid D=d)\\
    &=\sum_{k=1}^K \Pr(Z=1\mid C=k,D=d)\Pr(C=k\mid D=d)\\
    &=\sum_{k=1}^K \Pr(Z=1\mid C=k)\Pr(C=k\mid D=d)\\
    &= \sum_{k=1}^K \tilde{e}(k) \rho_d(k).
\end{align*}

\begin{figure*}[t!]
		\centering
        \includegraphics[width=0.8\textwidth]{figures/DiPPS_diagram.pdf} 
  \caption{Overview of \textit{DiPPS}.}
        \label{fig:diagram}
\end{figure*}

\Figref{fig:diagram} shows an overview of our method.
\subsection{Tasks}
The goal of \textit{DiPPS} is to approximate a probability distribution. Hence it is most natural to evaluate it w.r.t.\ a metric for the distance between distributions. As the metric we choose the Wasserstein distance of order 1. It measures how expensive it is to transform one distribution into the other, when the cost of transporting a given amount of probability mass between two points is the transported amount times the Euclidean distance between the two points. As compared to, e.g., the Kullback\hyp Leibler divergence, the Wasserstein distance does not require the two distributions to have the same support, which is important for us, because we work with finite samples from, in most cases, continuous distributions. We compute the Wasserstein distance using the R transport package \cite{transport2020}. Because the computation
--- which is only required for the evaluation, not for the deployment of our method ---
takes quite long, we only compute the distance for \(\varepsilon=1\).

We also evaluate \textit{DiPPS} on several potential downstream tasks. These tasks are the computation of mean, variance and median of each attribute for the entire dataset and for the non\hyp participant dataset. We measure performance in terms of the absolute value of the deviation from the ground truth, divided by the number of features in the dataset, to make different datasets comparable.

We repeat each experiment five times and report the means of the five repetitions. In addition, for \(\varepsilon=1\) we report the standard deviation of the Wasserstein distance (\Tabref{table:was}) and of the error when estimating mean, variance and median (Tables~\ref{table:mean_err}--\subref{table:median_err}).

\subsection{Datasets}
Throughout this subsection, we evaluate \textit{DiPPS} and its competitors on several datasets. We normalize all features to lie in the interval \([-1,1]\) to make the experiment results comparable. Each dataset is split according to one of its variables into participant and non\hyp participant data.

\xhdr{Web visits} \cite{Dua:2020}
This dataset contains web traces of different users of msnbc.com within a 24h period. The website visits are recorded at the level of 17 URL categories such as ``tech'' or ``business''. We convert this dataset into one that contains for each user whether they have visited a URL of a category zero times, one time or more than once. We cannot use the original unbounded counts, because the \textit{Laplace} and the \textit{Hybrid} mechanism that we compare with require attributes in a bounded range.\\
\textbf{Splitting}: We split the dataset on the ``bulletin board service'' (bbs) attribute. Users who used a bbs at least once are treated as participants, users who never used it are treated as non\hyp participants. The setting could be that of a provider of an Internet platform that is able to record the web visits of their users within the platform but no web visits outside of the platform.
The dataset contains 2k records for which the bbs value is non\hyp zero. From the remaining 988k records, we sample 2k at random as the non\hyp participants.

\xhdr{Credit cards} \cite{yeh2009comparisons}
The dataset contains information about credit card customers. This includes demographic information such as age and gender, and information about their payment history, that is, their repayment status, their previous payments and their bill statements in previous months. In total there are 24 attributes.\\
\textbf{Splitting}: We split the dataset according to the binary variable whether someone defaulted on their credit card payments or not. Participants are those who did not default, non\hyp participants those who did. The reasoning is that people who are struggling financially might be less willing to disclose information about their finances. There are 23.5k participants and 6.5k non\hyp participants.

\xhdr{Food} \cite{foodaps2016}
This data stems from the National Household Food Acquisition and Purchase Survey of the U.S.\ Department of Agriculture, in which data about food purchases of households was collected. We choose a subset of four attributes that are the answers to questions about whether respondents could afford enough food and the food they wanted in the past 30 days.\\
\textbf{Splitting}: For splitting we choose a fifth variable that indicates whether anyone in the households receives benefits from the Supplemental Nutrition Assistance Program, a governmental program to financially support food purchases. Those who receive benefits are participants, those who do not receive benefits are non\hyp participants. It could, e.g., be that people who are supported by a nutrition assistance program have to provide data about their food purchases to receive the benefits, while others do not. There are 1.5k participants and 3k non\hyp participants in the dataset.

\xhdr{Weather} \cite{zhang2017cautionary}
The dataset contains hourly air pollutant and weather data from twelve air\hyp quality monitoring sites in the Beijing area over a span of four years. The air pollutant attributes include PM2.5 and NO2 concentration, the weather attributes include temperature and precipitation. In total there are twelve attributes.\\
\textbf{Splitting}: We use Wanliu, a site in the city center of Beijing and Dingling, a site in a more rural area nearby. We do two kinds of experiments for this dataset: Once the city data is treated as the participant data and the rural data as the non\hyp participant data (``Weather 1''), and once it is the other way around (``Weather 1 Inverse''). One could imagine a scenario where an organization has a weather station in one region of a country, but would like to also gather weather information from a different region. The data collected from the weather station would be the participant data. Individuals from the other region of interest could donate weather data that they have collected in a differentially private way; this would then be the non\hyp participant data. To show that our results are robust w.r.t.\ the chosen sites, we repeat the experiments with the sites Nongzhanguan (city center) and Nongzhanguan (rural area); we denote those experiments with ``Weather 2'' and ``Weather 2 Inverse''. For each site there are between 30.5k and 33k records after the removal of missing data.


\subsection{Comparison Methods}
\label{sec:baselines}
To the best of our knowledge, \textit{DiPPS} is the first method for differentially private distribution estimation in the data setting that we consider. Therefore, there does not exist any method that we can directly compare with (apart from naive estimation, see below). Instead, we compare with a selection of more specialized methods that are limited to computing a small set of functions. As opposed to that, with \textit{DiPPS} any number of arbitrary functions of the data distribution can be computed. This is just to give an idea of the performance of \textit{DiPPS}; we do \emph{not} aim at improving upon the state\hyp of\hyp the\hyp art in these specialized tasks.

\xhdr{Naive}
Compute all functions on the participant data and ignore the non\hyp participant data. This measures the amount of participation bias in the data.

\xhdr{Propensity scores (\textit{PS})}
A strong ceiling: our method, but without the exponential mechanism, i.e., without a differential privacy guarantee.
Instead of sampling via the exponential mechanism, the non\hyp participants directly sample from the cluster distribution.

\xhdr{Laplace mechanism} \cite{dwork2006calibrating}
The \textit{Laplace} mechanism adds noise from a Laplace distribution to each non\hyp participant record and shares these noisy records with the data collector. We merge the resulting noisy dataset with the non\hyp noisy dataset of the participants. The \textit{Laplace} mechanism can be used to estimate mean and median.

\xhdr{Hybrid mechanism} \cite{wang2019collecting}
The \textit{Hybrid} mechanism is a recently introduced, state\hyp of\hyp the\hyp art mechanism for locally differentially private mean estimation. It works similar to the \textit{Laplace} mechanism, but instead of adding unbounded noise to the data, it returns a random value from a bounded interval based on its input. It combines the piecewise mechanism introduced in the same paper with the earlier Duchi mechanism \cite{duchi2018minimax}, and improves upon them in terms of worst\hyp case noise variance. While theoretically it could also be used to estimate the median, the \textit{Hybrid} mechanism is designed for estimating the mean and indeed performs much worse when estimating the median than even the naive method, and hence we only include it for comparison in the estimation of the mean.

\subsection{Hyperparameters of \textit{DiPPS}}
Unfortunately, we cannot choose the best hyperparameters for the downstream task, because in the real\hyp world setting we do not have access to the data of the non\hyp participants and thus no ground truth. We hence resort to choosing the number of PCA dimensions such that at least 80\% of the variance is retained, and the number of GMM components based on the GMM log\hyp likelihood using the elbow method \cite{thorndike1953belongs}.

\subsection{Results}
\input{tables/wasserstein_dist}
\xhdr{Wasserstein distance}
In \Tabref{table:was} we show the difference in Wasserstein distance with Euclidean norm cost between the non\hyp participant distribution estimated by \textit{PS} and \textit{DiPPS} (for \(\varepsilon=1\)) and true non\hyp participant distribution, and compare this with the distance between the participant distribution and the non\hyp participant distribution (\textit{Naive}). Note that we do not show the results for the distance between the estimated and true entire distribution, because they can be computed from the given results by multiplying the distance by \(|X^0|/(|X^0| + |X^1|)\). As a consequence, whenever one method outperforms the other for estimating the non\hyp participant distribution, it will also outperform that method for estimating the entire distribution.
We can see from the table that \textit{DiPPS} outperforms \textit{Naive} on all datasets, which shows that our method indeed successfully reduces participation bias as measured in Wasserstein distance. The Wasserstein distance is widely used, because it captures differences between distributions well. Hence a reduction in Wasserstein distance as achieved by \textit{DiPPS} will most likely be correlated with a reduction in participation bias in many downstream tasks. We confirm this hypothesis in the next subsection. To show the cost of privacy, we also include the non\hyp private \textit{PS} in the table.


\begin{figure*}[!htbp]
\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=0.998\textwidth]{figures/results_nonopt_complete.pdf}
  \subcaption{Error for non\hyp participant data.}
        \label{fig:stats_nonopt}
\end{subfigure}

\vspace*{1em}

\begin{subfigure}{\textwidth}
		\centering
        \includegraphics[width=0.998\textwidth]{figures/results_all_complete.pdf}
  \subcaption{Error for entire data.} 
        \label{fig:stats_all}
\end{subfigure}
\caption{\boldmath Mean absolute error (\(=2\ \times\) relative error) per attribute when estimating mean, variance and median of the non\hyp participant (\subref{fig:stats_nonopt}) and entire (\subref{fig:stats_all}) data of different datasets for different $\varepsilon$; lower is better. Averaged over five runs.\unboldmath}
\end{figure*}

\input{tables/error_all}

\xhdr{Statistics}
\Figref{fig:stats_nonopt} and \ref{fig:stats_all} show the results for estimating mean, variance and median of the non\hyp participant dataset and the entire dataset, respectively. Tables~\ref{table:mean_err}--\subref{table:median_err} contain the numerical values for estimating the non\hyp participant statistics at \(\varepsilon=1\). \textit{DiPPS} outperforms the naive estimation for \(\varepsilon \geq 1\) in almost all settings. This confirms our hypothesis that a smaller Wasserstein distance corresponds to a smaller error on downstream tasks. In most settings, \textit{DiPPS} outperforms the \textit{Laplace} mechanism. And even the \textit{Hybrid} mechanism is beaten by \textit{DiPPS} on two out of the five datasets.

We want to emphasize that our main goal is not to improve upon existing locally differentially private mechanisms that are specialized for certain tasks. While those are only suitable for computing one specific function, when using our method \emph{any} function of the dataset can be computed --- without decreasing the privacy guarantee when increasing the number of function computations ---, because \textit{DiPPS} estimates a data distribution instead of a specific function. Hence \textit{DiPPS} fulfills a much more holistic purpose than specialized mechanisms.

\xhdr{Influence of \boldmath\(\varepsilon\)\unboldmath}
In all plots it can be seen that the errors of the \textit{Laplace} and the \textit{Hybrid} mechanism decrease with increasing \(\varepsilon\). This is expected because the added noise gets smaller the larger \(\varepsilon\) is. However, \textit{DiPPS} does not necessarily behave the same; see, e.g., the variance estimation for the Weather 1 dataset. The reason is that the noise due to the exponential mechanism is more subtle. \(\varepsilon\) influences the sampling process, and the optimal \(\varepsilon\) depends on the data. For example, a very large \(\varepsilon\) would mean that in all but very few cases the class that has the highest probability for a point would get sampled, ignoring the probabilities of the other classes. If all data points have probability close to \(1\) for a single class, this would be desirable. However, if there is, e.g., one class that has probability considerably larger than \(0\) for many points, but for none of the points is the class with the highest probability, this class might not get sampled at all and thus its frequency would get underestimated.
The subtle influence of the exponential mechanism is also the reason why in some cases \textit{DiPPS} outperforms \textit{PS}. The distribution from which our datasets were sampled is most likely not an exact Gaussian mixture model. This means that the approximations of the propensity scores are not perfect. In some cases, the distortions introduced by the exponential mechanism correct part of this error and lead to slightly better results than when sampling directly from the cluster distribution.


\section{Introduction}
\label{sec:intro}
\input{010_intro}

\section{Related Work}
\label{sec:related}
\input{020_related}

\section{Problem Definition}
\label{sec:problem}
\input{030_problem}

\section{Proposed Solution}
\label{sec:solution}
\input{040_solution}

\section{Experiments}
\label{sec:experiments}
\input{050_experiments}

\section{Discussion}
\label{sec:discussion}
\input{060_discussion}

\section{Conclusion}
\label{sec:conclusion}
\input{070_conclusion}



