3
$\begingroup$

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$:

  1. 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?

  2. 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?

  3. 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).

|cite|improve this question
$\endgroup$
  • $\begingroup$ When I'm permitted to answer my own question (tomorrow), I'll post brief descriptions of two ad hoc methods I've considered, including the "truncate-and-transform" method used for the example. $\endgroup$ – Adam Hafdahl May 4 '13 at 4:25
  • $\begingroup$ The task of problem #2 is $\min_{\hat\theta_Y} \| f(\hat\theta_Y) - \theta_T \|$. Depending on what the set $\mathcal A$ (and therefore $f$) looks like, we can bring standard optimization techniques to bear here. If $f$ turns out to be convex and smooth, for example, then our optimization problem is convex, we can get a gradient at least through automatic differentiation if not in closed form, and then we can go to town (also solving problem 3 with standard bounds on the accuracy of our solutions). I don't have any great ideas on what kind of $\mathcal A$ would cause a convex $f$, though. $\endgroup$ – Dougal May 4 '13 at 4:47
  • $\begingroup$ I thought that truncation, a.k.a. censoring, of Normal distributions gives rise to problems addressed by probit regression. $\endgroup$ – DWin May 4 '13 at 5:38
  • $\begingroup$ Thanks, @Dougal. The main challenge is that for most sets $\mathcal{A}$ that I'm interested in, I don't know how to compute $f$ but can estimate it (e.g., sample from $Y$ for several values of $\theta_Y$). Could this be handled by stochastic versions of standard optimization techniques or by some type of MCMC approach? I'll look into your remark about $f$'s convexity, which might relate to the issue of invertibility in #1. $\endgroup$ – Adam Hafdahl May 4 '13 at 13:16
  • $\begingroup$ @DWin: I believe truncation is distinct from censoring (e.g., in tobit regression) and dichotomization (e.g., in probit regression), though all are forms of limiting a variable's information based on a cutpoint. Work on limited variables, missing data, or coarsened data might provide insight into my question. $\endgroup$ – Adam Hafdahl May 4 '13 at 14:23
2
$\begingroup$

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}$).

|cite|improve this answer
$\endgroup$
  • $\begingroup$ My PhD thesis was shorter ... though not much. $\endgroup$ – wolfies May 4 '13 at 17:10
  • $\begingroup$ @wolfies: If your thesis contains ideas relevant to my question, could you post its citation -- or just copy-and-paste the whole thing here? $\endgroup$ – Adam Hafdahl May 4 '13 at 18:44
  • 1
    $\begingroup$ Aha - at last - someone with a sense of humour :) $\endgroup$ – wolfies May 4 '13 at 18:45
  • $\begingroup$ This truncate-transfer method seems to work well for a doubly truncated univariate normal, for which $d = 1$ and $Y$ is truncated outside the real interval $\mathcal{A} = (a, b)$: Steps 2 and 3 can use simulation or this case's closed-form expression for $f$. Incidentally, in this case the post-truncation mean and variance have bounds $a < \eta < b$ and $\Omega < (b - a)^2 / 12$; more specifically, $\Omega$'s upper bound is highest when $\eta = (a + b) / 2$ and decreases to 0 as $\eta$ approaches $a$ or $b$. $\endgroup$ – Adam Hafdahl May 6 '13 at 20:01
1
$\begingroup$

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):

  1. 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.

  2. 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$.

  3. 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.

  4. 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$.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ This regression approach estimates $\mathrm{E}(\theta_Y | \theta_T)$, which suggests a Bayesian aproach: Specify a prior for $\theta_Y$, say $p(\theta_Y)$, and together with $\theta_Y$'s likelihood, say $L(\theta_Y | \theta_T) = p(\theta_T | \theta_Y)$, this yields $\theta_Y$'s posterior, say $p(\theta_Y | \theta_T)$. Using notation for the random variables $\Theta_Y$ and $\hat\theta_T$, we could write the posterior more explicitly as $p(\Theta_Y = \theta_Y | \hat\theta_T = \theta_T)$. I'm not sure whether this is valid or how to implement it. $\endgroup$ – Adam Hafdahl May 22 '13 at 21:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.