Pre-truncation moments for truncated multivariate normal Suppose the random variable $Y$ has a multivariate normal (MVN) distribution, and consider truncating $Y$ in some way to create $T$.  Given $T$'s mean and covariance matrix, I'd like to obtain $Y$'s mean and covariance matrix.
Below I state the above problem more precisely, ask three specific questions, give an example, and close with several remarks.
Suppose $Y$ is a random variable whose $d$-variate normal ($d$VN) distribution has mean $\mu = \mathrm{E}(Y)$ and covariance matrix $\Sigma = \mathrm{Cov}(Y)$, so $Y \sim \mathcal{N}_d(\mu, \Sigma)$.  Let $T = Y | (Y \in \mathcal{A})$ be a truncated version of $Y$ that takes only values of $Y$ in the set $\mathcal{A}$ and has mean $\eta = \mathrm{E}(T)$ and covariance matrix $\Omega = \mathrm{Cov}(T)$.  For instance, in the example below  $\mathcal{A}$ consists of vectors that form valid correlation matrices.  Collect $\mu$'s and $\Sigma$'s $d(d+3) / 2$ distinct parameters in the vector $\theta_Y$, and collect their counterparts from $\eta$ and $\Omega$ in $\theta_T$.  Here are three questions about obtaining $\theta_Y$ for a specified $\theta_T$:


*

*Under what conditions does a unique $\theta_Y$ exist for a specified $\theta_T$?  This may be a question of invertibility: If $f$ is a function such that $\theta_T = f(\theta_Y)$, under what conditions is $f$ invertible?

*In a condition from #1 above, how can we obtain $\theta_Y$ for a specified $\theta_T$, either exactly or by numerical approximation or simulation?  For instance, can this be framed as a problem of optimization or solving a nonlinear system of equations?

*If we approximate or simulate $\theta_Y$ in #2 above, how can we control the approximation's or estimate's accuracy?  For instance, we might want an estimate of $\theta_Y$ to be near $\theta_Y$ with high probability.
EXAMPLE: Consider a fairly simple version of the above problem where $d = 3$, so $Y \sim \mathcal{N}_3(\mu, \Sigma)$, and  $\mathcal{A}$ used for truncation is the set of all 3-vectors containing Fisher-z transformations of a valid Pearson-r correlation matrix's distinct elements (e.g., below the diagonal in row-major order).  For instance,  $\mathcal{A}$ contains [1.20 1.20 0.50] but not [1.20 1.20 0.40].  Let's specify $T$'s desired mean and covariance matrix -- with commas separating rows -- as follows, along with $\theta_T$ (e.g., $\eta$ followed by $\Omega$'s lower triangle in row-major order):


*

*$\eta$ = [0.40 0.60 0.80]

*$\Omega$ = [0.25 0.12 0.06, 0.12 0.16 0.06, 0.06 0.06 0.09]

*$\theta_T$ = [0.40 0.60 0.80 0.25 0.12 0.16 0.06 0.06 0.09]
An ad hoc method -- applying an iterative truncate-and-transform technique each of several large 3VN samples -- yields the following estimates of $Y$'s mean and covariance matrix and, hence, $\theta_Y$:


*

*$\hat \mu$ = [0.329 0.631 0.832]

*$\hat \Sigma$ = [0.292 0.101 0.043, 0.101 0.163 0.064, 0.043 0.064 0.095]

*$\hat \theta_Y$ = [0.329 0.631 0.832 0.292 0.101 0.163 0.043 0.064 0.095]
These estimates are meant to be within about 0.1% of their estimands (in $\theta_Y$) with probability at least .95.
REMARKS: Answering the above questions would probably be easier if $f$ in #1 above were available in closed form, but that's difficult for some truncation rules of interest (i.e., some sets  $\mathcal{A}$).  For now I'd mainly like an accurate solution to these questions, but I'd eventually like to create a fast implementation.  I'm also interested in extensions of this problem, such as working with a mixture of $d$VNs for $Y$ or finding $\mu$ and $\Sigma$ for specified $\mathrm{E}(G)$ and $\mathrm{Cov}(G)$, where $G = g(T)$ is a nonlinear function of $T$.  These questions arose while developing meta-analytic methods for functions of multivariate effect sizes, where between-studies heterogeneity in an effect-size parameter (e.g., correlation matrix) induces heterogeneity in a function of that parameter (e.g., path-model coefficients).  
 A: This truncate-transform method is the first of two ad hoc methods I've contrived.  I'd appreciate constructive criticism, questions, or other feedback about it to help improve it or avoid inappropriate use, and I'm certainly interested in better, more principled strategies -- preferably posted as other answers.
Before using this method, we specify $\eta$ and $\Omega$ (as described in the original question).  The gist of this method is to start with a large sample from a candidate distribution for $Y$ (i.e., based on guesses of $\mu$ and $\Sigma$) and transform this sample based on how well its truncated counterpart -- from a candidate distribution for $T$ -- matches the specified $\eta$ and $\Omega$.  In particular, we apply to the $Y$ candidate whatever linear transformation would make the $T$ candidate's mean $\eta$ and covariance matrix $\Omega$.  Upon convergence, this truncate-transform cycle yields estimates of $\mu$ and $\Sigma$.  Repeating this for several independent samples -- $y_0$ in Step 1 -- permits quantifying and controlling sampling error in the estimates.
Here's a more explicit procedure, omitting some details (e.g., choice of starting values, termination criteria):
1. Initialize: Choose initial values for $\mu$ and $\Sigma$, say $\mu_0$ and $\Sigma_0$, and draw a sample of size $N$ from $Y_0 \sim \mathcal{N}_d(\mu_0, \Sigma_0)$; call this $d$VN sample $y_0$, a $N \times d$ matrix.
2. Truncate: Retain $y_j$'s values that are in $\mathcal{A}$ as the truncated sample $t_j$, and compute $t_j$'s mean and covariance matrix; call these estimates $\eta_j$ and $\Omega_j$.  This amounts to sampling from $T_j = Y_j | (Y_j \in \mathcal{A})$.
3. Transform: Update the $d$VN sample using the linear transformation $y_{j+1} = (y_j - \mathbf{1}_N\otimes\eta_j)\mathbf{U}_j^{-1}\mathbf{U} + \mathbf{1}_N\otimes\eta$, where $\mathbf{1}_N\otimes\mathbf{v}$ is $N$ copies of $\mathbf{v}$ stacked vertically, and $\mathbf{U}_j$ and $\mathbf{U}$ are the upper triangular matrices from Cholesky decompositions of $\Omega_j$ and $\Omega$, respectively.  This amounts to sampling from $Y_{j+1} = (Y_j - \eta_j)\mathbf{U}_j^{-1}\mathbf{U} + \eta$.
4. Iterate: Repeat Steps 2 and 3 until $\eta_j$ and $\Omega_j$ are sufficiently near $\eta$ and $\Omega$, $\mu_{j+1}$ and $\Sigma_{j+1}$ are sufficiently near $\mu_j$ and $\Sigma_j$, or some combination of these.
I don't understand yet why this works or the conditions under which it works.  It seems to yield good estimates of $\mu$ and $\Sigma$ in some situations (e.g., mild truncation of $Y$, valid correlation matrices as $\mathcal{A}$) but not in others (e.g., severe truncation of $Y$, valid path-model coefficients as $\mathcal{A}$).
A: This multivariate-regression method is the second of two ad hoc methods I've contrived.  Again, I'd appreciate constructive criticism, questions, or other feedback about this method to help improve it or avoid inappropriate use, and I'm certainly interested in better, more principled strategies -- preferably posted as other answers.
Before using this method, we specify $\eta$ and $\Omega$ (as described in the original question), whose distinct elements can be collected in $\theta_T$.  The gist of this method is to draw several candidate values of $\theta_Y$, estimate $\theta_T$ for every such value using the same large $d$VN sample, and estimate $\theta_Y$ from a multivariate regression of $\theta_Y$ on $\theta_T$.  I think this essentially entails estimating an approximation to the inverse of the function $f$ mentioned in the original question.  Repeating it for several independent samples -- $z$ in Step 1 -- permits quantifying and controlling sampling error in the estimates.
Here's a more explicit procedure (e.g., choice of starting values, termination criteria):


*

*Draw a standard $d$VN sample of size $N$ from $Z \sim \mathcal{N}_d(\mathbf{0}, \mathbf{I}_d)$; call this $z$, a $N \times d$ matrix.

*Choose $M$ candidate values of $\theta_Y$, of which $\theta_{Y_i}$ contains distinct elements of $\mu_i$ and $\Sigma_i$, $i = 1, 2, \ldots, M$.

*Estimate $\theta_{T_i}$ from $\theta_{Y_i}$ by computing $y_i = z\mathbf{U}_i + \mu_i$, where $\mathbf{U}_i$ is the upper triangular matrix from a Cholesky decomposition of $\Sigma_i$; retaining values of $y_i$ in $\mathcal{A}$ as $t_i$; and computing $t_i$'s mean and covariance matrix.

*Regress the $M$ values of $\theta_Y$ on their $\theta_T$ counterparts, and use this regression's results to estimate $\theta_Y$ from the specified value of $\theta_T$.
This method requires several decisions, such as choosing values of $\theta_Y$ optimally (maybe iteratively, like in response-surface methodology), possibly transforming parameters for $\theta_Y$ or $\theta_T$ (e.g., $\Sigma$ or $\Omega$ as log variances and Fisher-z correlations), and specifying an appropriate regression model (e.g., low-order polynomial) or smoother.  I don't understand well how these choices or other factors affect performance.  The method seems to yield good estimates of $\mu$ and $\Sigma$ in some situations (e.g., mild truncation of $Y$, valid correlation matrices as $\mathcal{A}$) but not in others (e.g., severe truncation of $Y$, valid path-model coefficients as $\mathcal{A}$).
This method might be appropriate for more extensions of the original question than the truncate-transform method I described in a previous answer.  For instance, with this method we could specify parameters of a nonlinear function of $T$, say $G = g(T)$, instead of $\theta_T$, such as by specifying $\theta_G$ as distinct elements of $\mathrm{E}(G)$ and $\mathrm{Cov}(G)$, or we could specify a set of $d(d+3) / 2$ independent parameters from $\theta_Y$, $\theta_T$, and $\theta_G$.
