{'fraction_non_alphanumeric': 0.04446692009904789, 'fraction_numerical': 0.015380322951905974, 'mean_word_length': 4.31532116044119, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Semiparametric Multivariate Density Deconvolution via Stochastic Rotation of Replicates

We consider the problem of multivariate density deconvolution where the distribution of a random vector needs to be estimated from replicates contaminated with conditionally heteroscedastic measurement errors. We propose a conceptually straightforward yet fundamentally novel and highly robust approach to multivariate density deconvolution by stochastically rotating the replicates toward the corresponding true latent values. We also address the additionally significantly challenging problem of accommodating conditionally heteroscedastic measurement errors in this newly introduced framework. We take a Bayesian route to estimation and inference, implemented via an efficient Markov chain Monte Carlo algorithm, appropriately accommodating uncertainty in all aspects of our analysis. Asymptotic convergence guarantees for the method are also established. We illustrate the method's empirical efficacy through simulation experiments and its practical utility in estimating the long-term joint average intakes of different dietary components from their measurement error contaminated 24-hour dietary recalls.


Introduction
Many practical applications require the estimation of the unknown density of a vectorvalued random variable x. The variable x, however, may not be observed precisely, rather surrogate replicates w contaminated with measurement errors u may only be available. The replicates w are then generated from a convolution of the density of x and the density of the measurement errors u, and the problem of estimating the density of x from available contaminated measurements w becomes a problem of multivariate deconvolution of densities. Such problems routinely arise in diverse application areas, including especially in nutritional epidemiology (Bennett et al., 2017;Keogh et al., 2020;Shaw et al., 2020), where the problem of estimating long-term average intakes of different dietary components from their errorcontaminated 24-hour recalls is of fundamental importance.
This article proposes a robust approach to multivariate density deconvolution in the presence of conditionally heteroscedastic errors u from an unknown probability law, relying on the idea of probabilistically rotating the replicates w toward the underlying latent x in a statistically principled manner.
Throughout this article, for random vectors s and t, we denote the marginal density of s, the joint density of (s, t), and the conditional density of s given t, by the generic notation f s , f s,t and f s|t , respectively. Likewise, for univariate random variables s and t, the corresponding densities are denoted by f s , f s,t and f s|t , respectively.
The literature on density deconvolution is really vast (Carroll et al., 2006;Buonaccorsi, 2010). The early literature focused primarily on univariate problems with a single contaminated measurement for each subject and the measurement errors independently and identically distributed according to some known probability law f u , often normal. Deconvoluting kernel-based approaches have been studied by Stefanski and Carroll (1990); Devroye (1989); Fan (1991a,b), among others. See also Madrid-Padilla et al. (2018); Newton (2002). The distribution of measurement errors is, however, rarely known in practice. Robust deconvolution methods with the unknown aspects of the error density estimated using replicated proxies w for the unknown values of x have thus been considered (Li and Vuong, 1998;Diggle and Hall, 1993;Delaigle et al., 2008, and others). A Bayesian likelihood-based approach with an unknown but symmetric unimodal density f x has recently been developed in Su et al. (2020).
The assumption of independence of u from x is also often highly impractical, especially in nutritional epidemiology applications, where patterns of conditional heteroscedasticity can be very prominently seen. Bayesian hierarchical frameworks and associated Markov chain Monte Carlo (MCMC) based computational machinery have recently been shown to provide powerful tools for solving complex deconvolution problems under more realistic scenarios, including when the measurement error distribution can be asymmetric, heavy-tailed, conditionally heteroscedastic, etc. (Staudenmayer et al., 2008;Sarkar et al., 2014Sarkar et al., , 2018Sarkar et al., , 2021. In their seminal work, Staudenmayer et al. (2008) assumed the errors u to be normally distributed but allowed the variability of u to depend on x, employing positive mixtures of B-splines to flexibly characterize both f x and the conditional variability var(u|x). Sarkar et al. (2014) further relaxed the assumption of normality of u, employing flexible mixtures of normals (Escobar and West, 1995;Frühwirth-Schnatter, 2006) to model both f x and f u|x .
In stark contrast to the univariate setting, the multivariate problem has garnered little attention in the literature. Masry (1991); Youndjé and Wells (2008); Comte and Lacour (2013); Turlach (2009, 2010); Eckle et al. (2017) considered scenarios with errors u from a known probability law, independent from x. Bovy et al. (2011) obtained a Bayesian maximum-a-posteriori estimate of f x modeled by flexible mixtures of multivariate normal kernels, assuming f u to be multivariate normal, independent from x with subjectspecific but known covariance matrices. Utilizing the flexibility of Bayesian hierarchical frameworks, Sarkar et al. (2018Sarkar et al. ( , 2021 developed robust multivariate deconvolution methods, relaxing the restrictive assumptions of known error probability laws, homoscedasticity, and independence from x, etc. Sarkar et al. (2018) modeled f x and f u|x using flexible mixtures of multivariate normals whereas Sarkar et al. (2021) adopted a complementary approach, modeling the marginals f x and f u |x first and then building the joint distributions f x and f u|x by modeling the dependence structures using Gaussian copulas.
The focus of this article is also on multivariate deconvolution with conditionally heteroscedastic measurement errors from an unknown distribution in the presence of replicated proxies for each subject. To that end, we propose a novel approach to multivariate density deconvolution that assumes the replicates w to be generated by first stochastically rotating the underlying true x and then stochastically stretching or contracting their lengths. This is achieved by multiplying each x first with an orthogonal rotation matrix Q and then with a scalar length adjustment factor r. Going a significant step further, we also accommodate conditional heteroscedasticity by allowing the distributions of both the rotation matrices and the length-adjusting factors to flexibly depend on the latent true x's. The conditional distributions f u|x , and hence f w|x , are then obtained as novel functions of the Q's. For the main density of interest f x , we adopt a copula-based approach with the marginals modeled by flexible mixtures of truncated normals with shared atoms as in Sarkar et al. (2021). We take a Bayesian route to estimation and inference, implemented via an efficient MCMC algorithm, appropriately accommodating uncertainty in all aspects of our analysis. We illustrate our method's empirical efficacy through simulation experiments. Its practical utility is demonstrated in nutritional epidemiology applications in estimating the long-term joint average intakes of different dietary components from their measurement error-contaminated 24-hour dietary recalls.
Traditionally, the literature on deconvolution almost exclusively assumes the measurement errors to be additive. In Section 3.2 of this paper, we show that our rotation-based model can be reformulated as a classical additive model. What the rotation-based view does still is to provide a new perspective on measurement errors leading to a new way of constructing the likelihood function and resulting in new algorithms for inference.
Rotation-guided modeling of multivariate data is indeed getting increasing popularity in statistics (Hoff, 2009a,b;McCormick and Zheng, 2015;Mukhopadhyay et al., 2020;Song and Dunson, 2022). In the deconvolution literature, stochastic rotations have been proposed for directional data in Kim (1998), where the latent objects of primary interest, as well as the observed replicates, were all orthogonal matrices, with additional theoretical insights into this setup provided in Kim and Richards (2001). Our work, however, is focused on Euclidean deconvolution problems in the presence of replicates contaminated with conditionally heteroscedastic measurement errors. The deconvolution problem we consider and the solution we propose are thus very different from Kim (1998).
Overall, this article makes several important contributions to the literature on multivariate density deconvolution -(a) we introduce a new framework for multivariate deconvolution via stochastic rotation of the error-contaminated replicates toward their underlying true values, (b) additionally, we also address the significantly challenging problem of accommodating conditionally heteroscedastic errors in this newly introduced framework, and (c) we introduce HMC-based advanced MCMC methods to the deconvolution problem, significantly improving computational efficiency.
The rest of this article is organized as follows. Section 2 presents some important preliminary results used in the construction of our likelihood function. Section 3 details our proposed stochastic rotation-based approach to multivariate deconvolution, including likelihood construction, prior specification, and outline of posterior computation. Section 4 presents the results of some simulation experiments, illustrating the method's empirical performances. Section 5 presents the results produced by the proposed method applied to the problem of estimating the true long-term average intakes of different dietary components from their measurement error-contaminated 24-hour recalls. Section 6 contains concluding remarks. Substantive additional details are presented in the supplementary materials.

Geometry of Vector Rotations
To motivate our modeling framework, we first discuss some geometric properties of vector rotations. For any two vectors a ∈ R d and b ∈ R d , there exists a d × d orthonormal rotation matrix Q ab such that a a 2 = Q ab b b 2 , where · 2 stands for the Euclidean norm. We can thus establish that for the scalar s = a 2 / b 2 . Geometrically speaking, the orthonormal matrix Q ab rotates the unit vector b b 2 towards the unit vector a a 2 and the scalar s takes care of the change in magnitude due to this transformation from b to a. Let a = a/ a 2 and b = b/ b 2 . Then, using the Householder reflection result (Householder, 1958) for total internal reflection under Snell's law, a solution for Q ab is However, since it is based on a single Householder transformation, the solution is always symmetric, hence making it difficult to impose any distributional assumption. We thus consider a different solution based on a frequently used technique in numerical analysis to compute the orthogonal component Q in the QR-decomposition of a matrix (Stewart, 1980). The exact expression of Q ab is established in the following theorem.
Theorem 1. For any i = 1, . . . , d, let a i = a + a 2 e i and b i = b + b 2 e i , where e i is the unit vector with 1 at the i th place. Then, for (1). The result above is crucial in computing the likelihood function of our proposed model. The proof is based on some results from Euclidean geometry. If the singular value decomposition of ab T is UΣV T with the singular values in Σ = diag(σ 1,1 , . . . , σ d,d ), then a solution for Q ab is the orthogonal matrix UV T . The matrix ab T has only one non-zero singular value, a 2 b 2 , which we can assume to be σ i,i without any loss of generality. Based on a similar reflection argument as before, we can then compute U and V as functions of a and b explicitly. Specifically, we take a i = a + a 2 e i , the bisector of the angle between a and the i th unit vector e i on the reflecting surface. Then the Householder transformation matrix based on this bisector, , is a possible solution for U as this Householder reflection will transform a to − a 2 e i and thus Ua = − a 2 e i . It can be verified easily by noting that For the rest of the paper, without loss of generality, we use the result of Theorem 1 with the first unit vector e i = e 1 . Numerical experiments with other choices produced nearidentical results. In the next section, we use the result to compute the likelihood function of the replicates in our multivariate density deconvolution model introduced in the next section.

von-Mises Fisher and Matrix von-Mises Fisher Dstributions
The von-Mises Fisher distribution (vMF) and the matrix von-Mises Fisher distribution (MvMF) play important roles in our model construction and computation. We thus provide a brief description of these distributions here, starting with the vMF first, for easy reference.
The vMF distribution for a d-dimensional vector w is defined as where µ is a unit vector signifying the mean direction of w and c ≥ 0 stands for a scalar concentration parameter. The normalizing constant C d (c) = the modified Bessel function of the first kind at order ν. The maximum likelihood estimates for µ and c are obtained in Jupp and Mardia (1979). We next review the MvMF distribution for orthogonal matrices on the Stiefel manifold. The distribution was first proposed in Downs (1972) and thoroughly studied in Khatri and Mardia (1977). Since then, it has become a popular distributional choice for random matrices supported in the space of orthogonal matrices.
The MvMF distribution for a q × d matrix Q = ((Q 1 , 2 )), with q ≤ d and QQ T = I q , is defined as where etr(·) = exp{trace(·)}, F is q × d dimensional parameter matrix, O q,d is the set of all q × d orthogonal matrices (i.e., the q × d Stiefel manifold), and M (F) is the normalization constant. To keep the description simple, we only consider the case directly relevant to our deconvolution problem, namely q = d with a diagonal F, and review some key results from Khatri and Mardia (1977). The MGF of Q is given by E{etr(TQ T )} = M (F+T) M (F) . Let T 1 = diag(t 1,1 , . . . , t 1,d ) and T 2 = diag(t 2,1 , . . . , t 2,d ) be matrices such that t 1,ell = 1/t 2, for all = 1, . . . , d. Thus, the distributions of Z = T 1 QT 2 and Q are identical when F is diagonal. Hence, This implies E(Q 1 , 2 ) = 0 for all k = and E(Q 1 , 2 Q 1 , 2 ) = 0 for either 1 = 2 or 1 = 2 . To compute the expectation of the diagonal entries, we can take the derivative of the MGF with respect to T. Hence, we get E(Q 1 , 1 ) =

Deconvolution via Stochastic Rotation
Our main objective is to estimate the density of a d-dimensional vector x. However, we do not have accurate measurements of x. For each unobserved x i , we instead have m i replicated proxies w i,j 's contaminated with some error where i = 1, . . . , n and j = 1 . . . , m i , and m i ≥ 3 for all i. Each observation w i,j may be viewed as a point in the Cartesian co-ordinate system. A representation of w i,j can then be obtained in terms of its norm w i,j 2 and its direction from the origin w i,j / w i,j 2 . In this paper, our characterization of w i,j takes inspiration from the above representation and the relation in (1). Specifically, our model for the replicates w i,j conditional on the underlying true x i is Here r i,j 's are scalars with E(r i,j |x i ) = 1 and var(r i,j |x i ) = exp{s 2 ( x i 2 /d)} − 1, Q i,j = ((Q , ,i,j ))'s are d × d rotation matrices, and B k (x) are B-spline bases spanning the interval [A, B] (de Boor, 1978). The d × d dimensional scaling matrices C i = ((C , ,i )) ensures that the replicates w i,j 's are stochastically centered around the corresponding true x i . Specifically, As seen in Section 2, for a diagonal F i , the Euclidean expectation of Q i,j is also diagonal. Based on the results of Khatri and Mardia (1977), we specifically have E(Q , ,i,j |x i ) = The matrices F i also determine how far the corresponding Q i,j 's are allowed to vary around I d , larger values of κ (·) inducing greater concentration of the Q i,j 's around I d . Throughout the paper, we often keep the Rewriting the model as C −1 i w i,j = r i,j Q i.j x i and following Theorem 1 in Section 2, the solutions for Q i,j and r i,j are G C −1 i w i,j , x i and C −1 i w i,j 2 / x i 2 , respectively. The conditional likelihood of the replicates is therefore given by Although M (F i ) is computationally intractable, depending on the magnitudes of the In the alternative view adopted in this article, w is generated by first rotating x to Qx and then stretching (left panel) or contracting (right panel) it to w. The 3-dimensional case is shown here, the shaded spheres representing the set of all vectors that can be generated by orthogonal rotations of x.
elements in F i , efficient approximations are available in Khatri and Mardia (1977). In most practical applications, including ours, the replicates w i,j are expected to lie in the general direction of the underlying x i , which implies that the Q i,j 's shouldn't deviate too far from the identity matrix I d . Throughout this article, the functions κ i (·)'s are thus assumed to take large values which result in a higher concentration of Q i,j around I d . Thus, it is reasonable to consider the following approximation from Khatri and Mardia (1977) for large diagonal In what follows, for brevity, sometimes we also use the notation Also, when we say F i is large, we mean its diagonal entries f i are all large. Following Sarkar et al. (2021), we model the joint density f x of x using a Gaussian copula with component-wise univariate marginals characterized by flexible mixtures of truncated normals, truncated to their common support [A, B]. Specifically, we let where y x = (y x,1 , . . . , y x,d ) T with F x, (x ) = Φ(y x, ) for each = 1, . . . , d, where F x, is the cumulative distribution function (cdf) corresponding to f x, . Our model for the correlation matrix R x considers the spherical coordinate representation of Cholesky factorization. In this representation, the correlation matrix R x is written as R . . , m. The rest of the entries are all zeros. Thus, V is a lower triangular matrix. The angles ζ m,s are supported on [0, π] for s < m − 1 and the ζ m,m−1 's are supported on [0, 2π].
The truncated support [A, B] for the marginals is consistent with modeling conditional heteroscedasticity later in Section 3.1 by mixtures of B-splines which by definition have bounded local supports spanning a finite interval. Having a common unit free support [A, B] for all components also greatly simplifies assigning priors on the parameters of the mixture components as well as the choices of these prior hyper-parameters. This is often easy to achieve in practice via simple linear transformations of the originally observed proxies. See, e.g., Section S.5. in the supplementary material.

Conditional Heteroscedasticity Characterization
There are two different ways we accommodate conditional heteroscedasticity in the measurement errors in our model -(a) by allowing the distribution of length adjustment factors r i,j to depend on the underlying latent x i via the function s 2 ( x i 2 /d), and (b) by allowing the κ parameters to vary flexibly as functions of the corresponding latent x ,i . The part s 2 ( x 2 /d) accommodates our expectation that larger adjustments are needed for larger true latent vectors x; whereas the functions κ (x ,i )'s accommodate the expectation that larger deviations of w ,i,j should be allowed around larger values of x ,i . More specifically, we have cov(w i, is the covariance matrix for the diagonal entries in Q i,j . From equation (2.11) of Khatri and Mardia (1977), we have For a large F i , relying on the approximation of M (F i ), we can compute the entries V , ,i as To simplify notation, let us denote cov(w i,j |x i ) also by S i = ((S , ,i )). The other approximations are then Due to the latter component, a large f ,i may not always ensure low conditional variability in w ,i,j . If, however, smaller x ,i 's correspond to larger values of f ,i , then we have smaller conditional variability in w ,i,j . As discussed in detail in Sarkar et al. (2018), by the very nature of such problems, the conditional variability of each component w ,i,j should depend primarily on the corresponding latent component x ,i . Interestingly, however, unlike Sarkar et al. (2018), in our model the expected componentspecific conditional variances S , ,i of w ,i,j involve contributions from all entries of x i . This is not surprising as our characterization of the joint distribution of w rely on two conditionally independent distributions, namely the distributions of its norm and its direction given the unobserved x. And the norm and direction of a vector are functions of all of its entries. Primary dependence of S , ,i on x ,i is, however, still accommodated via the functions κ (x ,i ) and the resulting diagonally dominant nature of V i . Staudenmayer et al. (2008); Sarkar et al. (2014Sarkar et al. ( , 2018Sarkar et al. ( , 2021 also showed that the variability in w ,i,j usually increases with x ,i , especially in dietary recall data. It would have thus been ideal to have the functions κ (x ,i )'s be non-increasing in x ,i and the function s 2 ( x 2 /d) increasing in x . Although we have not imposed such shape constraints explicitly in our model, the estimated functions do exhibit such behaviour in all our simulated and real data applications. Even greater flexibility in the conditional covariance structure may be obtained by taking mixtures of MvMF distributions instead of one single MvMF. Such extensions will, however, be pursued elsewhere.

Connections with Classical Models
Our proposed model can be rewritten as a classical additive measurement error model as The additive measurement error u i,j is thus the vector joining the tips of x i and w i,j ( Figure 2). However, unlike previous works on conditionally varying measurement errors such as Sarkar et al. (2018Sarkar et al. ( , 2021, the entries of w i,j and hence those of u i,j are allowed to depend on all components of x i . The proposed model also has a resemblance with classical multiplicative measurement error models w i,j = x i • u i,j where • denotes element-wise product and the errors u i,j are distributed independently of x i with E( u i,j ) = 1. In our model, however, we have In both specifications, the scalar part r i,j is shared by all components of the measurement errors u i,j and u i,j . However, further interactions between the entries of x i and the entries of u i,j or u i,j may be observed depending on the concentration of Q i,j around I d . For example, when Q i,j = I d for all j, there is only an element-wise effect of x i 's on the errors. However, as Q i,j moves away from I d , there is more inter-component associations between the x i 's and the associated errors u i,j or u i,j .
In the following subsection, we show that a classical additive normal measurement error model can be obtained as a meaningful limiting case of our rotation-based model.

Limiting Equivalence with a Normal Error Model
To reduce notation complexity, we initially fix r i,j = 1 for all i, j in this section. The We recall that F i = diag(f i ) and we have argued that in practice its entries are expected to be large. We now study the limiting distribution of each coordinate of u Without any loss of generality, we set t = (t 1 , 0, . . . , 0). We have E{exp i C i will be diagonally dominant and its off-diagonal entries will be very small. Hence, the eigenvalues of F i + x i t T H 1/2 i C i will be its diagonal entries. Applying the approximation from (5), we then have applying Taylor series expansion. Since C 1,1,i ≈ 1 1 − k =1 1 2(f 1,i +f k,i ) , we can simplify the above expression further to exp(h Hence, each coordinate of u i,j marginally approaches to Normal(0, 1). Thus, in a limiting sense, the th coordinate of w i,j r i,j marginally reduces to a normally distributed classical measurement error model with mean x ,i and variance h −1 ,i .

Error Distribution Generating Function
To characterize the distribution of u i,j , we can compute the moment generating function (MGF) of (u i,j |x i ) using the MGF result in equation (2.7) of Khatri and Mardia (1977) as where M (·) is the normalizing constant whose approximation is given in (5) for large F i . To study the tail properties of each component in u i,j , we can consider different choices of t. Specifically, for studying the th component, the choice is t = (0, . . . , t , . . . , 0) T . Assuming F i is large, the approximation in (5) gives us Explicit characterization of the MGF beyond the above expression is difficult to obtain. We see, however, that the MGF might be undefined for positive valued t for the heavy-tailed lognormal distribution on the r i,j 's. Other light-tailed choices, such as a gamma distribution with the rate parameter modeled as a function of x i 2 , may also be considered. In this article, however, we focus on the log-normal.

Bayesian Inference
While the main idea underlying our rotation-based approach to multivariate deconvolution is clearly statistical paradigm generic, in this article, we adopt a Bayesian route to estimation and inference. The rest of this section discusses prior specification, posterior computation and posterior convergence in such settings.

Prior Specification
Since the functions s 2 (·) and κ (·)'s are strictly positive, we put truncated normal priors, truncated to [0, ∞), on the associated B-spline coefficients: . We assign Gamma priors on the inverse variances σ −2 . For the component specific parameters (µ k , σ 2 k ) of the truncated normal mixtures characterizing the marginals f x , we assign semi-conjugate independent priors µ k ∼ Normal(µ 0 , σ 2 0 ) and σ −2 k ∼ Ga(a 0 , b 0 ). We assign independent Dirichlet priors on the mixture probabilities π = (π ,1 , . . . , π ,K ) T ∼ Dir(α/K, . . . , α/K), and independent uniform priors on the polar angles ζ m,j ∼ Unif[0, π] and ζ m,m−1 ∼ Unif[0, 2π]. The choice for the hyper-parameters are a s = b s = a κ, = b κ, = 0.1. The hyper-parameter choices for µ s,k and µ κ, ,k are obtained based on the estimated x i 's from the univariate sampler. The hyper-parameters µ 0 and σ 0 from the copula model of x are also set based on these univariate estimates. The other two hyper-parameters a 0 and b 0 are set as a 0 = b 0 = 1. The univariate sampler follows the additive model of Sarkar et al. (2021). Details are provided in Section 3.4.2.

Posterior Computation
Our inference is based on samples drawn from the posterior using an MCMC algorithm. The joint log posterior distribution of the model parameters is given by It is possible to efficiently calculate the derivatives of the above likelihood with respect to β s , µ ,k , σ ,k , etc. These parameters may therefore be updated using HMC algorithms (Neal, 2011;Betancourt and Girolami, 2015;Betancourt, 2017). HMC has been shown to draw posterior samples much more efficiently than traditional random walk Metropolis-Hastings in complex Bayesian hierarchical models (Betancourt and Girolami, 2015) by more efficiently exploring the target distribution under local correlations among the parameters.
A short review of HMC sampling is provided in Section S.2 in the supplementary materials for easy reference. The x i 's, and the parameters specifying the R x of f x (x i ), can be updated using random walk Metropolis-Hastings steps. The parameters β κ are updated using adaptive Metropolis-Hastings. Other parameters have closed-form full conditionals and can be straightforwardly updated. Details are deferred to Section S.5 in the supplementary materials.

Simulation Experiments
In this section, we evaluate the performance of our proposed method, referred to in the tabular result summaries below as the DeStoR method, short for 'Deconvolution via Stochastic Rotation'. We compare with the method of Sarkar et al. (2021), referred to as the SPMC method here following the first letters of the authors' last names, which was shown to vastly outperform the only two previously existing multivariate deconvolution methods accommodating heteroscedastic measurement errors, namely the pseudo-Bayesian two-stage method of Zhang et al. (2011) and the multivariate mixture normal based method of Sarkar et al. (2018). We consider two different simulation scenarios -(a) a well-specified case; and (b) a mis-specified case. For both scenarios, we consider the same f x to generate the true latent x i 's but use different choices for f w|x to generate the replicates w i,j 's. Our choice for the well-specified case conforms to our the proposed formulation for f w|x in (4) in Section 3. The mis-specified case is designed to evaluate the robustness of the proposed method to deviations from model assumptions and considers additive measurement errors generated from a Gaussian copula model with a mean restricted mixture of normals for the marginals from Sarkar et al. (2021) to produce the replicates w i,j 's.
While our proposed method scales well to much higher dimensional problems, we consider a relatively low d = 3 dimensional problem here as the computation of the joint density on a 3 dimensional grid remains manageable and the results for 3 2 = 3 bivariate marginals can also be conveniently graphically summarized. We generate the true x ,i 's for = 1, . . . , d as follows. We (a) first sample . This way, the marginal distributions are mixtures of truncated normal distributions and hence can take widely varying shapes while the correlation between different components is R x . We set x, ,k = 0.75 2 for all , k. To generate w i,j for the well-specified case based on (4), we need to generate Q i,j and r i,j given x i that are generated in the previous step. We use the R package rstiefel (Hoff and Franks, 2019) to generate Q i,j with concentration parameters κ (x ,i )'s where κ (x ,i ) = 60/x ,i . The scalars r i,j 's are generated from log-normal distribution with mean = s 2 ( x i 2 /d)/2 and variance = s 2 ( While generating w i,j for the mis-specified case based on Sarkar et al. (2021), we consider the additive model w i,j = x i + s 1 (x i ) • i,j , where s 1 (x i ) = x i /4 and the i,j 's are generated from a standard multivariate normal distribution.
The integrated squared error (ISE) of estimation of f x by f x is defined as are random samples from the density p 0 . We used the true densities f x for p 0 and the true values of the x i 's for the x m 's. For the univariate marginals, a Monte Carlo estimate of ISE is given by are a set of grid points on the range of x and ∆  Table 1 reports the median ISEs (MISEs) for estimating the trivariate joint densities and the univariate marginals obtained by our method and the method of Sarkar et al. (2021). The reported MISEs are all based on 100 simulated data sets. In the well-specified case, when the data-generating mechanism conforms to our proposed model, our method significantly outperformed Sarkar et al. (2021) in estimating the three-dimensional joint density and as well as all univariate marginals. In the mis-specified case, when the data generating mechanism conforms to the model of Sarkar et al. (2021), our method still performed competitively with Sarkar et al. (2021) for the three-dimensional joint density estimation problem and actually outperformed Sarkar et al. (2021) for some of the univariate marginal density estimation problems. We attribute this to our more efficient MCMC sampling of the posterior via HMC samplers. Figures 3 and 4 show the estimates of the univariate marginal and bivariate joint densities obtained by our method for the data set that produced the 25 percentile ISE and the Figure 3: Results for the simulated data with n = 1000 subjects and m i = 3 replicates per subject when true data generating process follows the structure (4) proposed in this article. The figures in the diagonal panels illustrate the true marginal densities f x, (x ) in red and the corresponding estimates produced by our method in blue. The figures in the off-diagonal panels depict the contour plots of the true two-dimensional marginals (upper triangular panels) and the corresponding estimates obtained by our method (lower triangular panels). Figure 4: Results for the simulated data with n = 1000 subjects and m i = 3 replicates per subject when true data generating process follows the additive model from Sarkar et al. (2021). The figures in the diagonal panels illustrate the true marginal densities f x, (x ) in red and the corresponding estimates produced by our method in blue. The figures in the off-diagonal panels depict the contour plots of the true two-dimensional marginals (upper triangular panels) and the corresponding estimates obtained by our method (lower triangular panels).
corresponding truths for the well-specified and the mis-specified cases, respectively. The estimates clearly provide excellent fits to the truths. Additional plots comparing the estimates of the univariate marginals produced by our method with those obtained by the method of Sarkar et al. (2021) are presented in Figures S.3 and S.4 in Section S.9 of the supplementary materials.
We also evaluate the numerical performance of the proposed method in some higher dimensional cases, namely d = 5 and d = 10. We maintain a similar structure for R x and set σ 2 x, ,k = 0.75 2 for all the univariate components as before. The mixture probabilities are also kept the same, as the univariate distributions are again assumed to be mixtures of three univariate normals. However, since it is difficult to set the mean parameter µ x in higher dimension explicitly, we set the components of µ x using random numbers generated from Unif(1, 10). We summarize the results in Table 2. For clarity of our presentation, here we only present the MISEs for the d-dimensional joint densities. The results show patterns similar to that previously observed for the d = 3 case.

Real Data Application
The estimation of the joint and marginal distributions of long-term average daily intakes of different dietary components is a fundamentally important problem in nutritional epidemiology. The long-term average daily intakes of the dietary components, x, can not, however, be directly measured. Data are thus often collected via nutritional surveys in the form of dietary recalls, the subjects participating in the study remembering and reporting the type and amount of food they consumed in the past 24 hours. The problem of estimating the joint consumption pattern of the dietary components from the contaminated 24-hour recalls then becomes a problem of multivariate density deconvolution. One such large-scale nutritional survey is the Eating at America's Table (EATS) study (Subar et al., 2001) conducted by the National Cancer Institute where n = 965 participants were interviewed m i = 4 times over the course of a year and their 24 hour dietary recalls w i,j were recorded. The goal is to estimate the joint and marginal consumption patterns of the underlying true daily intakes x i . In this article, we focus particularly on the average daily intakes of three minerals, namely, iron, magnesium, and sodium. Figure 5 presents the estimated univariate and bivariate marginals obtained by our proposed method and the method of Sarkar et al. (2021). The estimates of the univariate marginals produced by the two methods are quite similar. The estimated bivariate densities, however, although share some commonalities across the two methods, are not exactly the same. Specifically, the estimates produced by the method of Sarkar et al. (2021) are more strongly correlated than ours.
We conducted a formal model comparison exercise to identify which model fits the EATS data set better. Since the competing methods are both Bayesian, we compared them in terms of the Bayes factor (Kass and Raftery, 1995). The Bayes factor comparing the two methods is given by B = P (D (i) |DeStoR) P (D (i) |SPMC) . We calculated the numerator and denominator from the posterior samples using the harmonic mean identity of Neton and Raftery (1994). Based on the suggestions in Kass and Raftery (1995), the evidence in favor of our proposed DeStoR is 'decisive' as we obtained B > 100.

Discussion
In this article, we developed a novel method for multivariate density deconvolution in the presence of conditionally heteroscedastic measurement errors and the availability of replicated proxies for the unknown values of the variable of interest. Our proposed method relies on stochastically rotating the observed replicates toward the underlying true values and then stochastically adjusting their lengths to match the lengths of the true values. We took a Bayesian route to estimation and inference, implemented via an efficient MCMC algorithm. In synthetic numerical experiments, the proposed method showed excellent performance in recovering the true density of interest. The method's practical utility was demonstrated in a nutritional epidemiology application in estimating the joint distribution of the true average long-term intakes of three different dietary components from their measurement error-contaminated 24-hour recalls.
In Section S.6 of the supplementary materials, we show that the posterior of our model convergences to the true unknown density of the observed data w i as the sample size grows to infinity. There is substantial literature studying the convergence properties of the posteriors in ordinary density estimation problems where f w is directly modeled. Our results are, however, established under a conditionally heteroscedastic measurement error setting as described in (4), where f w is obtained by a convolution of the models for f x and f u|x . To our knowledge, our results are novel to the literature, especially under dependence between the measurement error u and the unknown true vector of interest x. In the context of density deconvolution, however, our results lack strong inferential merit as the notion of recovery is in terms of f w but not f x . In the future, we hope to study posterior consistency with a more appropriate notion of recovery, that of the density of interest f x , under density deconvolution with conditional heteroscedastic matrix vMF distributed measurement errors.
While this is our first work in this line of research, our proposed approach is already quite general, not only allowing completely unknown measurement error distributions but also accommodating unknown conditional heteroscedasticity patterns. We thus believe the work contributes a novel and significant addition to the existing sparse literature on multivariate density deconvolution.
Directions for methodological extensions and topics of our ongoing research include exploration of other distributions on the Stiefel manifold, including mixtures of MvMF, adaptations to regression problems with errors-in-covariates, etc.

Supplementary Material
The supplementary material presents brief reviews of copulas, and Hamiltonian Monte Carlo, and also the explicit formula of cubic B-splines for easy reference; discusses model identifiability; and details the choice of hyper-parameters, the MCMC algorithm used to sample from the posterior, and some convergence results for the posterior and their proofs. The supplementary material also presents some additional figures summarizing the results of the simulation experiments. R programs implementing the deconvolution methods developed in this article are included in the supplementary material. The EATS data analyzed in Section 5 can be accessed from National Cancer Institute by arranging a Material Transfer Agreement. A simulated data set, generated using the estimates produced by our method for the EATS data set, and a 'readme' file providing additional details are also included in the supplementary material. guidance document on measurement error and misclassification of variables in observational epidemiology: Part 2 -more complex methods of adjustment and advanced topics. Statistics in Medicine, 39, 2232-2263. Song, H. and Dunson, D. B. (2022. Curved factor analysis with the ellipsoid-gaussian distribution. arXiv preprint arXiv:2201.08502 .  Supplementary material presents brief reviews of copula basics, Hamiltonian Monte Carlo, and cubic B-splines to make the article relatively self-contained. Supplementary material also discusses model identifiability; and presents details of the MCMC algorithm we designed to sample from the posterior; some convergence results for the posterior and their proofs; and some additional figures. Separate files additionally include a synthetic data set, generated using the estimates produced by our method for the EATS data analyzed in Section 5 of the main paper, and R programs implementing the proposed multivariate density deconvolution model developed in this article.

S.1 Gaussian Copula
The literature on copula models is enormous. See, for example, ??? and the references therein. For easy reference, we provide a brief review of the basics here.

S.2 Hamiltonian Monte Carlo
In this section, we present a brief review of the Hamiltonian Monte Carlo (HMC) sampling algorithm for easy reference. Our presentation of HMC is inspired from Neal (2011) and Betancourt and Girolami (2015). Let us assume that we want to draw samples of q from π(q). The HMC algorithm adds an additional momentum variable p and considers following joint density π(p, q) = π(p|q)π(q).
The Hamiltonian H(p, q) of a system is defined as the sum total of the kinetic energies and the potential energies of all the particles in the system. Liouville's theorem on phase space equation from statistical mechanics states that the probability density of phase space is time-invariant along the trajectory. Based on that, it can be shown that for multivariate normal initial distribution the joint density π(p, q) ∝ exp{−H(p, q)}.

Thus we get
where K(p|q) is the kinetic energy and U (q) is the potential energy. The first step of the HMC algorithm is to draw the momentum variable p such that p ∼ π(p|q). The changes in p and q over time are governed by the following Hamiltonian's equations dq dt = + ∂H ∂p = + ∂K ∂p , dp dt = − ∂H ∂q = − ∂K ∂q − ∂U ∂q .
Computational implementation of these equations requires discretization of time, with small step size and evaluate the states at t = , 2 , . . . , so on. Neal (2011) showed the superiority of the leapfrog method in producing excellent solution to a system of differential equations. In most applications, including our own, p is assumed to be independent of q and is generated from a Normal(0, 1) density. The leapfrog method then proceeds as follows We iteratively compute the states at t = , 2 , 3 , . . . , L , where L stands for the number of leapfrog steps. Neal (2011) showed that if we take two choices of L and such that L 1 1 = L 2 2 , the solution due to larger L and smaller works better. In our application, we keep L fixed and tune within our MCMC implementation to achieve an acceptance rate between 0.6 to 0.9 which exhibits good mixing and efficiency in posterior sampling. After running the iteration from L-many steps, we obtain an updated q * and p * . The acceptance probability for the new q * is where U (q) stands for the negative log-posterior and K(p) = 1 2 p 2 is the Gaussian kinetic energy term.

S.4 Model Identification
The additive zero mean error formulation of our model in Section 3.3 of the main paper readily establishes the identifiability of f x when the joint and conditional distributions of w, x are all bounded, the characteristic function of (x|w) is non-vanishing everywhere and m i ≥ 3 replicates are available for each subject i (?). Broadly speaking, if the density f x|w varies with x, its characteristic function does not vanish. With sufficient variability of the density of x|w, the observations w have enough information to allow the recovery of the density of x (Sarkar et al., 2018).

S.5 Posterior Computation
We now discuss our MCMC algorithm to draw posterior samples of the model parameters.
We shall provide full conditionals for all the parameters. In addition, we also provide explicit expressions of the derivatives for the parameters that are updated using HMC.

Choice of Hyper-parameters and MCMC Initial Values
The starting values of some of the parameters for the multivariate problem are determined by first fitting the univariate sub-model of Sarkar et al. (2021). We describe the hyperparameter choices and the initial values for the sampler for the marginal univariate models first. Unless otherwise explicitly specified, the prior hyper-parameter choices for similar model components for the multivariate model remain the same as that used for the univariate models. We only detail the sampling steps for the multivariate method. To make the recalls for all the components to be unit free and have shared support, we transformed the recalls as w ,i,j = 20 × w ,i,j max{w ,i,j } . The latent x ,i 's can then be safely assumed to lie in [0, 10], greatly simplifying model specification and hyper-parameter selection.
For the univariate samplers, we then used the subject-specific sample means w ,1:n as the starting values for x ,1:n . The appropriate number of mixture components in a mixture model depends on the flexibility of the component mixture kernels as well as on the specific demands of the particular application at hand. With appropriately chosen mixture kernels, univariate mixture models with 5-10 components have often been found to be sufficiently flexible. Similar observations have also been made for penalized mixtures of splines (?). Detailed guidelines on selecting the number of mixture components for the specific context of deconvolution problems can be found in Sections S.1 and S.6 in the supplementary material of Sarkar et al. (2018). Based on these guidelines and extensive numerical experiments, we used 10 equidistant knot points for the B-splines supported on [A, B] = [0, 10] for modeling the functions s 2 ( x /d) and κ (x ). Thus there are K s = K κ = 13 many basis functions. We used K = 10 mixture components for the truncated normal mixtures modeling their densities. For the Dirichlet prior hyper-parameter, we set α = 1/K. The hyper-parameters for the smoothness-inducing parameters are set to be mildly informative as a ξ = a β = a ϑ = 10, b ϑ = b β = b ξ = 1. Introducing latent mixture component allocation variables c x,1:d,1:n , c ,1:d,1:N and c2 ,1:d,1:N , we can write the univariate sub-model of Sarkar et al. (2021) as . . . , d, and ( ,i,j |c , ,i,j = k, c2 , ,i,j = t, µ , ,k,t , σ 2 , ,k,t ) ∼ Normal( ,i,j |µ , ,k,t , σ 2 , ,k,t ), = 1, . . . , d.
The mixture labels c x, ,i 's, and the component-specific parameters µ x, ,k 's and σ x, ,k 's are initialized by fitting a k-means algorithm with k = K. The parameters of the distribution of scaled errors are initialized at values that correspond to the special standard normal case.
The initial values of the smoothness-inducing parameters are set at σ 2 ϑ, = σ 2 ξ, = σ 2 ξ, = 0.1. The associated mixture labels c , ,i,j 's are thus all initialized at c , ,i,j = 1. The initial values of ϑ 's are obtained by maximizing We now discuss how we set the initial values of the sampler for the multivariate method. The starting values of the x ,i 's were all set at the corresponding estimates returned by the univariate samplers. The initialization of the B-spline coefficients of κ (·)'s and s 2 (·) were done based on the estimated x ,i 's from the univariate sampler as described below. We first normalize w i,j as w i,j = w i,j / w i,j 2 . Similarly, we also obtained x i 's. Then, for each individual i, we assume w i,j ∼ vMF(f i ), where vMF stands for the von Mises-Fisher distribution (Jupp and Mardia, 1979) and obtain maximum likelihood estimates f i = ( f 1,i , . . . , f d,i ) T using the R package movMF (?). Based on these estimates we set κ ,i (x ,i ) = f ,i / x ,i . Similarly for each individual, we calculate the variance of v i = log w i,j 2 x i 2 , j = 1, . . . , m i and set s 2 ( x i 2 /d) = v i . Using the R package nnls (?), we then fit the non-negative least squares to obtain the B-spline coefficients β κ, and β s for cubic splines with 10 knots. For β κ, , the first input in nnls is the matrix M of dimension n × 13, where M ,i,i = B i (x ,i ) and the second input is {κ ,i (x ,i ) : i = 1, . . . , n}. Similarly, to estimate β s , the corresponding first matrix input M s is of dimension n × 13, where M s,i,i = B i ( x i 2 /d) and the vector input is {s 2 ( x i 2 /d) : i = 1, . . . , n}.
We set the number of shared atoms of the mixture models for the densities f x, at K = 10. The shared atoms of the mixtures of truncated normals for the marginal densities f x, are initialized by iteratively sampling them from their posterior full conditionals 500 times, keeping the x ,i 's fixed at their estimated initial values.
Finally, the parameters specifying R x were set at values that correspond to the special case R x = I p . In our sampler for the multivariate problem, we first update the parameters specifying the different marginal densities using a pseudo-likelihood that ignores the contribution of the copula. The parameters characterizing the copula and the latent x i 's are then updated using the exact likelihood function conditionally on the parameters obtained in the first step. The results of ? suggest that such two-stage approach performs just as good as joint estimation procedures, validating their use for computational simplicity. We then update the parameters of the marginal densities again and so forth. Lastly, we update the B-spline coefficients involved in our measurement error model using an HMC sampler.

MCMC Iterations
Our sampler for the multivariate model iterates between the following steps.
1. Updating the parameters specifying f x, , = 1, . . . , d: We modelled the marginal densities of the components = 1, . . . , d using mixtures of truncated normals with shared atoms. We update the mean and variance of the mixture normal components using HMC. The mixing probabilities and mixture indicators are updated using the full conditional conjugate posterior distributions.
Specifically, the full conditional of π x, ,k is given by where n x, (k) = n i=1 1(c x, ,i = k) as before. The full conditional of c x, ,i is given by a standard multinomial. The full conditional of µ x,k is given by which gives the following negative log-likelihood where Φ(·) stands for the standard normal cdf and the derivative is given by Similarly, the full conditional of σ 2 x,k is given by which gives the following negative log-likelihood , and the derivative with respect to σ 2 x,k is given by These parameters are updated by HMC with Gaussian kinetic energy as described in Section S.2 in the supplementary materials.
2. Updating the B-spline coefficients specifying s 2 and κ for = 1, . . . , d: These parameters are also updated using HMC sampler. We provide below the associated negative log-likelihood functions and the corresponding derivatives.
(a) The negative log-likelihood to update the B-spline coefficients β s is And the derivative is (b) The negative log-likelihood to update the B-spline coefficients β κ, is We consider adaptive random-walk MH sampling for β κ, using a multivariate normal proposal with a covariance matrix being set based on the generated posterior samples adaptively.

3.
Updating the values of x: The full conditionals for x i are given by where F x, (x ,i |ζ) = Φ(y x, ,i ). These full conditionals do not have closed forms. Metropolis-Hastings (M-H) steps with independent truncated normal proposals for each component are used within the Gibbs sampler.
4. Updating the parameters specifying the copula: We have F x, (x ,i |ζ) = Φ(y x, ,i ) for all i = 1, . . . , n and = 1, . . . , d. Conditionally on the parameters specifying the marginals, y x,1:d,1:n are thus known quantities. We plug-in these values and use that (y x,i |R x ) ∼ MVN d (0, R x ) to update R x . The full conditionals of the parameters specifying R x do not have closed forms. We use random walk M-H steps to update the d(d − 1)/2 set of polar angles. Since polar angles have bounded support, they are proposals generated from truncated normals, supported according to the respective polar angles. The angles {ζ m,m−1 : m = 1, . . . , d} are supported in [0, 2π], whereas other angles are supported in [0, π]. We discretized the set of possible values of each ζ m,s with M = 41 equidistant grid points covering its support. A new value ζ s,new is proposed at random from the set comprising the current value and its two neighbors. Their proposed values are accepted with probabilities min{1, a(ζ m,s,new )/a(ζ m,s )}, where With carefully chosen initial values and proposal densities for the M-H steps, we were able to achieve quick convergence for the MCMC samplers. We also automatically tune the step lengths for all the Metropolis proposals to achieve acceptance rates within a pre-specified range, namely, 0.4-0.5 for the M-H steps and 0.6-0.9 for HMC. For our proposed method, 5, 000 MCMC iterations were run in each case with the initial 3, 000 iterations discarded as burn-in. The remaining samples were further thinned by a thinning interval of 5. We programmed in R. With n = 1000 subjects and m i = 3 proxies for each subject, on an ordinary laptop, 5, 000 MCMC iterations required approximately 4 hours to run.

S.6 Computational Cost and Complexity
We also determine the time complexity of the proposed method in the increasing d regime via simulations. Specifically, we run the proposed method for three choices of d, namely d = 3, 5, 10, and 30 simulation replications. We repeat the same steps for SPMC as well. For a fair comparison, we run both of the two methods for 2000 MCMC iterations. Our final computation times are averaged over all 30 replications. Subsequently, we fit the model log(t) on log(d) to compute a such that t = O(d a ), where t is the computation time when the data is d-dimensional. The results are shown in Figure S.2. The solution of a turns out to be 1.49 for our method. For SPMC, the order of computation turns out to be 1.51 which is almost the same as our proposed method. However, our proposed method converges faster and requires fewer MCMC steps, likely due to the use of more efficient gradient-based MH proposals for some of the parameters.

S.7 Posterior Consistency
We study the posterior consistency of the proposed model as the sample size increases to infinity. We make following assumption for true distribution of f 0,x with F 0,x, (x ) = Φ(y x, ), where F 0,x, is the cdf corresponding to f 0,x, , and the univariate densities f 0,x, (·)'s are supported in [A, B] such that f 0,x, (·) ∈ F where F is the class of square-integrable functions. Convex combinations of Gaussians are dense in L 2 meaning they can approximate any square-integrable function with arbitrarily small approximation error (Kostantinos, 2000). One may vary the support of f 0,x, with , for example, the support for x may be [A , B ]. However, choosing A = min A and B = max B , the proofs will remain the same with minor technical modifications. We present here the main results. To save space, detailed proofs are deferred to Section S.8 of the supplementary materials. Let θ = {κ (·) : = 1, . . . , d, s 2 (·), K κ , K s } {K, R x , f x, : = 1, . . . , d} = (θ , θ ) denote the complete set of parameters. The marginal density for w is f w (w) = X f w|x (w|x)f x (x)dx.
We have x ∈ X = [A, B] d . Also, let the replicates w ∈ W ⊆ R d . We have F 0,x (X ) = 1 and F x (X ) = 1. Let θ 0 = (θ 0 , θ 0 ) be the true value of θ. We make the following set of assumptions.
For the sake of generality, we also put priors on K κ , K s and K with probability mass functions given by where b κ1 , b κ2 , b s1 , b s2 , b 1 , b 2 > 0 and 0 ≤ b κ3 , b s3 , b 3 ≤ 1 for i = 1, 2. As special cases of the above, we can obtain the Poisson and geometric probability mass functions respectively for b κ3 , b s3 , b 3 = 1 or 0. However, for computational simplicity, we do not consider these priors while fitting the model but propose to tune K κ , K s and K based on the sample size.
Posterior consistency theory studies recovery of the 'true' parameter θ 0 with increasing sample size when the data is sampled from the distribution characterized by θ 0 . Our notion of recovery is based on following L 1 -distance metric d(θ 1 , θ 2 ) = f 1,w − f 2,w = |f 1,w − f 2,w |dw.
If we choose N 1n , M 1n , N 2n , M 2n , N 3n , M 3n , and M 4n as positive polynomial of n depending on the constants in above equations and choose m 2n and m 4n as negative polynomials of n, we are done.  replicates per subject when the true data generating process follows the additive model from Sarkar et al. (2021). We present the estimates from our proposed method (in black) and the method of Sarkar et al. (2021) (in green) along with the truth (in red).