text:
\begin{document}
\thispagestyle{empty}
\title{One Step to Efficient Synthetic Data}

\author{\name Jordan Awan \email jawan@purdue.edu \\
       \addr Department of Statistics\\
       Purdue University\\
       West Lafayette, IN 47907, USA
       \AND
       \name Zhanrui Cai \email zhanruic@andrew.cmu.edu \\
       \addr Department of Statistics and Data Science\\
       Carnegie Mellon University \\
	    Pittsburgh, PA 15213,  USA}

\editor{TBD}

\maketitle

	\begin{abstract}
        A common approach to synthetic data is to sample from a fitted model. We show that under general assumptions, this approach results in a sample with inefficient estimators and whose joint distribution is inconsistent with the true distribution. Motivated by this, we propose a general method of producing synthetic data, which is widely applicable for parametric models,  has asymptotically efficient summary statistics, and is both easily implemented and highly computationally efficient.
		Our approach allows for the construction of both partially synthetic datasets, which preserve certain summary statistics, as well as fully synthetic data which satisfy the strong guarantee of differential privacy (DP), both with the same asymptotic guarantees. We also provide theoretical and empirical evidence that the distribution from our procedure converges to the true distribution. Besides our focus on synthetic data, our procedure can also be used to perform approximate hypothesis tests in the presence of intractable likelihood functions.
	\end{abstract}

\begin{keywords}
  Synthetic data, differential privacy, indirect inference, parametric bootstrap, statistical disclosure control
\end{keywords}

	\section{Introduction}\label{s:intro}

	With the advances in modern technology, government and other research agencies are able to collect massive amounts of data from individual respondents. These data are valuable for scientific progress and policy research, but they also come with increased privacy risk \citep{lane2014privacy}. To publish useful information while preserving confidentiality of sensitive information, numerous methods of generating \emph{synthetic data} have been proposed (see \citet[Chapter 3]{hundepool2012statistical} for a survey). The goal of synthetic data is to produce a new dataset which preserves distributional properties of the original dataset, while protecting the privacy of the participating individuals. There are two main types of synthetic data: \emph{partially synthetic data}, which allows for certain statistics or attributes to be released without privacy while protecting the other aspects of the data, and \emph{fully synthetic data}, where all statistics and attributes of the data are protected.

	There has also been an increased interest in developing synthetic data with formal privacy guarantees, such as \emph{differential privacy} (DP). 	Differential privacy (DP) was proposed in \citet{dwork2006calibrating} as a framework to develop formally private methods. Methods which satisfy DP require the introduction of additional randomness, beyond sampling in order to obscure the effect of one individual on the output. Intuitively, DP ensures plausible deniability for those participating in the dataset.

		A common approach to synthetic data is that of	\citet{liew1985data}, which proposes to
		draw synthetic data from a fitted model.

		While \citet{liew1985data} do not incorporate formal privacy methods, this approach is often used to produce differentially private synthetic data.
	\citet{hall2013differential} develop DP tools for kernel density estimators, which can be sampled to produce DP synthetic data.

	\citet{machanavajjhala2008privacy} develop a synthetic data method based on a multinomial model, which satisfies a modified version of DP to accomodate sparse spatial data.

	\citet{mcclure2012differential} sample from the posterior predictive distribution to produce DP synthetic data, which is asymptotically similar to the \citet{liew1985data}.

	\citet{liu2016model} also use a Bayesian framework: first they produce DP estimates of the Bayesian sufficient statistics, draw the parameter from the distribution conditional on the DP statistics, and finally sample synthetic data conditional on the sampled parameter.

	\citet{zhang2017privbayes} propose a method of developing high-dimensional DP synthetic data which draws from a fitted model based on differentially private marginals.

	Another approach to synthetic data samples from a \emph{conditional distribution}, preserving certain statistics. The most fundamental perspective of this approach is that of \citet{muralidhar2003theoretical}, who propose drawing confidential variables from the distribution conditional on the non-confidential variables.

	\citet{burridge2003information} generate partially synthetic data, preserving the mean and covariance for normally distributed variables.
	This approach was extended to a computationally efficient version by \citet{mateo2004fast}, and  \citet{ting2005romm} give an alternative approach to preserving the mean and covariance by using random orthogonal matrix multiplication.

	There are also tools, often based in algebraic statistics, to sample conditional distributions preserving certain statistics for contingency tables. \citet{karwa2013conditional} give a survey of Markov Chain Monte Carlo (MCMC) techniques to sample conditional distributions. Another approach is sequential importance sampling, proposed by \citet{chen2006sequential}.
	\citet{slavkovic2010synthetic} use these techniques to generate synthetic contingency tables that preserve conditional frequencies.

	While not focused on the problem of synthetic data, there are other notable works sampling from conditional distributions. A series of works sample from distributions conditional on sufficient statistics focused on the application of hypothesis testing \citep{lindqvist2005monte,lindqvist2007conditional,lindqvist2013exact,lillegard1999exact,engen1997stochastic,lillegard1999exact,taraldsen2018conditional}. \citet{barber2020testing} recently gave a method of sampling conditional on an efficient statistic using the principle of asymptotic sufficiency; they showed that their method results in asymptotically valid $p$-values for certain hypothesis tests.

	In differential privacy, there are also synthetic data methods which preserve sample statistics.
	\citet{karwa2012differentially} generate DP synthetic networks from the beta exponential random graph model, conditional on the degree sequence.

	\citet{li2018privacy} produce DP high dimensional synthetic contingency tables using a modified Gibbs sampler.

	\citet{hardt2012simple} give a distribution-free algorithm to produce a DP synthetic dataset, which approximately preserves several linear statistics.

	While this paper will be focused on producing synthetic data for parametric models, there are several non-parametric methods of producing synthetic data, using tools such as multiple imputation \citep{rubin1993statistical,raghunathan2003multiple,drechsler2011synthetic}, regression trees \citep{reiter2005using,drechsler2008accounting}, and random forests \citep{caiola2010random}. Recently there has been success in producing differentially privacy synthetic data using generative adversarial neural networks \citep{jordon2018pate,triastcyn2018generating,xu2019ganobfuscator}.

	{\bf Our contributions and organization }

	The related work discussed above on synthetic data largely fits into one of two categories: 1) sampling from a fitted distribution or 2) sampling from a distribution conditional on sample statistics. Our first result, Theorem \ref{thm:TVnaive} shows that in very general settings, the first approach results in a sample with inefficient estimators, and whose distribution is ``inconsistent.'' In particular, we show that the joint distribution of the synthetic sample does not converge in total variation to the true joint distribution as the sample size increases. This result gives a strong indication that sampling from a fitted model is not ideal for synthetic data generation, whereas sampling conditional on certain sample statistics can overcome these issues.

	However, there are important limitations to the previous works which sample from a conditional distribution. First, the previous approaches are all highly specific to the model at hand, and require different techniques for different models. Second, many of the approaches are difficult to implement and computationally expensive, involving complex iterative sampling schemes such as MCMC.

The approach we propose in this paper also preserves summary statistics, but unlike previous methods it is applicable to a wide variety of parametric models, easily implemented, and highly computationally efficient. 	Our approach allows for the construction of both partially synthetic datasets, which preserve the summary statistics without formal privacy methods, as well as fully synthetic data which satisfy the strong guarantee of differential privacy (DP).

    More formally, given sensitive data $X_1,\ldots, X_n \iid f_\theta$ with unknown parameter $\theta$ and efficient estimator $\hat\theta$, we propose a method of producing a synthetic dataset $(Y_i)_{i=1}^n$, which is conditionally independent of $(X_i)_{i=1}^n$ given $\thx$ and satisfies $\thy= \thx + o_p(n^{-1/2})$.  We also provide both theoretical and empirical evidence suggesting that the joint distribution which results from our method converges to the true joint distribution, resulting in high quality synthetic data. In occasions where $\thx$ itself must be protected, we show that the proposed synthetic mechanism can be easily modified to satisfy the strong guarantee of differential privacy by first privatizing $\thx$.
	Our approach can be viewed as one step of an approximate Newton method, which aims to solve $\thy=\thx$.
	 Our approach is easily applied as the regularity conditions required for our asymptotics are similar to those required for the Central Limit Theorem of the maximum likelihood estimator (MLE), the computations only require efficient estimators for the parameters and the ability to sample the model, and the computational time is proportional to simply fitting the model.

	The rest of the paper is organized as follows: In Section \ref{s:background}, we review some terminology and set the notation for the paper. In Section \ref{s:bootstrap}, we prove the limitations of synthetic data generated by sampling from a fitted model, showing that the distribution is inconsistent and has inefficient summary statistics. We propose our one-step approach to synthetic data in Section \ref{s:asymptotics} and prove that it results in a sample which preserves an efficient estimator with asymptotically negligible error. In Section \ref{s:distribution}, we consider the distribution of the one-step synthetic sample, and prove that a related procedure results in a consistent sample, giving evidence that the one-step procedure itself is consistent. In Section \ref{s:simulations} we perform several simulation studies illustrating 1) the efficient estimators of the one-step synthetic data, 2) the distributional properties of the approach, 3) that the approach can give high quality DP synthetic, 4) the the one-step synthetic data can perform accurate and powerful hypothesis tests on models with intractible likelihood functions, and 5) the improved running time of our approach compared to a competing MCMC approach to synthetic data.  We end in Section \ref{s:conclusions} with some discussion. Proofs and technical details are postponed to Appendix \ref{s:proofs}.

	\section{Background and notation}\label{s:background}
	In this section, we review some background and notation that we use throughout the paper.

	For a parametric random variable, we write $X\sim f_\theta$ to indicate that $X$ has probability density function (pdf) $f_\theta$.

To indicate that a sequence of random variables from the model $f_\theta$ are independent and identically distributed (i.i.d.), we  write $X_1,\ldots, X_n \iid f_\theta$. We write $\ul X =(X_i)_{i=1}^n= (X_1,\ldots, X_n)^{\top}$.

	For $X\sim f_\theta$, we denote the \emph{score function} as $S(\theta,x) = \nabla_\theta \log f_\theta(x)$, and the \emph{Fisher Information} as $I(\theta) = \EE_\theta\left[S(\theta,X) S^\top(\theta,X)\right]$.  An estimator $\hat \theta: \mscr X^n \rightarrow \Theta$ is \emph{efficient} if for $X_1,\ldots, X_n \iid f_\theta$, we have $\sqrt n (\hat \theta(\ul X) - \theta) \overset d \rightarrow N(0, I^{-1}(\theta))$. We will often write $\thx$ in place of $\hat \theta(\ul X)$.

\section{Limitations of the parametric bootstrap for synthetic data}\label{s:bootstrap}
Sampling from a fitted model, also known as the parametric bootstrap, is a common approach to synthetic data.

However, the parametric bootstrap has considerable weaknesses when used in this setting, in that it results in significantly worse approximations to the true sampling distribution. In this section, we prove that the parametric bootstrap gives inefficient sample statistics and results in ``inconsistent'' synthetic data, where we show that the total variation distance between the true distribution and the parametric bootstrap approximation does not go to zero as $n\rightarrow \infty$. We also recall other limitations of the parametric bootstrap for hypothesis testing.

From the perspective of synthetic data, the goal is to produce a new dataset $\ul Y$ which is approximately equal in distribution to $\ul X$, where the approximation is measured by some metric or divergence. At minimum, we may expect that the distribution of $\ul Y$ approaches the distribution of $\ul X$ as the sample size $n$ grows. We begin with an example that show that the parametric bootstrap results in suboptimal asymptotics, calling into question the appropriateness of the parametric bootstrap for the generation of synthetic.

	\begin{ex}\label{s:normal}
		Suppose that $X_1\ldots, X_n \iid N(\mu,1)$. We use the esimator $\hat \mu(\ul X)= n^{-1} \sum_{i=1}^n X_i$ and draw $Z_1,\ldots, Z_n|\hat\mu(\ul X) \iid N\{\hat \mu(\ul X),1\}$.
		We can compute $\var(\hat \mu(\ul X)) = n^{-1} $, whereas $\var\{\hat \mu(\ul Z)\} = 2n^{-1}$.
		By using the synthetic data $\ul Z$, we have lost half of the effective sample size. We can also derive $(Z_1,\ldots, Z_n)^\top \sim N(\mu \ul 1, I_n + n^{-1} \ul 1\, \ul 1^\top)$, where $\ul 1 = (1,\ldots, 1)^\top$ is a vector of length $n$. It is easily calculated that the Hellinger distance between the distributions of $\ul X$ and $\ul Z$ is $h(\ul X,\ul Z) = \sqrt{1- 2^{3/4}(3)^{-1/2}} \approx .17$, regardless of the sample size $n$\footnote{Use the formula for the Hellinger distance between Gaussian distributions and Sylvester's Determinant Theorem \citep[p271]{pozrikidis2014introduction}.}. This demonstrates that the samples $\ul X$ and $\ul Z$ are fundamentally different, and that their distributions do not approach each other as $n$ grows.
	\end{ex}

	 Example \ref{s:normal} is in fact one instance of a very general phenomenon. In Theorem \ref{thm:TVnaive}, we show that when $\hat \theta(\cdot)$ is an efficient estimator, then $\thz$ is an inefficient estimator for $\theta$, and the distribution of $\ul Z$ is ``inconsistent'' in that the distributions of $\ul Z$ and $\ul X$ do not converge in total variation, as $n\rightarrow \infty$.

\begin{restatable}{thm}{thmTVnaive}\label{thm:TVnaive}
Let $X_1,\ldots, X_n\iid f_\theta$, such that there exists an efficient estimator $\hat\theta(\ul X)$ for $\theta$. Sample $Z_1,\ldots, Z_n |\hat\theta(\ul X) \iid f_{\hat\theta(\ul X)}$. Then
\begin{enumerate}
    \item $\sqrt n \{\hat\theta(\ul Z)-\theta\} \overset d \rightarrow N\{0, 2I^{-1}(\theta)\}$, whereas $\sqrt n \{\hat\theta(\ul X)-\theta\} \overset d \rightarrow N\{0, I^{-1}(\theta)\}$, and
    \item $\mathrm{TV}\{(X_1,\ldots, X_n),(Z_1,\ldots, Z_n)\}$ does not converge to zero as $n\rightarrow \infty$.
\end{enumerate}
\end{restatable}
The use of the total variation distance in Theorem \ref{thm:TVnaive} also has a hypothesis testing interpretation: if the parameter $\theta$ were known, then given either the sample $(X_1,\ldots, X_n)$ or $(Z_1,\ldots, Z_n)$, we can always construct a test to discern between the two distributions with power greater than its type I error. So, the samples $\ul X$ and $\ul Z$ never become indistinguishable. In short, Theorem \ref{thm:TVnaive} shows that the parametric bootstrap is not ideal for the generation of synthetic data.

    The parametric bootstrap also has limitations when applied to hypothesis testing problems. Indeed, \citet{robins2000asymptotic} show that if the test statistic $T(X)$ is asymptotically normal with mean $a(\theta)$, then the parametric bootstrap $p$-values are asymptotically conservative if $a(\theta)$ is not a constant in $\theta$; this results in tests with low power. On the other hand, if $a(\theta)$ is constant in $\theta$, then the $p$-values are asymptotically uniform under the null.

For complex models, it may be difficult to construct a test statistic with the necessary properties to ensure that the $p$-values are appropriately calibrated. On the other hand, \citet{barber2020testing} showed that using approximate co-sufficient sampling methods, similar to those in this paper, guarantee asymptotically uniform $p$-values under the null hypothesis, no matter which test statistic is used. In Section \ref{s:distribution}, we give evidence that our method also results in approximate co-sufficient samples, and demonstrate in Section \ref{s:DPtesting} that samples from our method have more accurate $p$-values and higher power compared to the parametric bootstrap.

	\section{One-step solution to synthetic data}\label{s:asymptotics}

    In this section, we present our synthetic data procedure and its sample statistic asymptotics in Theorem \ref{thm:onestep}. We also include a pseudo-code version of our approach in Algorithm \ref{alg:onestep}, to aid implementation.

    While sampling from a fitted model is a common approach to synthetic data, we saw in Theorem \ref{thm:TVnaive} that it results in inferior asymptotics of both the sample estimates as well as the joint distribution of the synthetic data. Our approach avoids the asymptotic problem of Theorem \ref{thm:TVnaive} by producing a sample $(Y_i)_{i=1}^n$ such that $\thy = \thx + o_p(n^{-1/2})$. Then marginally, the asymptotic distributions of $\thy$ and $\thx$ are identical. The intuition behind the method is that after fixing the ``seed,'' we search for a parameter $\theta^*$ such that when $(Y_i)_{i=1}^n$ are sampled from $f_{\theta^*}$, we have that $\thy =\thx+o_p(n^{-1/2})$.  	To arrive at the value $\theta^*$, we use one step of an approximate Newton method, described in Theorem \ref{thm:onestep}. Our method can also be viewed an one step of the iterative bootstrap, with a batch size of 1 \citep{guerrier2019simulation}.

    To facilitate the asymptotic analysis, we assume regularity conditions (R0)-(R3). (R1)-(R3) are similar to standard conditions to ensure that there exists an efficient estimator, which are relatively mild and widely assumed in the literature \citep{Serfling:1980,lehmann2004elements}.
    Assumption (R0) makes the notion of ``seeds'' precise in terms of measure theory, and also ensures that the transformation $X_\theta$ varies smoothly in $\theta$.

	\begin{itemize}
	\item   [(R0)] Let $(\Omega, \mscr F, P)$ be a probability space of the \emph{seed} $\omega$. Let $X_\theta:\Omega \rightarrow \mscr X$ be a measurable function which is continuously differentiable in $\theta$, where $(\mscr X,\mscr G)$ is a measurable space and $\theta \in \Theta \subset \RR^p$, where $\Theta$ is compact. We assume that there exists a measure $\mu$ on $(\mscr X,\mscr G)$ which dominates $PX^{-1}_{\theta}$ for all $\theta \in \Theta$. Then there exist  densities $f_\theta:\mscr X \rightarrow \RR^{\geq 0}$ such that $\int_A dPX^{-1}_\theta = \int_A f_\theta d\mu$ for all $A\in \mscr G$.
		\item [(R1)]  Let $\theta_0\in \Theta\subset \RR^p$ be the true parameter. Assume there exists an open ball $B(\theta_0)\subset \Theta$ about $\theta_0$, the model $f_\theta$ is identifiable, and that the set $\{x\in \mscr X\mid f_\theta(x)>0\}$ does not depend on $\theta$.
		\item [(R2)] The pdf $f_\theta(x)$ has three derivatives in $\theta$ for all $x$ and
		there exist functions $g_i(x)$, $g_{ij}(x)$, $g_{ijk}(x)$ for $i,j,k=1,\ldots, p$  such that for all $x$ and all $\theta \in B(\theta_0)$,
		\[\left| \frac{\partial f_{\theta}(x)}{\partial \theta_i}\right|\leq g_i(x),\qquad
		\left| \frac{\partial^2 f_{\theta}(x)}{\partial \theta_i\partial \theta_j}\right|\leq g_{ij}(x),\qquad
		\left| \frac{\partial^3 f_{\theta}(x)}{\partial \theta_i\partial \theta_j\partial \theta_k}\right|\leq g_{ijk}(x).\]
		We further assume that each $g$ satisfies $\int g(x) \ dx<\infty$ and $\EE_\theta g_{ijk}(X)<\infty$ for $\theta \in B(\theta_0)$.
		\item [(R3)] The Fisher Information matrix $I(\theta_0) = \EE_{\theta_0} [(\nabla_{\theta_0} \log f_{\theta_0}(X))(\nabla_{\theta_0} \log f_{\theta_0}(X))^\top]$ consists of finite entries, and is positive definite.
	\end{itemize}

Assumption (R0) tells us that to produce a sample $X\sim f_\theta$, we can first sample the seed $\omega \sim P$ and then transform the seed into $Y \defeq X_\theta(\omega)$. This procedure results in a sample equal in distribution: $X\overset d = Y$. Since $X_\theta$ is assumed to be continuously differentiable in $\theta$,  we have that the mapping $X_\theta(\cdot)$ is smooth as $\theta$ varies.

An idealized version of our algorithm works as follows: given $\thx$ computed from the original dataset, we first sample the seeds $\omega_1,\ldots, \omega_n\iid P$, and then while holding $\omega_1,\ldots, \omega_n$ fixed, solve for the value $\theta^*$ which satisfies:

\begin{equation}\label{eq:thetaStar}
\hat\theta\{X_{\theta^*}(\ul \omega)\} = \thx,
\end{equation}
thereby ensuring that the new sample has the same value of $\hat\theta$ as the original sample $\ul X$. Finally, we can produce our synthetic data $Y_i = X_{\theta^*}(\omega_i)$, which can be published for disclosure control or used for statistical inference. This idealized version has been previously employed when the statistic $\thx$ is a \emph{sufficient statistic}, and used mainly for the task of hypothesis testing \citep{lindqvist2005monte,lindqvist2007conditional,lindqvist2013exact,lillegard1999exact,engen1997stochastic,lillegard1999exact,taraldsen2018conditional}. However, when there is no low dimensional sufficient statistic, or when there does not exist an exact solution to Equation \eqref{eq:thetaStar}, such as in Example \ref{ex:noSolution}, the previous work is inapplicable.

\begin{ex}\label{ex:noSolution}
[No solution to Equation \eqref{eq:thetaStar}]
Suppose that $X_i = W_i + N_i$, where $W_i \iid \mathrm{Bern}(\theta)$ and $N_i \iid U(0,1)$ and $W_i \indep N_i$. Note that $\thx = n^{-1} \sum_{i=1}^n X_i - \frac 12 $ is an efficient estimator for $\theta$. Using the notation of (R0), $\Omega = U(0,1)^2$ and $X_\theta :
\Omega \rightarrow \RR$, defined by $X_\theta(u_1,u_2) = F^{-1}_{\theta}(u_1) + u_2$, where $F^{-1}_{\theta}(\cdot)$ is the quantile function for $\mathrm{Bern}(\theta)$. Then, for any value $\theta^* \in [0,1]$, $\thy = \hat\theta\left\{X_{\theta^*}(u^{(1)}_1,u^{(1)}_2),\ldots, X_{\theta^*}(u^{(n)}_1,u^{(n)}_{2})\right\}$ takes at most $n+1$ possible values:
$\left\{n^{-1} \sum_{i=1}^n u_2^{(i)} -\frac 12 + \frac{i}{n}\middle| i \in \{0,1,2,\ldots, n\}\right\}$. As such, given $\thx$ from an independent sample, the probability that $\thx$ is in this set is zero. From another perspective, the decimal value of $n\thy$ is completely determined by $n^{-1} \sum_{i=1}^n u_2^{(i)}$, and is not influenced by $\theta^*$.
\end{ex}

Based on Example \ref{ex:noSolution}, we cannot always expect to find an exact solution to Equation \eqref{eq:thetaStar}. In Theorem \ref{thm:onestep} and Algorithm \ref{alg:onestep}, we propose a one-step approximate solution which is based on the following intuition: Call $Z_i=X_{\thx}(\omega_i)$ and $Y_i = X_{\theta^*}(\omega_i)$, where $\theta^*$ is the exact solution to Equation \eqref{eq:thetaStar}. Since $\ul Z$ and $\ul Y$ use the same seeds, we expect that the differences $\thy-\theta^* = \thx-\theta^*$ and $\thz-\thx$ should be similar, especially for larger sample sizes. Solving $\thx-\theta^*=\thz-\thx$ for $\theta^*$ gives the solution $\theta^* = 2\thx-\thz$, which is exactly the one-step estimator of Theorem \ref{thm:onestep}.

The following Theorem shows that regardless of whether a solution to Equation \eqref{eq:thetaStar} exists, the one-step procedure preserves the efficient statistic up to $o_p(n^{-1/2})$.

\begin{restatable}{thm}{thmonestep}\label{thm:onestep}
Assume that (R0)-(R3) hold. Let $X_1,\ldots, X_n \iid f_{\theta_0}$ and let $\omega_1,\ldots, \omega_n\iid P$.

Set
$\theta^* =\mathrm{Proj}_{\Theta}\left(2 \thx-\thz\right)$, where  $\hat \theta$ is an efficient estimator and $(Z_i)_{i=1}^n = (X_{\thx}(\omega_i))_{i=1}^n$.

 Then for $(Y_i)_{i=1}^n = (X_{\theta_{new}}(\omega_i))_{i=1}^n$, we have  $\thy =  \thx + o_p(n^{-1/2})$.
\end{restatable}
\begin{comment}
\begin{proof}[Proof Sketch.]
The proof is based on two  expansions:
\[\thz - \thx = I^{-1}(\theta_0) \frac 1n \sum_{i=1}^n S(\theta_0,X_{\theta_0}(\omega_i)) + o_p(n^{-1/2})\] and
\[\thy=\theta^*+I^{-1}(\theta_0) \frac 1n\sum_{i=1}^n S(\theta_0, X_{\theta_0}(\omega_i)) + o_p(n^{-1/2}).\] As  $n\rightarrow \infty$, $\theta^*= 2\thx -\thz$ with probability one. Combining the results gives $\thy = \thx+o_p(n^{-1/2})$.
\end{proof}
\end{comment}

Theorem \ref{thm:onestep} shows that our one-step approach to synthetic data outperforms the parametric bootstrap in terms of the first result in Theorem \ref{thm:TVnaive}: whereas sampling from the fitted model results in estimators with inflated variance, the one-step approach gives a sample whose estimator $\thy$ is equal to $\thx$ up to an asymptotically negligible error of $o_p(n^{-1/2})$.

   \begin{comment}
    \begin{algorithm}\label{}
    \caption{One Step Synthetic Data Pseudo-Code}
    \scriptsize
    INPUT: Seed $\omega$, parametric family $\{f_\theta\mid \theta \in \Theta\}$, efficient estimator $\hat \theta(\cdot)$, and sample $X_1,\ldots, X_n \iid f_{\theta_0}$
    \begin{algorithmic}[1]
     \setlength\itemsep{0em}
     \STATE \texttt{set.seed}($\omega$) and a sample $Z_1,\ldots, Z_n \iid f_{\thx}$
      \STATE Choose $\theta^*\in \argmin_{\theta\in \Theta} \lVert \theta-( 2\thx-\thz)\rVert_2$
      \STATE \texttt{set.seed}($\omega$) and sample $Y_1,\ldots, Y_n \iid f_{\theta^*}$
    \end{algorithmic}
    OUTPUT: $Y_1,\ldots, Y_n$
    \label{alg:numerical}
    \end{algorithm}
    \end{comment}

\begin{rem}[Seeds]
 In the case of continuous real-valued random variables, we can be more explicit about the ``seeds.''
 Recall that for $U\sim U(0,1)$, $F^{-1}_{\theta}(U) \sim f_{\theta}$ where $F^{-1}_{\theta}(\cdot)$ is the quantile function. So in this case, the  distribution $P$ can be taken as $U(0,1)$, and $X_{\theta}(\cdot)$ can be replaced with $F^{-1}_{\theta}(\cdot)$.

When implementing the procedure of Theorem \ref{thm:onestep}, it may be convenient to use numerical seeds. For example in R, the command \texttt{set.seed} can be used to emulate the result of drawing $Z_i$ and $Y_i$ with the same seed $\omega_i$. In Algorithm \ref{alg:onestep}, we describe the one-step procedure in pseudo-code.
\end{rem}

    \begin{algorithm}[t]
    \caption{One-Step Synthetic Data Pseudo-Code in R}
    \label{alg:onestep}
        	\setstretch{1}
    	\small
\SetKwInput{KwInput}{Input}
\DontPrintSemicolon

  \KwInput{ Seed $\omega$, efficient estimator $\thx$, function \texttt{theta\_hat}(\ul y) to compute $\hat\theta(\ul y)$, function \texttt{rsample}($\theta$) which is continuously differentiable in $\theta$ which samples $n$ i.i.d. samples from $f_\theta$.\;}
  \texttt{set.seed}($\omega$)\;
 $\ul Z$ = \texttt{rsample}($\thx$)\;
 $\thz$ = \texttt{theta\_hat}($\ul Z$)\;

  $\theta^* = 2\thx-\thz$\;
  \If{$\theta^* \not\in \Theta$}
  {$\theta^*= \mathrm{Proj}_{\Theta}(\theta^*)$\;}
  \texttt{set.seed}($\omega$)\;
  $\ul Y$ = \texttt{rsample}($\theta^*$)\;
  \SetKwInput{KwOutput}{Output}
    \KwOutput{$Y_1,\ldots, Y_n$}
\end{algorithm}

	\begin{rem}[DP fully synthetic data]
	The one-step procedure results in a sample $\ul Y$, which is conditionally independent of $\ul X$, given $\thx$. This is called \emph{partially synthetic data} because all aspects of $\ul X$ are protected except for $\thx$. Partially synthetic data can be appropriate in some settings, but in others we may require $\ul Y$ to satisfy a stronger privacy guarantee such as differential privacy.
	We can easily use the one-step procedure to obtain a DP fully synthetic sample by using a DP efficient estimator as $\thx$. \citet{smith2011privacy} showed that under conditions similar to (R1)-(R3), there always exists such DP efficient estimators. By the post-processing property of DP the sample $\ul Y$ then has the same DP guarantee that $\thx$ has \citep[Proposition 2.1]{dwork2014algorithmic}. Theorem \ref{thm:onestep} is still valid when using a DP estimator for $\thx$ as only the efficiency of $\thx$ is used in the proof. In fact, the estimator $\hat\theta$ applied to the intermediate sample $\ul Z$ need not be the same one as used for $\thx$. For improved finite sample performance, it may be beneficial to use a non-private estimator for $\thz$ even when $\thx$ satisfies DP. See Section \ref{s:beta} where we investigate the performance generating of DP synthetic data using Algorithm \ref{alg:onestep}.
	\end{rem}

\section{Investigating the distribution of the one-step samples}\label{s:distribution}
In Section \ref{s:asymptotics}, we showed that the one-step approach to synthetic data solves one issue of the parametric bootstrap: whereas sampling from the fitted model results in statistics with inferior asymptotics, the one-step approach gives a sample with efficient estimators. The second problem of the parametric bootstrap was that the synthetic samples did not converge in total variation distance to the true joint distribution.

In this section, we provide results which give evidence that the distribution of the one-step samples approximates the true joint distribution.

Consider the one-step procedure: Let $\thx$ be given and draw $\omega_1,\ldots,\omega_n\iid P$. Set $Z_i = X_{\thx}(\omega_i)$ and call $\theta^* = 2\thx-\thz$. Finally, set $Y_i = X_{\theta^*}(\omega_i)$, which we know satisfies $\hat\theta \defeq \thy =\thx+o_p(n^{-1/2})$. Now, suppose that we knew the values $\theta^*$ and $\hat\theta$ beforehand and conditioned on them. The following Lemma shows that $Y_1,\ldots, Y_n\large|\{\theta^*,\hat\theta(\ul Y)=\hat\theta\} \approx f_{\theta^*}^n\{y_1,\ldots, y_n\mid \hat\theta(\ul y)=\hat\theta\}$. The approximation comes from modifying the procedure slightly. Given both $\theta^*,\hat\theta\in \Theta$, we produce $\ul Y^{\theta^*}$ as follows: sample $\omega_i\iid P$, and set $Y_i^{\theta^*} = X_{\theta^*}(\omega_i)$. For this modified procedure, we have exactly $Y_1^{\theta^*},\ldots, Y_n^{\theta^*}\large|(\hat\theta(\ul Y)=\hat\theta) \sim f_{\theta^*}^n(y_1,\ldots, y_n\mid \hat\theta(\ul y)=\hat\theta)$, which we prove in Lemma \ref{lem:conditional}.

\begin{restatable}{lemma}{lemconditional}\label{lem:conditional}
We assume (R0) and use the notation therein. Let $X_1,\ldots, X_n \iid f_{\theta_0}$ and let $\hat\theta,\theta^*\in \Theta$ such that there exists at least one solution in $\ul \omega$ to the equation  $\hat\theta\{X_{\theta^*}(\ul\omega)\}=\hat\theta$. Let $\omega_1,\ldots, \omega_n \iid P$ and call $Y_{i}^{\theta^*} = X_{\theta^*}(\omega_i)$. Then $Y_1^{\theta^*},\ldots, Y_n^{\theta^*} \Big|\left\{\hat\theta(\ul Y^{\theta^*})=\hat\theta\right\} \sim f_{\theta^*}^n\{y_1,\ldots, y_n \mid \hat\theta(\ul y) = \hat\theta\}$.
\end{restatable}

Lemma \ref{lem:conditional} suggests that the conditional distribution of $\ul Y$ is related to the conditional distribution of $\ul X$. Theorem \ref{thm:klGeneral} shows that when this is the case, the marginal distributions are also closely related. Theorem \ref{thm:klGeneral} demonstrates a phenomenon that efficient estimators are \emph{asymptotically sufficient}, a notion introduced by \citet{le1956asymptotic}. First, we need an additional assumption on the distribution of $\hat\theta(\ul X)$.

\begin{itemize}
\item [(R4)] Let $\thx$ be an efficient estimator of $\theta$.	Let $g_{\theta,n}(\hat\theta)$ be the marginal distribution of $\hat\theta(\ul X)$ based on the sample $X_1,\ldots, X_n \iid f_\theta$. Assume there exists functions $G_{ijk}(\hat\theta)$ such that  $\left| \frac{\partial^3}{\partial \theta_i\partial \theta_j\partial \theta_k} \log g_{\theta,n}(\hat\theta)\right|\leq nG_{ijk}(\hat\theta)$ for all $\hat\theta$, all $n\geq 1$,  and all $\theta \in B(\theta_0)$, where $\EE_{\hat\theta \sim \theta}G_{ijk}(\hat\theta)<\infty$.
\end{itemize}

\begin{restatable}{thm}{thmklGeneral}\label{thm:klGeneral}
Under assumptions (R0)-(R4), let $\theta\in \Theta$ and let $\theta_n$ be a sequence of values in $\Theta$. Let $\hat\theta(\cdot)$ be a randomized estimator based on a sample $X_1,\ldots, X_n \iid f_\theta$, with conditional distribution $g_n(\hat\theta\mid \ul X)$ and marginal distribution $g_{\theta,n}(\hat\theta)$. Then the KL divergence between the marginal distributions of $X_1,\ldots, X_n$ and $Y_1,\ldots,Y_n \sim f^n_{\theta_n}\{y_1,\ldots, y_n \mid \hat\theta(\ul y)=\hat\theta(\ul X)\}$, where $\hat\theta(\ul X)\sim g_{\theta,n}(\hat\theta)$ is
\begin{align*}
\mathrm{KL}\left(X_1,\ldots,X_n\middle|\middle| Y_1,\ldots, Y_n\right)

=o(n) \lVert \theta_n-\theta\rVert^2 +O(n) \lVert \theta_n-\theta\rVert^3.\end{align*}
In particular, if $\theta_n-\theta=O(n^{-1/2})$,

then the above KL divergence goes to zero as $n\rightarrow \infty$.
\end{restatable}

Noting that $2\thx-\thz =\theta+O_p(n^{-1/2})$, combining Lemma \ref{lem:conditional} with Theorem \ref{thm:klGeneral} along with the discussion at the beginning of this section suggests that the distribution of the one-step synthetic data approaches the distribution of the original $\ul X$ as the sample size increases giving consistent synthetic data. This suggests that the one-step synthetic data avoids the problem of property 2 of Theorem \ref{thm:TVnaive}, which showed that the parametric bootstrap resulted in inconsistent synthetic data.

In Section \ref{s:burr} we show that in the case of the Burr distribution, the one-step synthetic data is indistinguishable from the true distribution in terms of the K-S test, whereas for the parametric bootstrap the K-S test has significantly higher power, indicating that the one-step synthetic data is asymptotically consistent, whereas the parametric bootstrap synthetic data is inconsistent.

\begin{rem}[Application to Hypothesis Testing]
Because these results suggest that the one-step samples are approximately drawn from the conditional distribution of $\ul X\mid \thx$, the one-step synthetic data can potentially be used to produce approximate $p$-values for tests of goodness-of-fit. Suppose we want to test $H_0: \theta\in \Theta_0$. Recall that if there were a sufficient statistic for $\theta$ under $H_0$, conditioning on this statistic can be used to derive a valid hypothesis test. In fact, in many classical problems the optimal hypothesis tests are of this form \citep[Theorem 4.123]{schervish2012theory}. \citet{barber2020testing} showed that when sufficient statistics do not exist, an efficient estimator can be used to result in asymptotically valid hypothesis tests. Due to this, our one-step samples may be well suited to this problem.

We apply the one-step synthetic data to a hypothesis testing problem in Section \ref{s:DPtesting} showing using the one-step samples gives a test with more accurate type I errors and higher power than compared to using parametric bootstrap samples. This example provides additional numerical support for the consistency of the one-step samples.
\end{rem}

\begin{rem}[Testing based on Differential Private Statistics]
    One application where the one-step synthetic data may be particularly advantageous is in testing problems based on a differential private output. Suppose that the data analyst releases privatized summary statistics, and we want to use these statistics to test certain hypotheses about the generating distribution of the private data. However, despite the fact that many DP estimators have asymptotically negligible noise, standard asymptotic tests often have unacceptably poor performance with DP data \citep{wang2018statistical}. Furthermore, deriving the exact marginal likelihood of the DP output often requires evaluating a high dimensional and intractable integral \citep{williams2010probabilistic}, integrating over the space of possible databases, which makes the methods of \citet{barber2020testing} inapplicable. While parametric bootstrap tests may be applied, they do not guarantee accurate type I error rates, as discussed in Section \ref{s:bootstrap}. The one-step algorithm is nearly as easy to implement as the  parametric bootstrap, yet as we show in Section \ref{s:DPtesting} gives accurate $p$-values, even in moderate sample sizes.

\end{rem}

		\section{Examples and simulations}\label{s:simulations}
		In this section, we demonstrate the performance of the one-step algorithm on several simulations. In Sections \ref{s:burr} and \ref{s:log-linear}, we produce synthetic data from the Burr distribution as well as a log-linear model. In Section \ref{s:beta}, we produce differentially private synthetic data for the beta distribution. In Section \ref{s:DPtesting} we use our methods to derive an asymptotically valid test for the difference of proportions under differential privacy. In Section \ref{s:regression}, we compare the running time of the one-step approach compared to the MCMC approach of \citet{barber2020testing} when implemented on a linear regression problem.
		\subsection{Burr type XII distribution}\label{s:burr}

		The Burr Type XII distribution, denoted $\mathrm{Burr}(c,k)$, also known as the Singh–Maddala distribution, is a useful model for income \citep{mcdonald2008some}. The distribution has pdf $f(x)=ckx^{c-1}(1+x^c)^{-(k+1)}$, with support $x>0$. Both $c$ and $k$ are positive. The Burr distribution was chosen for our first simulation because 1) it is one-dimensional, allowing for the Kolmogorov-Smirnov (K-S) test to be applied, and 2) as it is not exponential family or location-scale,  the conditional distribution is non-trivial. The goal of the following simulation is to verify that the samples generated using Algorithm \ref{alg:onestep} are indistinguishable from the original unknown distribution, as tested via the K-S test.

			For the simulation, we set $c=2$ and $k=4$, and denote $\theta = (c,k)$. Let $\hat \theta_{MLE}$ be the maximum likelihood estimator (MLE).

			We draw $X_i \iid \mathrm{Burr}(2,4)$, $Z_i \iid \mathrm{Burr}\{\hat \theta_{MLE}(\ul X)\}$, and  $(Y_i)_{i=1}^n$ from Algorithm \ref{alg:onestep}.
			The simulation is conducted for $n\in \{100,1000,10000\}$ with results averaged over 10000 replicates for each $n$.

			We calculate the empirical power of the K-S test, comparing each sample with the true distribution $\mathrm{Burr}(2,4)$, at type I error $.05$. The results are presented in Table \ref{tab:burrKS}. We see that the $(X_i)$ have empirical power approximately $.05$, confirming that the type I error is appropriately calibrated. We also see that the K-S test using  $(Y_i)$ has power approximately $.05$, indicating that the empirical distribution of the one-step samples $(Y_i)$ is very close to the true distribution. On the other hand, we see that the K-S test with $(Z_i)$ has power $.15$, significantly higher than the type I error, indicating that the parametric bootstrap samples $(Z_i)$ are from a fundamentally different distribution than the $(X_i)$. This result is in agreement with Theorem \ref{thm:TVnaive} and the results of Section \ref{s:distribution}.

\begin{figure}
\begin{minipage}[c]{\linewidth}
\captionof{table}{Empirical power of the Kolmogorov-Smirnov test for $\mathrm{Burr}(2,4)$ at type I error $.05$. $(X_i)$ are drawn i.i.d from $\mathrm{Burr}(2,4)$, $(Z_i)$ are drawn i.i.d from $\mathrm{Burr}\{\thx\}$, and $(Y_i)$ are from Algorithm \ref{alg:onestep}. Results are averaged over 10000 replicates, for each $n$. Standard errors are approximately $0.0022$ for lines 1 and 3, and  $0.0036$ for line 2.}
		\label{tab:burrKS}
		\centering
		\begin{tabular}{llll}

	\phantom{x}	$n$:	& 100 & 1000 & 10000 \\
			\hline
			$(X_i)$  & 0.0471 & 0.0464 & 0.0503 \\
			$(Z_i)$  & 0.1524 & 0.1541 & 0.1493 \\
			$(Y_i)$  & 0.0544 & 0.0489 & 0.0485 \\

			\hline
		\end{tabular}
		\end{minipage}
	\end{figure}

\subsection{Log-linear model}\label{s:log-linear}

		This example is based on a dataset of of 68,694 passengers in automobiles and light trucks involved in accidents in the state of Maine in 1991. Table \ref{tab:seatbelt} tabulates passengers according to gender (G), location (L), seatbelt status (S), and injury status (I).

		As in \citet{agresti2003categorical}, we fit a hierarchical log-linear model based on all one-way effects and two-way interactions. The model is summarized in Equation \eqref{eq:loglinear}, where $\mu_{ijk\ell}$ represents the expected count in bin $i,j,k,\ell$.  The  parameter $\la_i^G$ represents the effect of Gender, and parameter $\la_{ij}^{GL}$ represents the interaction between Gender and Location. The other main effects and interactions are analogous.

		\begin{equation}\label{eq:loglinear}
		\log \mu_{ijk\ell} = \la + \la_i^G + \la_j^L +\la_k^S + \la_\ell^I
		+\la_{ij}^{GL}+\la_{ik}^{GS}+\la_{i\ell}^{GI}
		+\la_{jk}^{LS}+\la_{j\ell}^{LI}+\la_{k\ell}^{SI}\end{equation}

		For our simulations, we	treat the fitted parameters as the true parameters, to ensure that model assumptions are met. We simulate from the fitted model at sample sizes $n\in \{10^2,10^3,10^4,10^5\}$ and compare the performance in terms of the fitted probabilities for each bin of the contingency table. 		The results are plotted in Figure \ref{fig:log-linear}, with both axes on log-scale. The ``mean error'' is the average squared $\ell_2$ distance between the estimated parameter vector and the true parameter vector, averaged over 200 replicates. To interpret the plot, note that if the error is of the form $\mathrm{error} = cn^{-1}$, where $c$ is a constant, then $\log(\mathrm{error}) = c+(-1) \log(n)$. So, the slope represents the convergence rate, and the vertical offset represents the asymptotic variance. In Figure \ref{fig:log-linear}, we see that the curve for $\thy$, based on one-step samples,  approaches the curve for $\thx$, indicating that they have the same asymptotic rate and variance. On the other hand, the curve for  $\thz$, based on parametric bootstrap samples, has the same slope, but does not approach the $\thx$ curve, indicating that $\thz$ has the same rate but inflated variance.

		Recall that our procedure approximately preserves the sufficient statistics, similar to sampling from a conditional distribution. 	Previous work has proposed procedures to sample directly from conditional distributions for contingency table data. However, these approaches require sophisticated tools from algebraic statistics, and are computationally expensive (e.g., MCMC) \citep{karwa2013conditional}. In contrast, our approach is incredibly simple to implement and highly computationally efficient. Our approach is also applicable for a wide variety of models, whereas the techniques to sample directly from the conditional distribution often require a tailored approach for each setting.

	\begin{figure}
	\begin{minipage}[c]{\linewidth}
				\captionof{table}{Injury, Seat-Belt Use, Gender, and Location. Source: \citet[Table 8.8]{agresti2003categorical}. Originally credited to Cristanna Cook, Medical Care Development, Augusta, Maine.}
		\label{tab:seatbelt}
		\centering
		\begin{tabular}{lllrr}
			\multicolumn{3}{c}{\textbf{}} &   \multicolumn{2}{c }{\textbf{Injury}} \\
			\cmidrule(ll){4-5}
			\textbf{Gender}&\textbf{Location}&\textbf{Seatbelt}&\multicolumn{1}{c}{No}&\multicolumn{1}{c}{Yes}\\
			\midrule

			Female&Urban&No&7,287&996\\
			&&     Yes&11,587&759\\
			&Rural&No&3,246&973\\
			&&     Yes&6,134&757\\
			Male&Urban&No&10,381&812\\
			&&Yes&10,969&380\\
			&Rural&No&6,123&1,084\\
			&&Yes&6,693&513\\
			\bottomrule
		\end{tabular}
\end{minipage}
	\end{figure}

\subsection{Differentially private beta distributed synthetic data}\label{s:beta}

				\begin{figure}[ht]
		  \begin{subfigure}[t]{0.48\textwidth}

		\includegraphics[width=\textwidth]{Contingency.pdf}
		\caption{Simulations  corresponding to the log-linear model with two-way interactions from Section  \ref{s:log-linear}}
		\label{fig:log-linear}
		\end{subfigure}
		\hspace{.02\textwidth}
		\begin{subfigure}[t]{0.48\textwidth}
			\includegraphics[width=\textwidth]{DP_Beta2.pdf}
			\caption{Simulations for the beta distribution from Section \ref{s:beta}. $\thx$ is the MLE. $\thz$ and $\thy$ both satisfy 1-DP.}
			\label{fig:beta}
			\end{subfigure}
			\caption{Average squared $\ell_2$-distance between the estimated parameters and the true parameters on the log-scale. Averages are over 200 replicates for both plots. $\thx$ is from the true model, $\thz$ from the fitted model, and $\thy$ from Algorithm \ref{alg:onestep}.}
		\end{figure}

		In this example, we assume that $X_1,\ldots, X_n \iid \mathrm{Beta}(\alpha,\beta)$, where $\alpha,\beta\geq 1$, and our goal is to produce differentially private (DP) synthetic data. 	Often, to ensure finite sensitivity, the data are clamped to artificial bounds $[a,b]$, introducing bias in the DP estimate. Typically, these bounds are fixed in $n$, resulting in asymptotically negligible noise, but $O_p(1)$ bias. In Appendix \ref{s:betaProof} we show that it is possible to increase the bounds in $n$ to produce both noise and bias of order $o_p(n^{-1/2})$, resulting in an efficient DP estimator. In this section, we show through simulations that using this estimator along with Algorithm \ref{alg:onestep} results in a DP sample with optimal asymptotics.

		\begin{comment}
		Recall that  $n^{-1} \sum_{i=1}^n \log(X_i)$ and $n^{-1}\sum_{i=1}^n \log(1-X_i)$ are sufficient statistics for the beta distribution. We will add Laplace noise to each of these statistics to achieve differential privacy. However, the sensitivity of these quantities is unbounded. First we pre-process the data by setting
		$\twid X_i = \min\{\max(X_i,t),1-t\},$
		where $t$ is a threshold that depends on $n$. Then the $\ell_1$-sensitivity of the pair of sufficient statistics is
		$\Delta(t) =2n^{-1}\left|\log(t)-\log(1-t)\right|$.
		We add independent noise to each of the statistics from the distribution $\mathrm{Laplace}\{\Delta(t)/\epsilon\}$, which results in $\ep$-DP versions of these statistics. Finally, we estimate $\theta=(\alpha,\beta)$ by plugging in the privatized sufficient statistics into the log-likelihood function and maximizing over $\theta$. The resulting parameter estimate satisfies $\ep$-DP by post-processing.

		We must carefully choose the threshold $t$ to ensure that the resulting estimate is efficient. The choice of $t$ must satisfy $\Delta(t) = o(n^{-1/2})$ to ensure that the noise does not affect the asymptotics of the  likelihood function. We also require that both $P(X_i<t)=o(n^{-1/2})$, and $P(X_i>1-t)=o(n^{-1/2})$ to ensure that $\twid X_i =X_i + o_p(n^{-1/2})$, which limits the bias to  $o_p(n^{-1/2})$. For the beta distribution, we can calculate that $P(X_i<t) = O(t^\alpha)$ and $P(X_i>1-t) = O(t^\beta)$. Since we assume that $\alpha,\beta\geq 1$, so long as $t=o(n^{-1/2})$ the probability bounds will hold. Taking $t = \min[1/2,10/\{\log(n) \sqrt n\}]$ satisfies $t=o(n^{-1/2})$, and we estimate the sensitivity as
		\[\Delta(t)\leq 2n^{-1} \log(t^{-1})\leq 2n^{-1}\log\{\log(n)\sqrt n\} = O\{\log(n)/n\} = o(n^{-1/2}),\]
		which satisfies our requirement for $\Delta$. While there are many choice of $t$ which satisfy the requirements, our threshold was chosen to optimize the finite sample performance, so that the asymptotics could be illustrated with smaller sample sizes.
		\end{comment}

		For the simulation, we sample $X_1,\ldots, X_n \iid \mathrm{Beta}(5,3)$ for $n \in \{10^3,10^4,10^5,10^6\}$. We estimate $\thx$ with the MLE. Using $\epsilon=1$, we privatize the sufficient statistics as described in Appendix \ref{s:betaProof}, and obtain $\hat \theta_{DP}$ from the privatized log-likelihood function. We sample $Z_1,\ldots, Z_n \iid f_{\hat \theta_{DP}}$ and calculate the MLE $\thz$. We draw $(Y_i)_{i=1}^n$ from Algorithm \ref{alg:onestep} using $\hat\theta_{DP}$ in place of $\thx$.  In Figure \ref{fig:beta}, we plot the mean squared $\ell_2$ error between each estimate of $\theta$ from the true value $(5,3)$, over 200 replicates. From the log-scale plot, we see that $\hat\theta_{DP}$ and $\thy$ have the same asymptotic performance as the MLE, whereas $\thz$ has inflated variance. See Example \ref{s:log-linear} for more explanation of this interpretation.

		\subsection{DP two sample proportion test}\label{s:DPtesting}

    In this section, we illustrate how the one-step samples can be used to perform approximate hypothesis tests. We base the simulation on a problem in differential privacy, where we are given only access to DP summary statistics and are tasked with testing a hypothesis on the generating distribution for the (missing) private data. Such settings result in complex distributions which can be difficult to work with directly, making MCMC methods such as in \citet{barber2020testing} cumbersome and potentially intractable. We show through a simulation that the one-step samples give highly accurate $p$-values and improved power compared to the parametric bootstrap.

		Suppose we have two independent samples of binary data, one from a ``control population'' and another from a ``treatment population.'' We denote $X_1,\ldots, X_n\iid \mathrm{Bern}(\theta_X)$ as the control sample, and $Y_1,\ldots, Y_m \iid \mathrm{Bern}(\theta_Y)$ as the treatment sample.

		It was shown in \citet{awan2018differentially} and \citet{awan2020differentially} that the \emph{Tulap} distribution is the optimal DP mechanism to generate uniformly most powerful hypothesis tests and uniformly most accurate confidence intervals for Bernoulli data. Thus, to satisfy $\ep$-differential privacy, the data curators release the following noisy statistics: $\twid X = \sum_{i=1}^n X_i + N_1$ and $\twid Y = \sum_{i=1}^m Y_i + N_2$, where $N_1,N_2 \iid \mathrm{Tulap}\{0,\exp(-\ep),0\}$ and the sample sizes $m$ and $n$ are released without modification.

		Recall that $N_i \overset d=  G_1-G_2+U$, where $G_1,G_2 \iid \mathrm{Geom}\{1-\exp(-\ep)\}$ and $U\sim \mathrm{Unif}(-1/2,1/2)$. We can think of $\twid X$ and $\twid Y$ as noisy counts for the control and treatment groups, respectively.

		Based only on the privatized summary statistics $\twid X$ and $\twid Y$, we are to test  $H_0: \theta_X= \theta_Y$ versus $H_1: \theta_X\leq \theta_Y$. Without privacy, there exists a uniformly most powerful test, which is constructed by conditioning on the total number of ones: $\sum_{i=1}^n X_i + \sum_{j=1}^m Y_j$, a complete sufficient statistic under the null hypothesis. However, with the noisy counts, it can be verified that there is not a low-dimensional sufficient statistic. On the other hand, an efficient estimator for $\theta_X=\theta_Y$ under the null hypothesis is $\hat\theta(\twid X,\twid Y) = \min[\max\{(\twid X + \twid Y)/(m+n),0\},1]$. Note that deriving the exact distribution of $(\twid X,\twid Y)\mid \hat\theta(\twid X,\twid Y)$ is fairly complex, involving the convolution of distributions. However, the one-step method can easily produce approximate samples from this conditional distribution. In what follows, we use the one-step algorithm, given in Algorithm \ref{alg:onestep} and investigate the properties of a hypothesis test based on this conditional distribution in comparison with a parametric bootstrap test.

		Recall that without privacy, the uniformly most powerful test uses the test statistic $Y$, and threshold computed from the conditional distribution of $Y|X+Y$ under the null hypothesis. With privacy, we use the test statistic $\twid Y$, and compute the $p$-values based on $\twid Y \mid \hat \theta(\twid X,\twid Y)$. In particular, we compare the performance of this test versus the parametric bootstrap test, which uses the test statistic $\twid Y$ based on the approximate sampling distribution, which is the convolution of $\mathrm{Binom}\{m,\hat\theta(\twid X,\twid Y)\}$ and $\mathrm{Tulap}\{0,\exp(-\ep),0\}$.

		For the simulation, we use the sample size $n=m=200$, fix $\theta_X=0.3$, set the privacy parameter to $\ep=1$, and base the simulation on 10,000 replicates.  Under the null hypothesis, where $\theta_Y=\theta_X$, we plot the empirical cumulative distribution function (CDF) for the $p$-values of the  proposed test as well as for the parametric bootstrap test in Figure \ref{DP2prop}(a). Recall that a properly calibrated $p$-value will have the CDF of $U(0,1)$. We see that the empirical CDF of the $p$-values for the one-step test closely approximate the ideal CDF, whereas the $p$-values of the parametric bootstrap test are overly conservative.

		Next, we study the power of the one-step test versus the parametric bootstrap test in Figure \ref{DP2prop}(b). For this simulation, we set $n=m=200$, fix $\theta_X=.3$, and set $\ep=1$. We vary the value of $\theta_Y$ along the $x$-axis by increments of $.001$ and plot the empirical power of the two tests, averaged over 10,000 replicates for each value of $\theta_Y$. We see that the one-step test offers a considerable increase in power over the bootstrap test.

		\begin{figure}[t]
		\begin{subfigure}[t]{0.48\textwidth}
			\includegraphics[width=\textwidth]{Pvalue_200_200_0.3_pvalues.pdf}
		\caption{Under the null  hypothesis, the empirical cumulative distribution of both tests. $\theta_X=\theta_Y=.3$. Results based on 10,000 replicates.  }
			\label{fig:pvalue}
			\end{subfigure}
				\hspace{.02\textwidth}
		  \begin{subfigure}[t]{0.48\textwidth}

		\includegraphics[width=\textwidth]{Power_200_200_0.3_powers.pdf}
		\caption{Empirical power at type I error $.05$. $\theta_X=.3$, $\theta_Y$ is varied along the $x$-axis. Results averaged from 10,000 replicates for each value of $\theta_Y$.}
		\label{fig:power}
		\end{subfigure}
		\caption{Simulations for the DP two sample proportion test of Section \ref{s:DPtesting}. In red is the one-step test, and in blue is the parametric bootstrap test. Sample sizes are $n=m=200$, privacy parameter is $\ep=1$, and type I error is $.05$.}
		\label{DP2prop}
		\end{figure}

		\subsection{Running time versus MCMC for linear regression}\label{s:regression}
		In this section, we compare the running time of Algorithm \ref{alg:onestep} to an MCMC method of sampling the conditional distribution, proposed by \citet{barber2020testing}. We base the comparison on a linear regression problem.

		We use a similar setting as Simulation 1 in \citet{barber2020testing}. We sample $Y_i \iid N(Z_i \beta,1)$, where $\beta$ and $Z_i$ are 5-dimensional, with the true (unknown) parameter $\beta^\top = (.2,.2,.2,.2,.2)$, and $Z^\top_i \iid N(0,I)$. For synthetic data, we condition on an efficient estimator $\hat\beta$ and the matrix $Z$. The approach of \citet{barber2020testing} requires adding a small amount of independent noise $w \sim N(0, \sigma^2 I)$ to the score function before minimizing. To ensure that the estimator is still efficient, we set $\sigma = 1/n$. According to \citet{barber2020testing}, the pdf of the conditional distribution is proportional to
		\begin{equation}\label{eq:condPDF}f_{\hat\beta}(Y|Z) \exp\left\{-\frac{\lVert \nabla\ell(\hat\beta\mid Y,Z)\rVert^2}{2\sigma^2d}\right\} \det\left\{\nabla^2
		\ell(\hat\beta \mid Y,Z)\right\},\end{equation}
		where $\ell(\beta\mid Y,Z)$ is the log-likelihood function. This unnormalized density can then  be used in an MCMC sampler. While this is a relatively simple model, sampling from the conditional distribution is still non-trivial as we are tasked with sampling from an $n$-dimensional space where the mass is concentrated around an $n-5$ dimensional subspace.

		For each $n\in\{ 2^6, 2^7, \dots, 2^{14}\}$, we simulate a sample $(Y,Z)$ of the appropriate dimension. To do this, we run 100 rounds of a one-at-a-time MCMC method with Gaussian proposals, which consists of 100*n one-dimensional update steps, using Equation \eqref{eq:condPDF} and step-size of .1. We then run 1000 replicates of the one-step method of Algorithm \ref{alg:onestep}. In Figure \ref{fig:runningTime}, we average the replicates to get the average time for one round of the MCMC, as well as the average run tixme of Algorithm \ref{alg:onestep}.

				\begin{figure}[t]
		    \begin{subfigure}[t]{0.48\textwidth}
			\includegraphics[width=\textwidth]{runningTime.pdf}
		\caption{The `$\times$'s represent the computed datapoints, whereas the lines are interpolations. }
			\label{fig:realTime}
			\end{subfigure}
				\hspace{.02\textwidth}
		  \begin{subfigure}[t]{0.48\textwidth}
		  	\includegraphics[width=\textwidth]{log_runningTime.pdf}
		\caption{Both axes are log-scale. In orange is a line with slope 2, and in green is a line of slope 1.}
		\label{fig:logTime}
		\end{subfigure}
		\caption{Linear regression simulation of Section \ref{s:regression}. Running time for one round of a one-at-a-time MCMC in red, and average time for one run of Algorithm \ref{alg:onestep} in blue.}
		\label{fig:runningTime}
		\end{figure}

		In Figure \ref{fig:runningTime}a, we see that even a single round of the MCMC is much more computationally expensive than a complete run of Algorithm \ref{alg:onestep}. Through additional simulations at smaller sample sizes ($n\approx 100$), we found that approximately $35,000$ rounds of the MCMC algorithm were required to produce a sample with less than $.01$ correlation with the initial state, when the step size was optimized to have between $10-20\%$ acceptance rate. Thus, to produce a single sample at $n=16384$, that is minimally correlated with the original $Y$, the MCMC method is expected to take 295610s $\approx$ 82hours. On the other hand, Algorithm \ref{alg:onestep} only takes 0.002538s to generate a sample at $n=16384$. If these simulation methods are to be used to produce $p$-values, then thousands of samples are needed. Furthermore, in real-world datasets $n$ and $p$ may be far larger than those considered in this simulation, making the MCMC implementation even more impractical.

		In Figure \ref{fig:runningTime}b, we have another plot of the running time of the MCMC procedure as well as Algorithm \ref{alg:onestep}, but with both axes on the log-scale. As discussed in Section \ref{s:log-linear}, this sort of plot allows us to infer the computational complexity. In green is a line of slope 1, which represents the computational complexity $O(n)$, and we see that the running time of Algorithm \ref{alg:onestep} has approximately the same slope as this line. In orange is a line of slope 2, which represents the computation complexity of $O(n^2)$, and we see that the running time of one round of the MCMC procedure has approximately the same slope as this line. Thus, not only is the observed running time of Algorithm \ref{alg:onestep} much smaller than for the MCMC, but the computation complexity is reduced from $O(n^2)$ to $O(n)$.

		\begin{rem}[Computational Complexity]
		    The computation complexity of Algorithm \ref{alg:onestep} is expected to be $O(n)$ for many common models, as it requires simulating samples of size $n$, and computing an efficient estimator for the parameter, which can often be done in $O(n)$ time. On the other hand, MCMC procedures are expected to take at least $O(n^2)$ time, as they require at least $O(n)$ steps to converge (consider for instance, the one-at-a-time MCMC) as well as evaluating the likelihood at each step, which often takes at least $O(n)$ time. This intuition is confirmed by the results in Figure \ref{fig:runningTime}.

		    This analysis assumes that the likelihood
		    can be evaluated in $O(n)$ time. With latent variables or missing-data, the likelihood may be far more expensive, or even
		    intractable.
		\end{rem}

	\section{Discussion} \label{s:conclusions}
	In this paper, we proposed a simple method of producing synthetic data from a parametric model, which approximately preserves efficient statistics ensuring optimal asymptotics.  We also provided evidence in Section \ref{s:distribution} that the one-step synthetic data results in a distribution which approaches the true underlying distribution. Both of these properties are in contrast to the common approach of sampling from a fitted model, which we showed results in inefficient estimators and inconsistent synthetic data. Our one-step approach is also widely applicable to parametric models, requiring mild regularity conditions, and is both easily implemented and highly computationally efficient. Our approach also allows for both partially synthetic data, as well as differentially private fully synthetic data by incorporating a DP efficient estimator.

     Besides synthetic data, there is also promise for using the one-step approach for hypothesis tests as well. \citet{barber2020testing} showed that one can produce powerful and accurate hypothesis tests by conditioning on an efficient statistic for the null model. However, their approach is likelihood based and requires MCMC methods for implementation. In problems with high-dimensional latent variables such an approach is inapplicable, for example in differential privacy where the entire private database is  latent.  The one-step approach offers a computationally efficient alternative, which we demonstrated in Section \ref{s:DPtesting} gives a more accurate and powerful test that the parametric bootstrap. We also showed in Section \ref{s:regression} that the one-step method has a much improved computational complexity compared to the method of \citet{barber2020testing}, even in the simple case of linear regression.

    In Section \ref{s:distribution}, we studied the distributional properties of a modified version of the one-step method, which suggested that the one-step synthetic data converges to the true joint distribution as the sample size increases. 	We also saw in Section  \ref{s:burr} that in the case of the Burr distribution, the Kolmogorov-Smirnov test cannot distinguish between our synthetic sample and the true distribution which generated the original sample, supporting the conjecture that the one-step sample is consistent. In future work, we propose to formally prove that the total variation distance between the one-step samples and the true distribution goes to zero as the sample sizes goes to infinity. This will provide additional theoretical justification for using the one-step method for synthetic data and hypothesis tests.

    While our approach was focused on parametric models, similar to Theorem \ref{thm:TVnaive}, there is a  loss of ``efficiency'' when sampling from a non-parametric model as well. In the spirit of our approach in this paper, it may be advantageous to sample from a non-parametric model conditional on the sample having ``similar'' estimates as the original data. An interesting future direction would be to formalize and investigate this direction.

    \section*{Acknowledgements}
    This work was supported in part by NSF grants SES-1853209 awarded to The Pennsylvania State University. The authors are very grateful to Roberto Molinari, Matthew Reimherr, and Aleksandra Slavkovi\'c for several helpful discussions, as well as feedback on early drafts.

\appendix

	\section{Background on differential privacy}\label{s:dp}

	In this section, we review the basics of differential privacy (DP), which was originally proposed in \citet{dwork2006calibrating} as a framework to develop methods of preserving privacy, with mathematical guarantees. Intuitively, the constraint of differential privacy requires that for all possible databases, the change in one person's data does not significantly change the distribution of outputs. Consequently, having observed the DP output, an adversary cannot accurately determine the input value of any single person in the database. Definition \ref{def:dp} gives a formal definition of DP. In Definition \ref{def:dp}, $h: \mscr X^n \times \mscr X^n \rightarrow \ZZ^{\geq 0}$ represents the \emph{Hamming metric}, defined by $h(\ul x,\ul x') =\# \{i \mid x_i\neq x'_i\}$.

	\begin{definition}
		[Differential privacy: \citet{dwork2006calibrating}]\label{def:dp}
		Let the \emph{privacy parameter} $\ep>0$ and the sample size $n\in \{1,2,\ldots\}$ be given. Let $\mscr X$ be any set, and $(\mscr Y, \mscr S)$  a measurable space. Let $\mscr M = \{M_{\ul x} \mid \ul x\in \mscr X^n\}$ be a set of probability measures on $(\mscr Y, \mscr S)$, which we call a \emph{mechanism}. We say that $\mscr M$ satisfies \emph{$\ep $-differential privacy} ($\ep$-DP) if $M_{\ul x}(S) \leq e^\ep M_{\ul x'}(S)$ for all $S\in \mscr S$ and all $\ul x,\ul x' \in \mscr X^n$ such that $h(\ul x,\ul x')=1$.
	\end{definition}

	An important property of differential privacy is that it is invariant to post-processing. Applying any data-independent procedure to the output of a DP mechanism preserves $\ep$-DP \citep[Proposition 2.1]{dwork2014algorithmic}. Furthermore, \citet{smith2011privacy} demonstrated that under conditions similar to (R1)-(R3), there exist efficient DP estimators for parametric models. Using these techniques, the one-step procedure can produce DP synthetic data by using a DP efficient statistic.

	\begin{comment}	\begin{cor}[Differentially private synthetic data]\label{cor:dp}
		Assume that (R0)-(R3)  hold. Let $X_1,\ldots, X_n \iid f_\theta$, and let $\omega_1,\ldots, \omega_n\iid P$.
		Set $\theta^* \in \argmin_{\theta\in \Theta} \left\lVert \theta-\left( 2\thx-\hat\theta(\ul Z)\right)\right\rVert_2$,
		where  $\thx$ is an $\ep$-DP efficient estimator, $\hat \theta$ is efficient, and  $(Z_i)_{i=1}^n =(X_{\thx}(\omega_i))_{i=1}^n$. Releasing $(Y_i)_{i=1}^n = (X_{\theta^*}(\omega_i))_{i=1}^n$  satisfies $\ep$-DP, and $\hat \theta(\ul Y)=\thx+o_p(n^{-1/2})$.
	\end{cor}\end{comment}

	\begin{rem}
		Besides Definition \ref{def:dp}, there are many other variations of differential privacy, the majority of which are relaxations of Definition \ref{def:dp} which also allow for efficient estimators. For instance, approximate DP \citep{dwork2006our}, concentrated DP \citep{Dwork2016:ConcentratedDP,Bun2016ConcentratedDP}, truncated-concentrated DP \citep{bun2018composable}, and Renyi DP \citep{Mironov2017} all allow for efficient estimators. On the other hand, local differential privacy \citep{kasiviswanathan2011,Duchi2013FOCS} in general does not permit efficient estimators and would not fit in our framework.  For an axiomatic treatment of formal privacy, see \citet{Kifer2012axiomatic}.
	\end{rem}

	One of the earliest and simplest privacy mechanisms is the \emph{Laplace mechanism}.

	Given a statistic $T$, the Laplace mechanism adds independent Laplace noise to each entry of the statistic, with scale parameter proportional to the \emph{sensitivity} of the statistic. Informally, the sensitivity of $T$ is the largest amount that $T$ changes, when one person's data is changed in the dataset.

	\begin{prop}[Sensitivity and Laplace: \citet{dwork2006calibrating}]\label{prop:laplace}
		Let the privacy parameter $\ep>0$ be given, and let $T: \mscr X^n \rightarrow \RR^p$ be a statistic. The \emph{$\ell_1$-sensitivity} of $T$ is
		$\Delta_n(T) = \sup \lVert T(\ul x) - T(\ul x')\rVert_1,$
		where the supremum is over all $\ul x,\ul x'\in \mscr X^n$ such that $h(\ul x,\ul x')=1$.
		Provided that $\Delta_n(T)$ is finite,
		releasing the vector $\left\{T_j(\ul x) + L_j\right\}_{j=1}^p$ satisfies $\ep$-DP, where
		$L_1,\ldots, L_p \iid \mathrm{Laplace}\left\{\Delta_n(T)/\ep\right\}$.
	\end{prop}

	\section{Deriving an DP efficient estimator for the beta distribution}\label{s:betaProof}

		We assume that $X_1,\ldots, X_n \iid \mathrm{Beta}(\alpha,\beta)$, where $\alpha,\beta\geq 1$, and our goal is to produce differentially private (DP) synthetic data. 	Often, to ensure finite sensitivity, the data are clamped to artificial bounds $[a,b]$, introducing bias in the DP estimate. Typically, these bounds are fixed in $n$, resulting in asymptotically negligible noise, but $O_p(1)$ bias. However, we show that it is possible to increase the bounds in $n$ to produce both noise and bias of order $o_p(n^{-1/2})$, resulting in an efficient DP estimator. We show through simulations that using this estimator along with Algorithm \ref{alg:onestep} results in a DP sample with optimal asymptotics. While we work with the beta distribution, this approach may be of value for other exponential family distributions as well.

		Recall that  $n^{-1} \sum_{i=1}^n \log(X_i)$ and $n^{-1}\sum_{i=1}^n \log(1-X_i)$ are sufficient statistics for the beta distribution. We will add Laplace noise to each of these statistics to achieve differential privacy. However, the sensitivity of these quantities is unbounded. First we pre-process the data by setting
		$\twid X_i = \min\{\max(X_i,t),1-t\},$
		where $t$ is a threshold that depends on $n$. Then the $\ell_1$-sensitivity of the pair of sufficient statistics is
		$\Delta(t) =2n^{-1}\left|\log(t)-\log(1-t)\right|$.
		We add independent noise to each of the statistics from the distribution $\mathrm{Laplace}\{\Delta(t)/\epsilon\}$, which results in $\ep$-DP versions of these statistics. Finally, we estimate $\theta=(\alpha,\beta)$ by plugging in the privatized sufficient statistics into the log-likelihood function and maximizing over $\theta$. The resulting parameter estimate satisfies $\ep$-DP by post-processing.

		We must carefully choose the threshold $t$ to ensure that the resulting estimate is efficient. The choice of $t$ must satisfy $\Delta(t) = o(n^{-1/2})$ to ensure that the noise does not affect the asymptotics of the  likelihood function. We also require that both $P(X_i<t)=o(n^{-1/2})$, and $P(X_i>1-t)=o(n^{-1/2})$ to ensure that $\twid X_i =X_i + o_p(n^{-1/2})$, which limits the bias to  $o_p(n^{-1/2})$. For the beta distribution, we can calculate that $P(X_i<t) = O(t^\alpha)$ and $P(X_i>1-t) = O(t^\beta)$. Since we assume that $\alpha,\beta\geq 1$, so long as $t=o(n^{-1/2})$ the probability bounds will hold. Taking $t = \min[1/2,10/\{\log(n) \sqrt n\}]$ satisfies $t=o(n^{-1/2})$, and we estimate the sensitivity as
		\[\Delta(t)\leq 2n^{-1} \log(t^{-1})\leq 2n^{-1}\log\{\log(n)\sqrt n\} = O\{\log(n)/n\} = o(n^{-1/2}),\]
		which satisfies our requirement for $\Delta$. While there are other choices of $t$ which would satisfy the requirements, our threshold was chosen to optimize the finite sample performance, so that the asymptotics could be illustrated with smaller sample sizes.

	\section{Proofs and technical lemmas}\label{s:proofs}

	For two distributions $P$ and $Q$ on $\RR^k$, the \emph{Kolmogorov-Smirnov distance} (KS-distance) is $\mathrm{KS}(P,Q) = \sup_{R\text{ rectangle}}|P(R)-Q(R)|$, where the supremum is over all axis-aligned rectangles. If $F_P$ and $F_Q$ are the multivariate cdfs of $P$ and $Q$, then $\lVert F_P-F_Q\rVert_\infty \leq \mathrm{KS}(P,Q)\leq 2^k\lVert F_P-F_Q\rVert_\infty$, so convergence in distribution is equivalent to convergence in KS-distance \citep{smith2011privacy}. By definition, we have that $\mathrm{TV}(P,Q)\geq \mathrm{KS}(P,Q)$.

	\thmTVnaive*
	\begin{proof}
	First we will establish the asymptotic distribution of $\thz$. Recall that by efficiency, we know that $\sqrt n \{\thx - \theta\} \overset d \rightarrow N\{0,I^{-1}(\theta)\}$ and $\sqrt n [\thz - \EE\{ \thz\mid \thx\}]\mid \thx \overset d \rightarrow N[0,I^{-1}\{\thx\}]$. Then by Slutsky's theorem, we have that $\sqrt n [\thz-\EE \{\thz\mid \thx\}] \overset d \rightarrow N\{0, I^{-1}(\theta)\}$. We can easily compute that $\cov[\thz - \EE\{\thz\mid \thx\}, \thx]=0$ using the law of total covariance. So, we have that $\sqrt n [\thz - \EE\{\thz\mid \thx\}+\thx-\theta] \overset d \rightarrow N\{0, 2I^{-1}(\theta)\}$. We also know that  $\sqrt n [\EE\{\thz\mid \thx\} - \thx] = o_p(1)$, since $\EE\{\thz\mid \thx\} = \thx+o_p(n^{-1/2})$. Together, we have that $\sqrt n \{\thz-\theta\} \overset d \rightarrow N\{0, 2I^{-1}(\theta)\}$.

\begin{align}
\mathrm{TV}\left(\ul X,\ul Z\right)&\geq \mathrm{TV}\left\{\sqrt n(\thx-\theta),\sqrt n (\thx-\theta)\right\}\label{eq:process}\\
&\geq \mathrm{KS}\left\{\sqrt n(\thx-\theta),\sqrt n (\thx-\theta)\right\}\label{eq:ks}\\
&\geq\mathrm{KS}\left[N\{0,I^{-1}(\theta)\},N\{0,2I^{-1}(\theta)\}\right]-\mathrm{KS}\left[\sqrt n (\thx-\theta),N\{0, I^{-1}(\theta)\}\right]\\
&\phantom{=} -\mathrm{KS}\left[\sqrt n(\thz-\theta),N\{0,2I^{-1}(\theta)\}\right]\label{eq:triangle}\\
&= \mathrm{KS}\left[N\{0,I^{-1}(\theta)\},N\{0,2I^{-1}(\theta)\}\right]+o(1)\label{eq:convD}\\
&\geq \Phi\left\{-\frac{\sqrt{\log(4)}}{\sqrt 2}\right\}-\Phi\left\{-\sqrt{\log(4)}\right\}+o(1)\label{eq:lower}\\
&\geq .083+o(1)

\end{align}
where \eqref{eq:process} is by the data processing inequality, \eqref{eq:ks} uses the KS-distance as a lower bound on total variation, \eqref{eq:triangle} applies two triangle inequalities since KS-distance is a metric, and \eqref{eq:convD} uses the asymptotic distributions of $\thx$ and $\thz$.

To establish \ref{eq:lower}, consider the following. Denote $\sigma^2 = (I^{-1}(\theta))_{1,1}$. Then consider the sequence of rectangles $R_i = \{x\in \RR^k \mid -(i+1)\sigma\leq x_1\leq -\sqrt{\log(4)}\sigma, \text{ and }-i\leq x_j\leq i, \forall j\neq 1\}$. Note that $R_i\subset R_{i+1}$ and $\bigcup_{i=1}^\infty R_i = \{x\in \RR\mid x_1\leq -\sqrt{\log(4)} \sigma\}$. Denote by $P$ the probability measure for $N\{0,I^{-1}(\theta)\}$ and $Q$ the probability measure for $N\{0,I^{-1}(\theta)\}$. Then

\begin{align*}
    \mathrm{KS}\left[N\{0,I^{-1}(\theta)\},N\{0,2I^{-1}(\theta)\}\right]
    &\geq \lim_{i\rightarrow \infty} |P(R_i)-Q(R_i)|\\
    &=\left|\lim_{i\rightarrow \infty} P(R_i)-\lim_{i\rightarrow\infty}Q(R_i)\right|\\
    &=\left|P\left(\bigcup_{i=1}^\infty R_i\right)-Q\left(\bigcup_{i=1}^\infty R_i\right)\right|\\
    &= \Phi\left\{-\frac{\sqrt{\log(4)}}{\sqrt 2}\right\}-\Phi\left\{-\sqrt{\log(4)}\right\}\\
&\geq .083.
\end{align*}

\end{proof}

For the following proofs, we will overload the $\frac{d}{d\theta}$ operator when working with multivariate derivatives. For a function $f:\RR^p\rightarrow \RR$, we write $\frac{d}{d\theta} f(\theta)$ to denote the $p\times 1$ vector of partial derivatives $(\frac{\partial}{\partial \theta_j} f(\theta))_{j=1}^p$. For a function $g:\RR^p\rightarrow \RR^q$, we write $\frac{d}{d\theta} g(\theta)$ to denote the $p\times q$ matrix $(\frac{\partial}{\partial \theta_j} g_k(\theta))_{j,k=1}^{p,q}$.

	Lemmas \ref{lem:score} and \ref{lem:fisher} are used for the proof of Theorem \ref{thm:onestep}. 		Parts 1 and 2 of Lemma \ref{lem:score} can be rephrased as the following: $\hat \theta$ is efficient if and only if it is consistent and $n^{-1}\sum_{i=1}^n S(\hat \theta,X_i) = o_p(n^{-1/2})$. The third property of Lemma \ref{lem:score} is similar to many standard expansions used in asymptotics, for example in \citet{van2000asymptotic}. However, we require the expansion for arbitrary efficient estimators, and include a proof for completeness.

	\begin{lemma}\label{lem:score}
		Suppose $X_1,\ldots, X_n\iid f_{\theta_0}$, and assume that (R1)-(R3) hold. Let $\hat \theta$ be an efficient estimator, which is a sequence of zeros of the score equations. Suppose that $\twid \theta$ is a $\sqrt n$-consistent estimator of $\theta_0$. Then
		\begin{enumerate}
			\item If $n^{-1} \sum_{i=1}^n S(\twid \theta, X_i) = o_p(n^{-1/2})$, then $\twid \theta - \hat \theta = o_p(n^{-1/2})$.
			\item If $\twid \theta$ is efficient, then $n^{-1}  \sum_{i=1}^n S(\twid \theta, X_i) = o_p(n^{-1/2})$.
			\item If $\twid \theta$ is efficient, then $\twid \theta = \theta_0 + I^{-1}(\theta_0) n^{-1} \sum_{i=1}^n S(\theta_0,X_i) + o_p(n^{-1/2})$.
		\end{enumerate}
	\end{lemma}
	\begin{proof}
		As $\twid \theta$ and $\hat \theta$ are both $\sqrt n$-consistent, we know that $\twid \theta -\hat \theta = O_p(n^{-1/2})$. So, we may consider a Taylor expansion of the score function about $\twid \theta = \hat \theta$.
		\begin{equation}\begin{split}\label{eq:efficient}
		n^{-1} \sum_{i=1}^n S(\twid \theta, X_i) &= n^{-1} \sum_{i=1}^n S(\hat \theta,X_i) + \left\{\frac{d}{d \hat\theta} n^{-1} \sum_{i=1}^n S( \hat\theta,X_i)\right\}(\twid \theta - \hat \theta) + O_p(n^{-1}) \\
    		&= 0+\left\{ \frac{d}{d\hat\theta} n^{-1} \sum_{i=1}^n S(\hat\theta, X_i) +  O_p(n^{-1/2})\right\}(\twid \theta - \hat \theta)\\
		&=\left\{-I(\theta_0) + o_p(1)\right\}(\twid \theta - \hat \theta),
		\end{split}\end{equation}
		where we used assumptions (R1)-(R3) to justify that 1) the second derivative is bounded in a neighborhood about $\theta_0$ (as both $\hat \theta$ and $\twid \theta$ converge to $\theta_0$), 2) the derivative of the score converges to $-I(\theta_0)$ by \citet[Theorem 7.2.1]{lehmann2004elements} along with the Law of Large Numbers, and 3) that $I(\theta_0)$ is finite, by (R3).

		To establish property 1, note that the left hand side of Equation \eqref{eq:efficient} is $o_p(n^{-1/2})$ implying that $\left(\twid \theta - \hat \theta\right) = o_p(n^{-1/2})$. For property 2, recall that by \citet[Page 479]{lehmann2004elements}, if $\twid \theta$ and $\hat \theta$ are both efficient, then $\left(\twid \theta - \hat \theta\right) = o_p(n^{-1/2})$. Plugging this into the right hand side of Equation \eqref{eq:efficient} gives $n^{-1} \sum_{i=1}^n S(\twid \theta, X_i) = o_p(n^{-1/2})$, establishing property 2.

		For property 3, we consider a slightly different expansion:
		\begin{align*}
		    o_p(n^{-1/2})& = n^{-1} \sum_{i=1}^n S(\twid \theta, X_i)\\
		&= n^{-1} \sum_{i=1}^n S(\theta_0,X_i) + \frac d{d\theta_0} n^{-1}\sum_{i=1}^n S(\theta_0,X_i) (\twid \theta - \theta_0) + O_p(n^{-1}),\\
		&=n^{-1} \sum_{i=1}^n S(\theta_0,X_i) +\{- I(\theta_0)+o_p(1)\}(\twid \theta - \theta_0) + O_p(n^{-1})
		\end{align*}
		where we used property 2 for the first equality, expanded the score about $\hat\theta = \theta_0$ for the second, and justify the $O_p(n^{-1})$ by (R2). By (R1)-(R2) and Law of Large Numbers along with \citet[Theorem 7.2.1]{lehmann2004elements}, we have the convergence of the derivative of score to $-I(\theta_0)$. By (R3), $I(\theta_0)$ is invertible. Solving the equation for $\twid \theta$ gives the desired result.
	\end{proof}

		\begin{lemma}\label{lem:fisher}
		Assume that (R0)-(R3) hold, and let $\omega_1,\ldots, \omega_n \iid P$. Then
		\begin{align*}
		n^{-1} \sum_{i=1}^n \frac{d}{d\theta}  S\{\theta,X_\theta(\omega_i)\}= o_p(1).
		\end{align*}
	\end{lemma}

	\begin{proof}
	First we can express the derivative as
	\[ n^{-1} \sum_{i=1}^n \frac{d}{d\theta}S\{\theta,X_\theta(\omega_i)\}
	=n^{-1} \sum_{i=1}^n \left[ \frac{d}{d\alpha} S\{\alpha,X_{\theta}(\omega_i)\} + \frac{d}{d\alpha}S\{\theta,X_\alpha(\omega_i)\}\right]\Big|_{\alpha=\theta}.\]
	The result follows from the Law of Large Numbers, provided that
	\[\EE_{\omega \sim P}\left[ \frac{d}{d\alpha} S\{\alpha,X_{\theta}(\omega)\} + \frac{d}{d\alpha}S\{\theta,X_\alpha(\omega)\}\right]\Big|_{\alpha=\theta}=0.\]  The expectation of the first term is $-I(\theta)$, by \citet[Theorem 7.2.1]{lehmann2004elements}. For the second term, we compute
		\begin{align*}
		\EE_{\omega \sim P} \frac{d}{d\alpha} S\{\theta, X_\alpha(\omega)\}\Big|_{\alpha=\theta}
		&=\int_\Omega \frac{d}{d\alpha} S\{\theta,X_\alpha(\omega)\} \Big|_{\alpha=\theta}\ dP(\omega)\\
		&= \int_{\mscr X} \frac{d}{d\alpha}S(\theta,x)  f_\alpha(x)\Big|_{\alpha=\theta}\ d\mu(x)\\
		&= \int_{\mscr X} S(\theta,x)\left\{ \frac{d}{d\alpha} f_\alpha(x)\Big|_{\alpha=\theta}\right\}^\top \ d\mu(x)\\
		&= \int_{\mscr X} S(\theta,x) \left\{\frac{\frac{d}{d\theta} f_\theta(x)}{f_\theta(x)}\right\}^\top f_\theta(x) \ d\mu(x)\\
		&=\int_{\mscr X} S(\theta,x)S^\top(\theta,x) f_\theta(x) \ d\mu(x) \\
		&=\EE_{X\sim \theta} \left\{S(\theta,X)S^\top(\theta,X)\right\}\\
		&=I(\theta).
		\end{align*}
	\end{proof}

\thmonestep*
\begin{proof}
We expand $\thz$ about $\thx$ using part 3 of Lemma \ref{lem:fisher}:
\begin{equation}\label{eq:score0}
    \thz=\thx + I^{-1}\{\thx\} n^{-1} \sum_{i=1}^n S\{\thx,X_{\thx}(\omega_i)\} + o_p(n^{-1/2})
    \end{equation}

   The score can be expanded about $\thx=\theta_0$:
\begin{align*}
&n^{-1} \sum_{i=1}^n S\{\thx,X_{\thx}(\omega_i)\}\\
&=n^{-1} \sum_{i=1}^n S\{\theta_0,X_{\theta_0}(\omega_i)\} + \left[\frac{d}{d\twid\theta} n^{-1} \sum_{i=1}^n S\{\twid\theta,X_{\theta^*}(\omega_i)\}\right]\{\thx-\theta_0\}\\
&=n^{-1} \sum_{i=1}^n S\{\theta_0,X_{\theta_0}(\omega_i)\} + o_p(1)O_p(n^{-1/2}),
\end{align*}
where $\twid\theta$ is between $\thx$ and $\theta_0$; by Lemma \ref{lem:fisher}, we justify that the derivative is $o_p(1)$.

Combining this derivation along with the fact that $I^{-1}\{\thx\}= I^{-1}(\theta_0)+o_p(1)$ by the continuous mapping theorem, we have the following equation:

  \begin{equation}\label{eq:score1}
    \thz =\thx + I^{-1}(\theta_0) n^{-1} \sum_{i=1}^n S\{\theta_0,X_{\theta_0}(\omega_i)\} + o_p(n^{-1/2}).
\end{equation}

Using the same techniques, we do an expansion for $\thy$ about $\theta^*=2\thx-\thz$:
\begin{align}
    \thy &= \theta^* + I^{-1}(\theta^*) n^{-1} \sum_{i=1}^n S\{\theta^*, X_{\theta^*} (\omega_i)\} + o_p(n^{-1/2})\\
    &= \theta^* + I^{-1}(\theta_0) n^{-1} \sum_{i=1}^n S\{\theta_0, X_{\theta_0}(\omega_i)\} + o_p(n^{-1/2})\label{eq:substitute}\\
    &=\theta^* + \{\thz-\thx\} + o_p(n^{-1/2})\label{eq:combine}\\
    &=\thx+o_p(n^{-1/2}),\label{eq:final}\\
\nonumber
\end{align}
where line \eqref{eq:substitute} is a similar expansion as used for equation \eqref{eq:score0}, in line \eqref{eq:combine} we substituted the expression from \eqref{eq:score1}, and line \eqref{eq:final}
uses the fact that as $n\rightarrow \infty$, $\theta^*=2\thx-\thz$ with probability tending to one. Indeed, since $2\thx-\thz$ is a consistent estimator of $\theta_0$, we have that as $n\rightarrow \infty$, $ P\{2\thx-\thz\in \Theta\}\geq P\{2\thx-\thz \in B(\theta_0)\}\rightarrow 1$.
\end{proof}

\lemconditional*
\begin{proof}
For a fixed $\theta \in \Theta$, for $\omega\sim P$, the random variable $Y = X_{\theta}(\omega)$ is distributed with probability measure $P X_{\theta}^{-1}$: for any measurable set $E$, $P(Y\in E) = P X_{\theta}^{-1}(E)$. We denote by $P^n_\Omega$ the joint probability measure on $\Omega^n$, and $(P X_{\theta}^{-1})^n$ the joint probability measure on $\mscr X^n$.

Given $\theta^* \in \Theta$, our goal is to derive the probability distribution of the random variables $X_{\theta^*}(\omega_1),\ldots,$ $X_{\theta^*}(\omega_n)$ conditioned on

the event that $\{\omega_1,\ldots,\omega_n\mid \hat\theta\{X_{\theta^*}(\omega_i)\}= \hat\theta\}$.
However, this event may have zero probability.  Instead, we will condition on $S_{\hat \theta,\theta^*}^\delta  = \{\omega_1,\ldots, \omega_n\mid \hat \theta\{X_{\theta^*}(\omega)\} \in B_\delta(\hat \theta)\}$, where $B_\delta(\hat \theta) = \{\theta  \mid \lVert \hat\theta - \theta\rVert\leq \delta\}$, which has positive probability. At the end, we will take the limit as $\delta\rightarrow 0$ to derive the desired distribution.

 Let $E\subset \mscr X^n$ be a measurable set. Then
 \begin{align*}
     &P\{X_{\theta^*}(\omega_1),\ldots, X_{\theta^*}(\omega_n) \in E \mid \omega_1,\ldots, \omega_n \in S_{\theta^*,\hat\theta}^\delta\}\\
     &=P(\omega_1,\ldots, \omega_n \in X^{-1}_{\theta^*} E\mid \omega_1,\ldots, \omega_n \in S^{\delta}_{\theta^*,\hat\theta})\\
     &= \frac{P^n( X^{-1}_{\theta^*}E \cap S_{\theta^*,\hat\theta}^\delta)}{P^n (S_{\theta^*,\hat\theta}^\delta)}\\
     &= \frac{(P X_{\theta^*}^{-1})^n (E\cap X_{\theta^*} S^{\delta}_{\theta^*,\hat\theta})}{(P X_{\theta^*}^{-1})^n(X_{\theta^*} S_{\theta^*,\hat\theta}^\delta)},
 \end{align*}
 where we used the definition of conditional probability and the fact that $X_{\theta^*}^{-1} X_{\theta^*} S_{\theta^*,\hat\theta}^\delta = S_{\theta^*,\hat\theta}^\delta$.

 This last expression shows that $X_{\theta^*}(\omega_1),\ldots, X_{\theta^*}(\omega_n)$ conditioned on $\ul \omega \in S_{\theta^*,\hat\theta}^\delta$ is distributed as $f^n_{\theta^*}\{y_1,\ldots, y_n \mid \hat\theta(\ul y) \in B_\delta(\hat\theta)\}$.
  This derivation is valid for all $\delta>0$. Taking the limit as $\delta\rightarrow 0$ gives the desired formula:
  \[Y^{\theta^*}_1,\ldots, Y^{\theta^*}_{n} \Big | \hat\theta(\ul Y^{\theta^*})=\hat\theta(\ul X) \sim f_{\theta^*}^n\{y_1,\ldots, y_n \mid \hat\theta(\ul y) = \hat\theta(\ul X)\}.\]

\end{proof}

  \thmklGeneral*
  \begin{proof}
  We will consider for the purposes of this proof that $\thx$ is a randomized statistic, so we can also write $g(\hat\theta\mid \ul X)$ to represent the distribution of $\thx$ given $\ul X$ (which we assume does not depend on $\theta$, since   $\thx$ is a statistic). Any deterministic statistic can be expressed as a limit of randomized statistics, where the noise due to randomness goes to zero. For example, \citet{barber2020testing} consider statistics which are solutions to noisy score equations, where the noise is normally distributed and $o_p(n^{-1/2})$. While the distributions $g$ depend on $n$, we will suppress this dependence for notational simplicity. We can then express the desired KL divergence as follows:

    First, by the data processing inequality, we can add in the random variable $\hat\theta(\ul X) = \hat\theta(\ul Y)$ to get an upper bound on the KL divergence. We then have closed formulas for the joint distributions $\{X_1,\ldots,X_n,\hat\theta(\ul X)\}$ and $\{Y_1,\ldots, Y_n, \hat\theta(\ul X)\}$.
    \begingroup
\allowdisplaybreaks
  \begin{align}
  \mathrm{KL}&\left\{X_1,\ldots,X_n\middle|\middle| Y_1,\ldots, Y_n\right\}\\
  &\leq
  \mathrm{KL}\left\{X_1,\ldots,X_n,\hat\theta(\ul X)\middle|\middle| Y_1,\ldots, Y_n,\hat\theta(\ul X)\right\}\\
  &=\mathrm{KL}\left[f_\theta^n\{\ul x\mid \hat\theta(\ul x)\}g_{\theta}\{\hat\theta(x)\} \middle|\middle| f_{\theta_n}^n\{\ul x\mid \hat\theta(\ul x)\}g_{\theta}\{\hat\theta(\ul x)\}\right]\\

    &=\EE_{\hat\theta \sim g(\cdot\mid \ul X)}\EE_{\ul X\sim f_\theta} \log \left\{ \frac{f_\theta(\ul X\mid\hat\theta) g_\theta(\hat\theta) }
    {f_{\theta_n}(\ul X\mid \hat \theta) g_\theta(\hat\theta) } \right\}\label{eq:kl2}\\
    &=\EE_{\hat\theta \sim g(\cdot\mid \ul X)}\EE_{\ul X\sim f_\theta} \log \left\{ \frac{f_\theta (\ul X)g(\hat\theta\mid \ul X) }
    {f_{\theta_n}(\ul X\mid \hat \theta) g_\theta(\hat \theta)}\right\}\label{eq:kl3}\\
    &=\EE_{\hat\theta \sim g(\cdot\mid \ul X)}\EE_{\ul X\sim f_\theta} \log \left[\frac{f_\theta(\ul X) g(\hat\theta\mid \ul X)}{\{f_{\theta_n}(\ul X) g(\hat\theta\mid \ul X)/g_{\theta_n}(\hat\theta)\}  g_\theta(\hat\theta)}\right]\label{eq:kl4}\\
         &=\EE_{\hat\theta \sim g(\cdot\mid \ul X)}\EE_{\ul X\sim f_\theta} \log \left\{\frac{f_{\theta}(\ul X) {g(\hat\theta\mid \ul X)}}{f_{\theta_n}(\ul X){g(\hat\theta\mid \ul X)}}\right\}
      + \EE_{\hat\theta \sim g(\cdot\mid \ul X)}\EE_{\ul X \sim f_{\theta}} \log \left\{\frac{g_{\theta_n}(\hat\theta)}{g_{\theta}(\hat\theta)}\right\}\label{eq:kl5}\\
      &=-\EE_{\ul X\sim f_\theta} \log \left(\frac{f_{\theta_n}(\ul X) }{f_{\theta}(\ul X)}\right) + \EE_{\hat\theta \sim g_\theta} \log \left(\frac{g_{\theta_n}(\hat\theta)}{g_{\theta}(\hat\theta)}\right)\label{eq:kl6},
  \end{align}
  \endgroup

 Above, line \eqref{eq:kl2} simply applies the definition of KL divergence, and line \eqref{eq:kl4} uses the definition of conditional distribution.

At this point, we need to compute the two expectations of line \eqref{eq:kl6}, and show that everything cancels except for an $o_p(1)$ term.

We write $\ell(\theta\mid \ul x) = \sum_{i=1}^n \log f_\theta(x_i)$. Using our assumptions, we can expand $\ell(\theta_n\mid \ul x)$:
\begin{align*}
    \ell(\theta_n\mid \ul x) &= \ell(\theta\mid \ul x) + (\theta_n - \theta)^\top \nabla \ell(\theta\mid \ul x) + \frac 12 (\theta_n - \theta)^\top \nabla^2 \ell(\theta\mid \ul x) (\theta_n- \theta)\\
    &\phantom{=}+\frac{1}{6} \xi^*\sum_{i,j,k} (\theta_n-\theta)_i(\theta_n-\theta)_j(\theta_n-\theta)_k \sum_{s=1}^ng_{ijk}(x_s),
    \end{align*}
    where $|\xi^*|<1$ and $g_{ijk}(x)$ is an upper bound for $\left|\frac{\partial^3 \ell(\theta\mid \ul x)}{\partial \theta_i\theta_j\theta_k}\right|$ for a ball $\theta \in N(\theta)$ which exists by (R3). These expansions are based on those from \citet{Serfling:1980}.
Applying $\EE_{\ul X \sim f_\theta}$ to this derivation gives
\begin{equation}\label{eq:klF}
\begin{split}
    \EE_{\ul X\sim f_\theta} \log \left\{ \frac{f_{\theta_n}(\ul X)}{f_\theta(\ul X)}\right\}
    &= 0-\frac {n}{2} (\theta_n - \theta)^\top I(\theta)(\theta_n - \theta) \\
    &\phantom{=}+ O(1)\frac{n}{6} \sum_{i,j,k} \{\EE g_{i,j,k}(x)\} (\theta_n-\theta)_i(\theta_n- \theta)_j (\theta_n - \theta)_k,\\
    &= \frac{-n}{2}  (\theta_n - \theta)^\top I(\theta)(\theta_n - \theta) + O(n) \lVert \theta_n-\theta\rVert^3
\end{split}
\end{equation}
where the first term is zero as the expected value of the score function is zero by (R3),
the second term uses \citet[Theorem 7.2.1]{lehmann2004elements} and (R3). The $O(1)$ factor in the third term is based on the fact that $|\xi^*|\leq 1$. Finally, note that $\sum_{i,j,k}\{\EE g_{i,j,k}(x)\}(\theta_n-\theta)_i(\theta_n-\theta)_j(\theta_n-\theta)_k\leq p^3 \sup_{i,j,k} \{\EE g_{i,j,k}(x)\} \lVert \theta_n-\theta\rVert^3_\infty = O(1) \lVert \theta_n-\theta\rVert^3$. Note that all norms are equivalent in $\RR^p$, so they can be interchanged up to a factor of $O(1)$.

Next, we will derive a similar formula for $\log g_{\theta^*}(\hat \theta)$:
\begin{equation}
\begin{split}
    \log g_{\theta_n}(\hat\theta)
    &=\log g_{\theta}(\hat\theta)
    +\nabla \log g_\theta(\hat\theta)(\theta_n - \theta)
    +\frac12 (\theta_n- \theta)^\top \nabla^2 \log g_\theta(\hat\theta) (\theta_n - \theta)\\
    &\phantom{=}+\frac n6\xi_2^* \sum_{i,j,k} (\theta_n-\theta)_i(\theta_n-\theta)_j(\theta_n-\theta)_k  G_{i,j,k}(\hat\theta),
    \end{split}
    \label{eq:expansionG}
\end{equation}
where $|\xi_2^*|\leq 1$. In order to apply the expectation $\EE_{\hat\theta \sim \theta}$ to this equation, we will first show $\EE_{\hat\theta\sim \theta} \nabla \log g_\theta(\hat\theta)=0$ and $\EE_{\hat\theta\sim \theta} \nabla^2 \log g_\theta(\hat\theta)=-nI(\theta)+o(n)$.

\begin{align*}
    \left\{\EE_{\hat\theta \sim \theta} \nabla \log g_\theta(\hat\theta)\right\}_j
    &= \int \left\{\frac{\partial}{\partial \theta_j} \log g_\theta(\hat\theta)\right\} g_\theta(\hat\theta)\ d\hat\theta\\
    &= \int \frac{\partial}{\partial \theta_j} g_\theta(\hat\theta)\ d\hat\theta\\
    &=\int \frac{\partial}{\partial \theta_j} \int_x f_\theta(x) g(\hat\theta\mid x) \ dx\ d\hat\theta\\
    &= \int_{\hat\theta} \int_x \frac{\partial}{\partial \theta_j} f_\theta(x)  g(\hat\theta\mid x) \ dx\ d\hat\theta\\
    &=\frac{\partial}{\partial \theta_j} \int\int f_\theta(x) g(\hat\theta\mid x) \ dx\ d\hat\theta\\
    &=0,
\end{align*}
where we use the assumption (R3) that $\left|\frac{\partial}{\partial \theta_j}f_\theta(x)\right|$ is  bounded above by an integrable function and the dominated convergence theorem to interchange the derivative and the integral.

Next we work on the second derivative:

\begin{align*}
    \left\{\EE_{\hat\theta \sim \theta} \nabla^2 \log g_\theta(\hat\theta)\right\}_{j,k}
    &= \int \left\{\frac{\partial^2}{\partial \theta_j\partial \theta_k} \log g_\theta(\hat\theta)\right\} g_\theta(\hat\theta) \ d\hat\theta\\
    &=\int \frac{g_\theta(\hat\theta) \frac{\partial^2}{\partial \theta_j\partial\theta_k} g_\theta(\hat\theta) - \frac{\partial}{\partial \theta_j} g_\theta(\hat\theta) \left\{\frac{\partial}{\partial \theta_k} g_\theta(\hat\theta)\right\}^\top }{g_\theta^2(\hat\theta)} g_\theta(\hat\theta) \ d\hat\theta\\
    &=0-\EE_{\hat\theta\sim \theta}\left\{ \nabla \log g_\theta(\hat\theta) \nabla^\top \log g_\theta(\hat\theta)\right\}_{j,k},\\
\end{align*}
where we used (R3) along with the dominated convergence theorem to set the first term equal to zero. We see that $\EE \nabla^2 \log g_\theta(\hat\theta) = -I_{\thx}(\theta)$, where $I_{\thx}$ represents the Fisher information of the random variable $\hat\theta(\ul X)\sim g$. It is our current goal to show that $I_{\thx}(\theta) = nI(\theta)+o(n)$, where $I(\theta)$ is the Fisher information for one sample $X\sim f_\theta$. First note that by the data processing inequality \citep{zamir1998proof}, $I_{\hat\theta(\ul X)}(\theta)\leq I_{X_1,\ldots, X_n}(\theta) = nI(\theta)$, where the inequality represents the positive-definite ordering of matrices. Next, we need to find a matching lower bound. By the Cramer Rao lower bound, we have that
\[ \{I_{\hat\theta(\ul X)}(\theta)\}^{-1} \leq \var\{\hat\theta(\ul X)\}+o(1/n),\]
where $\var\{\thx\}$ is the covariance matrix of the random variable $\thx$, and we used the fact that $\thx$ is asymptotically unbiased. By the efficiency of $\hat\theta(\ul X)$, we have that
\[\var\{\hat\theta(\ul X)\} = n^{-1} I^{-1}(\theta) + o(1/n).\]
We then have
\begin{align*}
    I_{\hat\theta(\ul X)}(\theta) &\geq \left\{n^{-1} I^{-1}(\theta) + o(1/n)\right\}^{-1}
    = n\{I^{-1}(\theta) + o(1)\}^{-1}
    =n \{I(\theta) + o(1)\},
\end{align*}
where for the last equality, we use the following matrix identity:
\[(A+B)^{-1} = A^{-1} - A^{-1} B(A+B)^{-1},\]
where we set $A = I^{-1}(\theta)$ and $B=o(1)$.

Combining our results, we have that

\[\EE_{\hat\theta \sim \theta} \nabla^2 \log g_\theta(\hat\theta) = -I_{\hat\theta(\ul X)}(\theta) = n\{-I(\theta) + o(1)\}.\]
Finally, applying the expectation to equation \eqref{eq:expansionG}, we have
\begin{equation}\label{eq:g}
\begin{split}
\EE_{\hat\theta \sim \theta}\log \left\{\frac{g_{\theta_n}(\hat\theta)}{g_\theta(\hat\theta)}\right\}
&=0 - \frac{n}{2} (\theta_n - \theta) \{I(\theta) + o(1)\}(\theta_n-\theta)\\
&\phantom{=}+ O(1) \frac n6 \sum_{i,j,k} \{\EE G_{i,j,k}(\hat\theta)\} (\theta_n-\theta)_i(\theta_n-\theta)_j(\theta_n-\theta)_k\\
&=\frac{-n}{2}(\theta_n-\theta)^\top I(\theta)(\theta_n-\theta) + o(n)\lVert \theta_n-\theta\rVert^2 + O(n) \lVert \theta_n-\theta\rVert^3.
\end{split}
\end{equation}
Combining equations \eqref{eq:klF} and \eqref{eq:g}, we have
\begin{align*}
\mathrm{KL}&\left(X_1,\ldots,X_n\middle|\middle| Y_1,\ldots, Y_n\right)\\
&\leq\mathrm{KL}\left[f_\theta^n\{\ul x\mid \hat\theta(\ul x)\}g_{\theta}\{\hat\theta(x)\} \middle|\middle| f_{\theta_n}^n\{\ul x\mid \hat\theta(\ul x)\}g_{\theta}\{\hat\theta(\ul x)\}\right]\\

    &=-\EE_{\ul X\sim f_\theta} \log \left\{\frac{f_{\theta_n}(\ul X) }{f_{\theta}(\ul X)}\right\}
      + \EE_{\hat\theta \sim g_\theta(\hat\theta)} \log \left\{\frac{g_{\theta_n}(\hat\theta)}{g_{\theta}(\hat\theta)}\right\}\\
      &=o(n)\lVert \theta_n-\theta\rVert^2 + O(n) \lVert \theta_n-\theta\rVert^3.
\end{align*}

\end{proof}

\begin{comment}
\begin{proof}[Lemma 2]
First note that for all $\theta\in \Theta$, we have that $\hat\theta\{X_{\theta}(\ul \omega)\} = \theta + O_p(n^{-1/2})$, by the efficiency of $\hat\theta$. In particular, we have that $\hat\theta\{X_{\theta^*}(\ul \omega)\} = \theta^* + O_p(n^{-1/2})$.
We also have that $\hat\theta(X_{\theta^*}(\ul \omega)) = \hat\theta(X) = \theta + O_p(n^{-1/2})$. Combining these facts, we have that $\theta^* = \theta + O_p(n^{-1/2})$.
\end{proof}

\end{comment}

\bibliography{bibliography.bib}

\end{document}

config: arxiv