Measurement error models: from nonparametric methods to deep neural networks

We use a fully connected feed-forward neural network (FNN) to model the regression function $f$. With $\theta$ denoting parameters of this FNN, we can write $f=f_{\theta}$. Next, we use a normalizing flow (Tabak and Turner,, 2013) to represent the prior distribution $p_{X}$. A normalizing flow is a sequence of transformations $g_{1},g_{2},\ldots,g_{m}$, where each $g_{j}$ is an invertible mapping from $\mathbb{R}^{d}$ to $\mathbb{R}^{d}$. Let $V\in\mathbb{R}^{d}$ be a random vector whose density has a simple analytic form (e.g., $V\sim{\cal N}(0,I_{d})$). We model the density of $X$ by assuming $X=g^{-1}(V)$, where $g=g_{m}\circ g_{m-1}\circ\ldots\circ g_{1}$. Let $J(v)$ denote the Jacobian of the mapping $g$. The change of variables formula implies that

$$p_{X}(x)=p_{V}(v)\cdot|\det J(v)|\Big{|}_{v=g(x)}.$$

By parametrizing each $g_{j}$ via neural network layers, this idea has been applied to density estimation (Papamakarios et al.,, 2017) and variational inference (Rezende and Mohamed,, 2015). We adopt the specific flow proposed by Dinh et al., (2015), called non-linear independent components estimation (NICE). NICE parametrizes each invertible transformation $x=g_{j}(v)$ as follows: $x_{I_{1}}=v_{I_{1}}$ and $x_{I_{2}}=v_{I_{2}}+h(v_{I_{1}})$, where $I_{1}\cup I_{2}$ is a partition of $\{1,2,\ldots,d\}$ and $h$ is a neural network with $|I_{1}|$ input units and $|I_{2}|$ output units. Such a mapping has a trivial inverse: $v_{I_{1}}=x_{I_{1}}$ and $v_{I_{2}}=x_{I_{2}}-h(x_{I_{1}})$. Furthermore, its Jacobian always has a unit determinant, regardless of the choice of $h$. It thus follows that

$$p_{X}(x)=p_{V}(g_{\gamma}(x)),\qquad\mbox{where $\gamma$ denotes parameters of NICE}.$$

(3)

Combining (3) with the MEM in (1), we can write down the joint density of $(W,Y,X)$:

$$p(w,y,x;\theta,\gamma,\sigma^{2})=p_{V}(g_{\gamma}(x))\cdot p_{U}(w-x)\cdot p_{\epsilon}(y-f_{\theta}(x);\sigma^{2}).$$

(4)

Here, $p_{U}(\cdot)$ is the density of measurement errors, which is known; $p_{V}(\cdot)$ is a simple analytical density, which is typically chosen as the $d$-dimensional standard Gaussian; and $p_{\epsilon}(\cdot)$ is the normal density of noise in $y$. For a moment, we assume $\sigma^{2}$ is known and write this joint density as $p_{\theta,\gamma}(w,y,x)$. Later, in our main algorithm, $\sigma^{2}$ will be estimated jointly with other parameters.

Now that the MEM is parametrized by $(\theta,\gamma)$, where $\theta$ contains parameters of interest and $\gamma$ contains nuisance parameters, we estimate $(\theta,\gamma)$ by maximizing the marginal log-likelihood of the observed variables:

$$L^{*}(\theta,\gamma|w,y)\equiv\log\int p_{\theta,\gamma}(w,y,x)dx.$$

(5)

This objective function involves an integral over $x$ and cannot be evaluated analytically.

The variational inference is a common tool for parameter estimation and inference in such latent variable models (Jordan et al.,, 1999). The variational auto-encoder (VAE) (Kingma and Welling,, 2014; Rezende et al.,, 2014) is a variational inference approach via neural networks. Standard variational inference approaches approximate the intractable posterior distribution by a parametric family whose densities have simple analytical forms, whereas VAE achieves the same goal by an inference network (also called a recognition model), where the posterior distribution is approximated by a neural network transformation of multivariate normal densities. It has been widely recognized that VAE is capable of approximating complex posterior distributions. Burda et al., (2016) proposed an improvement of VAE by incorporating the idea of importance sampling, namely the importance weighted autoencoder (IWAE). We adapt IWAE to maximize the marginal log-likelihood in (5).

Let $x_{1},x_{2},\ldots,x_{K}$ be IID samples drawn from a proposal distribution $q(x|w,y)$. According to Jensen’s inequality, the marginal log-likelihood has a lower bound:

	$\displaystyle L^{*}(\theta,\gamma|w,y)$	$\displaystyle=\log\left(\mathbb{E}_{x_{1:K}\sim q(\cdot|w,y)}\left[\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\theta,\gamma}(w,y,x)}{q(x_{k}|w,y)}\right]\right)$
		$\displaystyle\geq\mathbb{E}_{x_{1:K}\sim q(\cdot|w,y)}\left[\log\left(\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\theta}(w,y,x_{k})}{q(x_{k}|w,y)}\right)\right]\;\equiv L(\theta,\gamma|w,y).$

This lower bound $L(\theta,\gamma|w,y)$ is called an evidence lower bound (ELBO). It is tighter if $q(x|w,y)$ is closer to the true posterior distribution of $X$ and/or $K$ is larger.

The idea of VAE is to use a neural network to model the proposal distribution and optimize it jointly with other parameters. However, VAE only corresponds to $K=1$. One can use the idea of importance sampling (e.g., see Liu, (2008)) with an appropriate sample size $K$ to tighten the lower bound. IWAE of Burda et al., (2016) combines these two ideas. We adopt IWAE and model $q(x|w,y)$ as a multivariate normal distribution ${\cal N}\bigl{(}\mu_{\phi}(w,y),\Sigma_{\phi}(w,y)\bigr{)}$, where the vector $\mu_{\phi}(w,y)\in\mathbb{R}^{d}$ and the diagonal matrix $\Sigma_{\phi}(w,y)\in\mathbb{R}^{d\times d}$ are both determined by another FNN whose parameters are denoted by $\phi$. This is equivalent to drawing samples $x_{k}$ from

$$x(\phi,z_{k})=\mu_{\phi}(w,y)+[\Sigma_{\phi}(w,y)]^{1/2}z_{k},\qquad\mbox{where}\quad z_{1},z_{2},\ldots,z_{K}\overset{iid}{\sim}{\cal N}(0,I_{d}).$$

(6)

Let $q_{\phi}(x|w,y)$ denote the density of ${\cal N}\bigl{(}\mu_{\phi}(w,y),\Sigma_{\phi}(w,y)\bigr{)}$ and let $x(\phi,z_{k})$ be the mapping in (6). We consider the following objective:

$$Q(\theta,\gamma,\phi|w,y)\equiv\mathbb{E}_{z_{1:K}\sim{\cal N}(0,I_{d})}\left[\log\left(\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\theta,\gamma}\big{(}w,y,x(\phi,z_{k})\big{)}}{q_{\phi}\big{(}x(\phi,z_{k})|w,y\big{)}}\right)\right].$$

(7)

The estimator is

$$\bigl{(}\hat{\theta}(w,y),\;{\hat{\gamma}(w,y)},\;\hat{\phi}(w,y)\bigr{)}=\mathrm{argmax}_{\theta,\gamma,\phi}Q(\theta,\gamma,\phi|w,y).$$

(8)

The neural network structure that implements this estimator is shown in Figure 1. The FNN parametrized by $\phi$ that generates samples from the proposal distribution is called “encoder.” The FNN parametrized by $\theta$ that approximates the regression function $f(x)$ and the normalizing flow parametrized by $\gamma$ that approximates the prior distribution of $X$ are both called “decoders.” The “encoder” takes data $(W,Y)$ as input, and outputs parameters $(\mu_{\phi},\Sigma_{\phi})$ for the proposal distribution. The “decoders” take as input the random samples generated from the proposal distribution, and outputs the estimated regression function $f_{\theta}(x)$ and density $p_{\gamma}(x)$ at those random samples of $X$.

We use a gradient ascent algorithm to solve (8). Given initial values $(\theta^{(0)},\gamma^{(0)},\phi^{(0)})$, the update at iteration $t$ is

$$(\theta^{(t+1)},\gamma^{(t+1)},\phi^{(t+1)})=(\theta^{(t)},\gamma^{(t)},\phi^{(t)})-\alpha_{t}\cdot\nabla Q(\theta^{(t)},\gamma^{(t)},\phi^{(t)}|w,y),$$

(9)

where $\alpha_{t}>0$ is the step size at iteration $t$, often called the learning rate. We use the adaptive learning rate suggested by Kingma and Ba, (2014), which is $\alpha_{t}=\alpha_{0}\sqrt{1-\alpha_{2}^{t}}/(1-\alpha_{1}^{t})$, for some default parameters $(\alpha_{0},\alpha_{1},\alpha_{2})$.

It remains to compute the gradient. Write for short $\beta_{k}=\frac{p_{\theta,\gamma}(w,y,x(\phi,z_{k}))}{q_{\phi}(x(\phi,z_{k})|w,y)}$, $1\leq k\leq K$. By direct calculations, we have

$$\nabla_{\theta,\gamma,\phi}Q(\theta,\gamma,\phi|w,y)=\mathbb{E}_{z_{1:K}\sim{\cal N}(0,I_{d})}\left[\sum_{k=1}^{K}\frac{\beta_{k}}{\sum_{\ell=1}^{K}\beta_{\ell}}\nabla_{\theta,\gamma,\phi}\log\left(\frac{p_{\theta,\gamma}(w,y,x(\phi,z_{k}))}{q_{\phi}(x(\phi,z_{k})|w,y)}\right)\right].$$

(10)

The Monte Carlo method approximates the gradients by replacing the expectation in (10) by realized values at a random sample of $z_{1:K}$. At each iteration, we draw $z_{1},z_{2},\ldots,z_{K}$ IID from ${\cal N}(0,I_{d})$ and estimate the gradient (10) by

$$\widehat{\nabla_{\theta,\gamma,\phi}Q}(\theta,\gamma,\phi|w,y)=\sum_{k=1}^{K}\frac{\beta_{k}}{\sum_{\ell=1}^{K}\beta_{\ell}}\nabla_{\theta,\gamma,\phi}\log\left(\frac{p_{\theta,\gamma}(w,y,x(\phi,z_{k}))}{q_{\phi}(x(\phi,z_{k})|w,y)}\right).$$

(11)

This plain gradient estimator works for approximating $\nabla_{\theta}Q$ and $\nabla_{\gamma}Q$, but not $\nabla_{\phi}Q$, where we recall that $\phi$ contains parameters in the encoder. The issue of $\widehat{\nabla_{\phi}Q}$ happens when $K\to\infty$. By theory of IWAE, the objective in (7), as a lower bound of the true marginal log-likelihood, converges to the true marginal log-likelihood as $K\to\infty$ (Burda et al.,, 2016). Therefore, a larger $K$ (i.e., more “importance samples”) should be favored. Unfortunately, if (7) is solved by a plain gradient ascent algorithm, the performance is often worse with an increased $K$. The reason is that $\widehat{\nabla_{\phi}Q}$ is too “noisy.” The true gradient $\nabla_{\phi}Q$ vanishes as $K\to\infty$; therefore, it is the bias $(\mathbb{E}_{z_{1:K}}[\widehat{\nabla_{\phi}Q}]-\nabla_{\phi}Q)$, not the true gradient $\nabla_{\phi}Q$, that contains the signal. Since the bias diminishes faster than the standard deviation, as $K$ increases, the signal-to-noise ratio in the gradient actually decreases. This is, however, not a problem for the gradient with respect to $(\theta,\gamma)$, as the true gradient $\nabla_{\theta,\gamma}Q$ does not vanish as $K\to\infty$. See Rainforth et al., (2018) for a rigorous theoretical justification.

Tucker et al., (2018) proposed a method to resolve the above issue. They derived an alternative expression of the gradient with respect to $\phi$:

		$\displaystyle\nabla_{\phi}Q(\theta,\gamma,\phi|w,y)$		(12)
	$\displaystyle=\;$	$\displaystyle\mathbb{E}_{z_{1:K}\sim{\cal N}(0,I_{d})}\left[\sum_{k=1}^{K}\biggl{(}\frac{\beta_{k}}{\sum_{\ell=1}^{K}\beta_{\ell}}\biggr{)}^{2}\cdot\frac{\partial}{\partial x_{k}}\log\left(\frac{p_{\theta,\gamma}(w,y,x_{k})}{q_{\phi}(x_{k}|w,y)}\right)\bigg{|}_{x_{k}=x(\phi,z_{k})}\cdot\frac{\partial}{\partial\phi}x(\phi,z_{k})\right],$		(13)

where $\beta_{k}$ is the same as in (10). We note that (12) cannot be derived from (10) via a simple chain rule. In fact, a direct use of the chain rule should yield a term related to $\frac{\partial}{\partial\phi}q_{\phi}(x_{k}|w,y)$, but there is no such term in (12). This alternative gradient formula motivates another Monte Carlo estimator for $\nabla_{\phi}Q$:

$$\widetilde{\nabla_{\phi}Q}(\theta,\phi|w,y)=\sum_{k=1}^{K}\biggl{(}\frac{\beta_{k}}{\sum_{\ell=1}^{K}\beta_{\ell}}\biggr{)}^{2}\cdot\frac{\partial}{\partial x_{k}}\log\left(\frac{p_{\theta,\gamma}(w,y,x_{k})}{q_{\phi}(x_{k}|w,y)}\right)\bigg{|}_{x_{k}=x(\phi,z_{k})}\cdot\frac{\partial}{\partial\phi}x(\phi,z_{k}),$$

(14)

which is called the doubly reparametrized gradient estimator. Compared with the plain gradient estimator (11), the variance of this estimator converges to zero at a faster rate, so that the signal-to-noise ratio still increases with $K$ (see Tucker et al., (2018) for a detailed explanation). The difference between (11) and (14) is analogous to the difference between the score function gradient estimator (or reinforce gradient estimator) and the reparametrization gradient estimator in Monte Carlo variational inference (Kingma and Welling,, 2014; Titsias and Lázaro-Gredilla,, 2014; Rezende et al.,, 2014). It is well-known that the reparametrization gradient estimator has a smaller variance.

Our main algorithm, called Neural Network for Measurement Error (NNME), is summarized in Algorithm 1.

We now provide a more detailed explanation of NNME. Note that the calculations in Sections 2.1-2.2 are for $n=1$. For a general $n$, write $w^{(n)}=\{w_{i}\}_{i=1}^{n}$ and $y^{(n)}=\{y_{i}\}_{i=1}^{n}$. At each epoch, we draw samples $\{z_{ik}\}_{1\leq i\leq n,1\leq k\leq K}$ IID from ${\cal N}(0,I_{d})$ and obtain $x_{ik}=x(\phi,z_{ik})$. Similar to (11) and (14), we have

	$\displaystyle\widehat{\nabla_{\theta,\gamma}Q}(\theta,\gamma,\phi|w^{(n)},y^{(n)})$	$\displaystyle=\sum_{i=1}^{n}\sum_{k=1}^{K}\frac{\beta_{ik}}{\sum_{\ell=1}^{K}\beta_{i\ell}}\nabla_{\theta,\gamma}\log\left(\frac{p_{\theta,\gamma}(w_{i},y_{i},x(\phi,z_{ik}))}{q_{\phi}(x(\phi,z_{ik})|w_{i},y_{i})}\right),$		(15)
	$\displaystyle\widetilde{\nabla_{\phi}Q}(\theta,\gamma,\phi|w^{(n)},y^{(n)})$	$\displaystyle=\sum_{i=1}^{n}\sum_{k=1}^{K}\biggl{(}\frac{\beta_{ik}}{\sum_{\ell=1}^{K}\beta_{i\ell}}\biggr{)}^{2}\left[\frac{\partial}{\partial x_{ik}}\log\left(\frac{p_{\theta,\gamma}(w_{i},y_{i},x_{ik})}{q_{\phi}(x_{ik}|w_{i},y_{i})}\right)\right]\frac{\partial}{\partial\phi}x(\phi,z_{ik}),$

where $\beta_{ik}=p_{\theta,\gamma}(w_{i},y_{i},x_{ik})/q_{\phi}(x_{ik}|w_{i},y_{i})$ and the expressions of $p_{\theta,\gamma}(w,y,x)$ and $q_{\phi}(x|w,y)$ are given by (4) and (6), respectively. The gradient computation involves calculation of $\nabla_{\theta}f_{\theta}$, $\nabla_{\gamma}g_{\gamma}$, and $(\nabla_{\phi}\mu_{\phi},\nabla_{\phi}\Sigma_{\phi})$, which are computed via the back propagation algorithm (McClelland et al.,, 1986; Hecht-Nielsen,, 1992). We also use the $L_{2}$-regularization trick to reduce numerical instability. We add a penalty $\lambda_{0}\|\theta\|^{2}+\lambda_{1}\|\phi\|^{2}+\lambda_{2}\|\gamma\|^{2}$ to the objective function. It is similar in spirit to the ridge regression and helps stabilize the numerical performance. This $L_{2}$ penalty changes nothing of the gradient ascent algorithm except for an additional term to each of the gradients $(\nabla_{\theta},\nabla_{\phi},\nabla_{\gamma})$; see Step 3 of Algorithm 1. The noise variance $\sigma^{2}$, which is assumed known in derivations of Sections 2.1-2.2, can be estimated along with training the neural network. At each epoch, we estimate $\sigma^{2}$ by a weighed sum of residuals:

$$\hat{\sigma}^{2}=\frac{1}{nK}\sum_{i=1}^{n}\sum_{k=1}^{K}(y_{i}-f_{\theta}(x_{ik}))^{2}\frac{\beta_{ik}}{\sum_{\ell=1}^{K}\beta_{i\ell}}.$$

(16)

No additional computing effort is needed, since $\beta_{ik}$ are already obtained in the computation of gradients.

Regarding the input of Algorithm 1, the measurement error density $p_{U}(\cdot)$ is supplied by the user, as in other measurement error methods. In the literature, it is common to assume that the measurement errors follow ${\cal N}(0,\sigma_{0}^{2}I_{d})$, where $\sigma_{0}^{2}$ is estimated from other data source or determined by prior knowledge. By default, we set the number of Monte Carlo samples as $K=50$ and the maximum number of epochs as Max_Epoch = 500. The neural network structure $T$ and the L2-penalty coefficients $(\lambda_{0},\lambda_{1},\lambda_{2})$ are chosen by a 5-fold cross validation. The validation loss is calculated as follows: For each test sample $(w_{i},y_{i})$, we draw $\{z_{ik}\}_{k=1}^{K}$ IID from ${\cal N}(0,I_{d})$ and compute $x_{ik}$ and $\beta_{ik}$ by plugging in the estimated parameters $(\hat{\theta},\hat{\gamma},\hat{\phi})$ and the testing $(w_{i},y_{i})$. These numbers are then plugged into (16) to give the validation loss. The implementation of gradient ascent is via Adam (Kingma and Ba,, 2014), with batch size equal to $\min\{512,n\}$. Adam requires three parameters $(\alpha_{0},\alpha_{1},\alpha_{2})$ to determine the learning rate (see the paragraph below (9)), where we set as the default values as $(0.001,0.9,0.999)$.

The method we propose here uses three neural networks, one for representing the regression function $f$, one for representing the (prior) distribution of $X$, and one for approximating the posterior distribution of $X$. If the prior distribution of $X$ is less complicated, we can model it with a parametric density $p_{\gamma}(x)$, where $\gamma$ now has a different meaning and represents parameters of the parametric density. Despite that the expression of $p_{\theta,\gamma}(w,y,x)$ changes, the objective (7) and the gradient ascent algorithm remain the same. It gives rise to a simplified version of NNME. This version works well if either the dimension is $1$ and/or the sample size is not very large. In such cases, using a parametric model for the distribution of $X$ leads to more stable numerical performance, even if the parametric model is misspecified (see Section 2.4). Throughout this paper, we still call this simplified version NNME but will clarify whether a neural network or parametric model is used for the distribution of $X$.

We compare NNME with alternative options of neural network algorithms for MEMs. In the first method (“NN”), we regress $y$ directly on $w$ using FNN, pretending that there is no measurement error. The second method is called the “Maximizing joint likelihood” method (“MJL”), which maximizes the joint log-likelihood of $(w,x,y)$. We use an FNN to approximate the MLE of $x$ given $(w,y)$ and another FNN to approximate the function $f(x)$. In comparison with NNME, the MJL approach is analogous to “optimal imputation,” which does not account for uncertainties in the imputation of $x_{i}$; see Appendix A for more details. In the third method (“VAE”), we apply the standard variational autoencoder, which corresponds to $K=1$ in the objective (7). To approximate the objective function at each iteration, the VAE requires multiple Monte Carlo samples, which are denoted as $z_{1},\ldots,z_{L}\sim{\cal N}(0,I_{d})$. Then, the objective in VAE is approximated by

$$\frac{1}{L}\sum_{\ell=1}^{L}\log\left(\frac{p_{\theta,\gamma}(w,y,x(\phi,z_{\ell}))}{q_{\phi}(x(\phi,z_{\ell})|w,y)}\right).$$

(17)

The implementation details are in Appendix A. The fourth method (“GA”, abbreviation for gradient ascent) maximizes (7) using the plain gradient estimator in (11). The last is our algorithm (“NNME”), where we maximize (7) with the doubly reparametrized gradient estimator. For all methods, we add L2 regularization on the weights in neural networks.

Below, we use two numerical examples to compare different neural network algorithms. We measure the performance by integrated squared error (ISE):

$$\int[f_{\hat{\theta}}(x)-f_{\theta}(x)]^{2}dx.$$

The three methods, VAE, GA and NNME, all require modeling of the prior distribtion on $X$. For a fair comparison, we use the same prior for all three methods, which is $2\cdot t_{3}$, where $t_{3}$ denotes a $d$-variate distribution whose each coordinate is an independent $t$-distribution with 3 degrees of freedom. We also use the same cross-validation procedure for selecting tuning parameters (including the neural network structure) in all the five methods.

Experiment 1: The univariate case. In this experiment, the true regression function is $f(x)=sin(\pi x)$. The data are generated from (1), where $x_{i}$’s are drawn from $unif(-2,2)$, the measurement error distribution is $N(0,\sigma_{0}^{2})$, and the noise standard deviations are $\sigma=\sigma_{0}=0.1$. The parameter $\sigma_{0}$ is supplied to all algorithms.

The results are shown in Figure 2 (a). From left to right, the sample size is $1000$, $2000$, and $5000$, respectively. Each of the three methods, VAE, GA and NNME, requires multiple Monte Carlo samples at each iteration. Note that

	number of MC samples
	$\displaystyle=\begin{cases}K\mbox{ (number of importance samples)},&\mbox{for GA and NNME},\\ L\mbox{ (number of repeated draws at each iteration)},&\mbox{for VAE},\end{cases}$

We have tried $K\in\{1,10,50,100\}$ and $L\in\{1,10,50,100,200\}$.

There are a few noteworthy observations. First, NN always has a higher error than the best of other methods. It suggests that ignoring measurement errors yields unsatisfactory performance. Additionally, MJL incurs a huge error when the sample size is $1000$ or $2000$, suggesting that the maximum joint likelihood approach only works for sufficiently large sample sizes (e.g., $n\geq 5000$). Second, we compare the three algorithms: VAE, GA and NNME, all aiming to maximize an ELBO of the marginal log-likelihood. The VAE objective is based on only one importance sample (i.e., $K=1$). It has a similar performance with NNME for $n=5000$, but its performance is worse than that of NNME for $n\in\{1000,2000\}$. Additionally, increasing $L$, the number of repeated draws at each iteration, yields no significant improvement. This suggests that IWAE (i.e., $K>1$) is indeed a better option. GA and NNME use different gradient ascent algorithms to solve the same IWAE objective. For $K=1$, the two have similar performances, but for larger $K$, GA is much worse than NNME. This shows that, without the doubly reparametrized technique in gradient estimation, the advantage gained by increasing the number of Monte Carlo samples is counterbalanced by the large variance in gradients. In contrast, increasing $K$ in NNME significantly improves the performance. For $n\in\{2000,5000\}$, NNME with $100$ importance samples achieves the best performance and is much better than the alternatives. For $n=5000$, MJL and VAE also perform well, and NNME performs similarly to them.

In Appendix A, we also investigate the estimated curve $f_{\hat{\theta}}(x)$ and the fitted values $\hat{x}$ for all methods, which help explain the reported ISEs. We relegate details to the appendix.

Experiment 2: The multivariate case. In this experiment, the true regression function is $f(x)=(x_{1}x_{2}+x_{3})^{2}$. We generate $\{x_{ij}\}_{1\leq i\leq N,1\leq j\leq 3}$ IID from $unif(-1,1)$, and let the measurement error distribution be ${\cal N}(0,\sigma_{0}^{2}I_{3})$, with $\sigma=\sigma_{0}=0.1$. Other settings are similar to those in Experiment 1.

The results are shown in Figure 2 (b). Compared with the univariate case, the performance of VAE is much worse because the ELBO of the marginal log-likelihood used by VAE becomes less tight for dimension higher than $1$. In (6), $\Sigma_{\phi}$ is set as a diagonal matrix, implying that a product measure is used to approximate the posterior distribution of $X$, which is restrictive for $d>1$. The two IWAE based methods, GA and NNME, help tighten the lower bound by increasing the number of importance samples, and perform better than VAE. In comparison, NNME is significantly better than GA due to the use of doubly reparametrized gradient estimator; this is the same as what we have observed in Experiment 1.

Experiment 3: Variants of NNME. In this experiment, we investigate a few variants of NNME that use different models for the prior (marginal) distribution of $X$. The true regression function is generated from a Gaussian process on $\mathbb{R}^{2}$ (see Appendix B for details). We simulate $\{x_{i}\}_{i=1}^{n}$ IID from a 2-component Gaussian mixture distribution:

$$x_{i}\overset{iid}{\sim}0.7\cdot{\cal N}\left(\begin{bmatrix}-0.4\\ 0.2\end{bmatrix},\;\begin{bmatrix}0.2^{2}&0\\ 0&0.3^{2}\end{bmatrix}\right)+0.3\cdot{\cal N}\left(\begin{bmatrix}0.2\\ 0.4\end{bmatrix},\;\begin{bmatrix}0.3^{2}&0\\ 0&0.2^{2}\end{bmatrix}\right).$$

(18)

We consider NNME with the following models on $X$: (a) True distribution of $X$ (benchmark). (b) A correct parametric model, i.e., assuming a 2-component Gaussian-mixture with unknown parameters. (c) A misspecified parametric model, e.g., a 4-component Gaussian-mixture or (after centering and scaling) a standard bivariate $t$-distribution. (d) The NICE model. The performance of every version is evaluated by the ISE, calculated on a uniform grid on $[-1,0.2]\times[-1,0.5]\cup[-0.5,1]\times[-0.2,1]$. The results are shown in Figure 3.

We consider settings where the sample size $n$ varies in $\{1000,4000,8000\}$ and the standard deviation $\sigma_{0}$ of measurement errors varies in $\{0.05,0.2\}$. Fixing a model for $X$, we also vary the structures of the two neural networks for representing $f$ (decoder) and approximating the posterior distribution of $X$ (encoder). Suppose the encoder has $\ell_{1}$ layers and the decoder has $\ell_{2}$ layers, with 32 nodes per layer; we consider four values of $(\ell_{1},\ell_{2})$, $\{(6,3)$, $(6,5)$, $(9,3)$, $(9,5)\}$, corresponding to the 4 box plots associated with each method in each panel of Figure 3.

When the measurement error is small (e.g., top panels of Figure 3), the performance of NNME is insensitive to the chosen model for $X$. In this case, the posterior distributions of $x_{i}$’s are primarily determined by the observed $w_{i}$’s, and the prior is uninformative. When the measurement error is large (e.g., bottom panels of Figure 3), if the sample size is large (e.g., $n\in\{4000,8000\}$), then using a correct parametric model (i.e., a 2-component Gaussian mixture) yields similar performance as using the true model. The NICE model yields slightly worse performance, but it does not require specifying any parametric model. Regarding the two misspecified parametric models, the 4-component Gaussian mixture performs as well as (sometimes even better than) the true model, but the (unimodal) $t$-distribution performs unsatisfactorily. If the sample size is small (e.g., $n=1000$), then the performances of different models are again similar. The results suggest a practical guideline for choosing the model for $X$: We can use either the NICE model or a parametric model that allows for enough modes.

Most classical methods for estimating MEMs start from a nonparametric regression method for the error-free case, and aim to address the bias caused by replacing $X$ with $W$ in the presence of measurement errors. Some notable approaches include the deconvolution approach (Fan and Truong,, 1993), the regression spline approach (Carroll et al.,, 1999; Jiang et al.,, 2018), and the simulation extrapolation approach (Carroll et al.,, 1999).

The deconvolution approach builds on local smoothing estimators for nonparametric regression in the error-free case. A local smoother takes the form $\hat{f}(x)=\sum_{i=1}^{n}\ell_{i}(x)y_{i}$, where $\ell_{i}(x)$ is determined by $x$ and $x_{1},x_{2},\ldots,x_{n}$. Given a kernel function $K(\cdot)$ and a bandwidth $h$, the Nadaraya-Watson estimator corresponds to $\ell_{i}(x)=K(\frac{x_{i}-x}{h})/[\sum_{j=1}^{n}K(\frac{x_{j}-x}{h})]$. In the presence of measurement errors, $\{K(\frac{x_{j}-x}{h})\}_{j=1}^{n}$ are unobserved. Fan and Truong, (1993) proposed replacing $K(\frac{x_{j}-x}{h})$ by $L(\frac{w_{j}-x}{h})$, such that $\mathbb{E}[L(\frac{w_{j}-x}{h})|x_{j}]=K(\frac{x_{j}-x}{h})$. The function $L(\cdot)$ is defined by

$$L(x)=\frac{1}{2\pi}\int e^{-\sqrt{-1}\,tx}\frac{\phi_{K}(t)}{\phi_{\epsilon}(t/h)}dt,\qquad\left\{\begin{array}[]{l}\phi_{K}(t):\mbox{Fourier transform of $K(\cdot)$},\\ \phi_{\epsilon}(t):\mbox{characteristic function of $\epsilon$}.\end{array}\right.$$

This method is inspired by a technique in density deconvolution problems (Carroll and Hall,, 1988; Stefanski and Carroll,, 1990) and, hence, called the deconvolution approach. This approach can be generalized to combine with other local smoothing estimators. For example, Delaigle et al., (2009) generalized the idea to the local polynomial estimator, where they proposed functions $\{L_{m}(\cdot)\}_{m\geq 0}$ such that $\mathbb{E}[(w_{j}-x)^{m}L_{m}(\frac{w_{j}-x}{h})|x_{j}]=(x_{j}-x)^{m}K(\frac{x_{j}-x}{h})$. In our numerical experiments, we include both the original deconvolution method (abbreviated as “deconv”) and the generalization to local polynomials (abbreviated as “lpoly”). See Table 2.

One advantage of the deconvolution approach is that it requires no structural assumption on $f(\cdot)$, except the smoothness. Another advantage is that it attains the minimax rate of convergence (Fan and Truong,, 1993). However, its practical performance crucially depends on the bandwidth selection. For density estimation with error-in-variables, Delaigle and Gijbels, 2004a ; Delaigle and Gijbels, 2004b proposed data-driven bandwidth selectors using cross-validation or bootstrap. There is relatively little work on bandwidth selection for nonparametric regression with error-in-variables. One exception is Delaigle and Hall, (2008), where they combined cross-validation with the simulation extrapolation idea (to be introduced below).

The regression spline approach approximates $f(\cdot)$ by spline bases. Different from the error-free case, fitting splines with measurement errors requires either estimating the posterior distribution of $X$ or calculating an unbiased score function. Carroll et al., (1999) used the truncated power basis and replaced $x_{i}^{m}$ and $(x_{i}-\xi_{j})_{+}^{m}$ by $\mathbb{E}(x_{i}^{m}|w_{1},\ldots,w_{n})$ and $\mathbb{E}[(x_{i}-\xi_{j})_{+}^{m}|w_{1},\ldots,w_{n}]$, where $\xi_{j}$’s are the fixed knots. To compute these conditional moments, they assumed that the marginal distribution of $X$ is a mixture of normals and developed a Gibbs sampling algorithm to obtain the posterior distribution of $X$. Jiang et al., (2018) used the B-spline basis and proposed an estimating equation method. Denote by $\beta$ the spline coefficients and let $S_{\beta}^{*}(W,Y,\beta)$ be the gradient of the log likelihood of $(W,Y)$, assuming a working density on $X$. This method finds a function $a(X,\beta)$ such that $\mathbb{E}\{S_{\beta}^{*}(W,Y,\beta)-\mathbb{E}^{*}[a(X,\beta)|W,Y]|X\}=0$, where $\mathbb{E}$ and $\mathbb{E}^{*}$ stand for the probability measure with $X$ following the true and working density, respectively. $\beta$ is estimated by the equation $\sum_{i=1}^{n}\{S^{*}_{\beta}(w_{i},y_{i},\beta)-\mathbb{E}^{*}[a(X,\beta)|w_{i},y_{i}]\}=0$. We abbreviate this method as “BSSP”; see Table 2.

A key assumption made by the spline approach is that the density of $X$ has a compact support (Hall and Qiu,, 2005). Under this assumption, the rate of convergence can be faster than the minimax rate: For Gaussian measurement errors, when $f(\cdot)$ has a continuous $k$th derivative, the minimax rate for the mean squared error is as slow as $\log(n)^{-k}$ (Fan and Truong,, 1993); but assuming a compact support for $X$, the spline approach can attain the standard nonparametric rate $n^{-\frac{2k}{2k+1}}$ (Jiang et al.,, 2018). The spline approach still requires choosing tuning parameters, which are the knots in spline construction. Following Eilers and Marx, (1996), we use a large number of equally spaced knots and add a penalty on the difference of coefficients of adjacent splines. It reduces to choosing the penalization parameter $\lambda$. Carroll et al., (1999) introduced a double-smoothing technique for selecting $\lambda$.

The simulation extrapolation (SIMEX) approach (Cook and Stefanski,, 1994) builds upon an arbitrary nonparametric regression method for the error-free case. It then replaces $w_{i}$ by $w_{i}(\lambda)=w_{i}+\sigma_{0}\sqrt{\lambda}\delta_{i}$, where $\delta_{i}$ is a ${\cal N}(0,1)$ variable independent of data and $\sigma^{2}_{0}$ is the variance of measurement errors. Let $\hat{f}(x;\lambda)$ be the estimated regression function by plugging in $\{w_{i}(\lambda)\}_{1\leq i\leq n}$ as covariates. In a parametric MEM, $\hat{f}(x;\lambda)$, as a curve of $\lambda$, gives a consistent estimator of the true $f(x)$ when extrapolated to $\lambda=-1$ (Cook and Stefanski,, 1994). This idea was generalized by Carroll et al., (1999) to nonparametric MEMs: It first computes $\hat{f}(x;\lambda)$ for a few values of $\lambda\geq 0$ using some nonparametric regression method, and then extrapolates to $\lambda=-1$. The extrapolation is often done by fitting a quadratic curve of $\lambda$.

SIMEX provides a convenient way to extend arbitrary nonparametric regression methods from the error-free case to the error-in-variable case. However, its consistency requires that the variance of measurement errors tends to zero. It was argued in Cook and Stefanski, (1994) that in a parametric MEM the simulation extrapolation with a quadratic extrapolant has an asymptotic bias of $O(\sigma_{0}^{6})$; hence, it is consistent only if $\sigma_{0}\to 0$ as $n\to\infty$. The use of SIMEX cannot avoid tuning parameter selection in the original nonparametric regression method. When the nonparametric regression method in SIMEX is the local polynomial estimator, Staudenmayer and Ruppert, (2004) introduced the empirical bias band width selector (EBBS) for bandwidth selection. EBBS estimates the bias of SIMEX for each fixed $h$ and selects $h$ that minimizes the mean-squared error.

All the above classical methods for fitting MEMs are primarily designed for dimension $d=1$. When $d>1$, although these methods have natural extensions, their implementations are much more challenging: the construction of kernels or splines, the computation of estimators, and the selection of bandwidth or knots, all become substantially more difficult. In our numerical experiments, we find it quite difficult to implement these methods when $d>1$; in particular, the computation of some methods is expensive even for a moderately large $n$. In contrast, our proposed neural network approach is a promising alternative, and can be cheaply implemented for a wide rage of $(d,n)$.

In numerical experiments we also include the kriging methods for Gaussian process regression as a competitor. Different from nonparametric regression, Gaussian process regression assumes that $f(x)$ is randomly generated from a stationary Gaussian process (SGP). It is common to use a radial basis exponential kernel in the SGP, i.e.,

$$Y=f(X)+\epsilon,\ \ f(x)\sim\mathrm{SGP}(0,K(\cdot,\cdot)),\ \ K(x,y)=\tau^{2}e^{-\beta\|x-y\|^{2}},\ \ \epsilon\sim{\cal N}(0,\sigma^{2}I_{n}).$$

(19)

Without measurement error, the best linear unbiased predictor of $f(x)$ at any given $x$ is $\hat{f}(x)=K(x,{\bf x}_{n})(K({\bf x}_{n},{\bf x}_{n})+\sigma^{2}I_{n})^{-1}{\bf y}_{n}$, where ${\bf y}_{n}=\{y_{i}\}_{1\leq i\leq n}$, $K(x,{\bf x}_{n})=\{K(x,x_{i})\}_{1\leq i\leq n}$, and $K({\bf x}_{n},{\bf x}_{n})$ is a $n\times n$ matrix whose $(i,j)$th entry is $K(x_{i},x_{j})$. The kriging method first estimates $(\sigma^{2},\tau^{2},\beta)$ by maximizing the likelihood and then plugs these parameters into the best linear unbiased predictor. In the presence of measurement errors, Cressie and Kornak, (2006) proposed a variant of the best linear unbiased predictor by replacing $K(x,{\bf x}_{n})$ by $\tilde{K}_{*}(x,{\bf w}_{n})$ and $K({\bf x}_{n},{\bf x}_{n})$ by $\tilde{K}({\bf x}_{n},{\bf x}_{n})$, where

$$\tilde{K}(w_{i},w_{j})=\begin{cases}\tau^{2}\frac{\exp\bigl{(}-\frac{\beta}{1+4\beta\sigma^{2}_{0}}\|w_{i}-w_{j}\|^{2}\bigr{)}}{(1+4\beta\sigma^{2}_{0})^{d/2}},&w_{i}\neq w_{j},\cr\tau^{2},&w_{i}=w_{j}\end{cases}\qquad\tilde{K}_{*}(x,w_{i})=\tau^{2}\frac{\exp(-\tfrac{\beta}{1+\beta\sigma^{2}_{0}}\|w_{i}-w_{j}\|^{2})}{(1+\beta\sigma^{2}_{0})^{d/2}}.$$

Same as before, $\sigma_{0}^{2}$ is the variance of measurement errors. The estimator of $f(x)$ is a plug-in version where $(\tau,\beta,\sigma)$ are estimated by maximizing a pseudo-likelihood. Although this method is not designed for nonparametric regression, we include it as a competitor of the neural network approach, especially for $d>1$. In our numerical experiments, we include both this kriging method that accounts for measurement errors (abbreviated as “KALE”), as well as the standard kriging for the error-free case (abbreviated as “KILE”).

We consider a 1-dimensional smooth regression function

$$f(x)=\frac{\sin((3x-1.5)\pi)}{1+4(6x-3)^{2}[\mathrm{sgn}(2x-1)+1]},$$

which was used in the simulation study of Berry et al., (2002). We generate data as in model (1), where $U\sim{\cal N}(0,\sigma_{0}^{2})$, $\epsilon\sim{\cal N}(0,\sigma^{2})$, and $X$ is generated from either a uniform distribution or a beta distribution.

We implemented all the methods in Table 2. Since the focus of this paper is the neural network approach, we prefer not to investigate tuning parameter selection for other methods in detail. To this end, we use the ideal tuning parameters for most classical methods, defined as the tuning parameters that minimize the true loss function (e.g., ISE), which biases slightly in favor of these classical methods. In Pspline and BSSP, we use cubic B-splines with equally spaced knots. We run the methods for different numbers of knots, from 5 to 30, and choose the one that minimizes the true loss function. Similarly, we run the methods for a range of values of the penalization parameter and select the ideal one. BSSP requires initial estimates of spline coefficients. We initialize by fitting B-splines on the observed data, ignoring measurement errors. In SIMEX, the original nonparametric regression method is the regression spline with truncated power bases. The selection of knots and penalization parameter, as well as initialization of coefficients, is similar to that for BSSP. In “deconv,” we set the bandwidth as $\sigma_{0}[\log(n)]^{-1/2}$, which is the theoretically optimal one for density deconvolution under Gaussian measurement errors (Stefanski and Carroll,, 1990). In “poly*” and “poly,” we set the polynomial degree as 1. We run the methods for 5 bandwidth values near $1.06\sigma_{0}n^{-1/5}$ and select the one that minimizes the true loss function.

The tuning parameter selection for NNME has been described in Section 2.3. In this simulation, we fix the values of tuning parameters, instead of using data-driven ones. We use 6 hidden layers in the decoder and 6 hidden layers in the encoder, where each layer has 32 nodes. We use a parametric model for the distribution of $X$, which is $2\cdot t_{3}$ (after centering and scaling), fix the L2-regularization parameters $\lambda_{0}$ and $\lambda_{1}$ as $10^{-5}$, and set $\lambda_{2}=0$.