My problem concerns the estimation of truncated multivariate normal distributions under constraints.
Let $X_1$ and $X_2$ two random variables following normal distributions $\mathcal{N_1}(m_1,\sigma_1^{2}) $ and $\mathcal{N_2}(m_2,\sigma_2^{2}) $ respectively.
The estimation of $X_1, X_2$ over a polyhedra is usually done with a Gibbs sampler.
My question is : How to proceed to estimate $X_1, X_2$ over a polyhedra when you add the constraint $X_1 = X_2$ ?
Thank you.
A general approach is explained here. Read section "3.2 Sampling from a constrained Gaussian". A quote from the paper:
A similar problem has previously been treated by Geweke [2], who proposes a Gibbs sampling procedure, that does not handle equality constraints and no more than N inequality constraints. Rodriguez-Yam et al. [6] extends the method in [2] to an arbitrary number of inequality constraints, but do not provide an algorithm for handling equality constraints. Here, we present a general Gibbs sampling procedure that handles any number of equality and inequality constraints.
Anyway, it is not simple literature.
An other, less involved way is to use rejection/MCMC sampling with a proposal algorithm that gives $(X_1, X_2)$ in the constrained space. A simple one would be $X_1 \sim \mathcal{U}(0,1)$ and $X_2 = X_1$. Under such space the denisty $f(x_1,x_2) \propto \mathcal{N}(x_1,x_2|\mu, \Sigma)$, which you can use within your sampling algorithm as an acceptance/rejection criterion. Note that $\mathcal{N}(x_1,x_2|\mu, \Sigma)$ under such space might be very (!!) small, so small in fact that it might result in zero over the whole space on your computer. So you might have to take advantage of the division $f(x_1,x_2)/f(x_1^*,x_2^*)$ in Metropolis-Hastings algorithm which in this case results in $exp(a - b^*)$.
I think the difficulty is worse than just finding the right algorithm. The first matter of business is deciding which conditional distribution you want to draw from, because they are non-unique.
I'm going to ignore the restricted support, which is a red herring. Instead, consider a bivariate normal distribution with correlation $\rho$: $$(X_1, X_2)^t \sim \mbox{N}((0,0)^t,\Sigma_{\rho}).$$
Now suppose we want to sample from this distribution, restricted to the line $X_1 = X_2$. The measure on this set is non-unique. It will be different depending on if you condition on $Z = 1$, with $Z = X_1/X_2$ versus $Z=0$, with $Z = X_1 - X_2$. I asked about this in my question here (and provide some references, unfortunately behind a pay-wall).
You start out with independent Gaussian random variables, but one you restrict to some rectangle the random variables are no longer independent and the issue above arises. I think it is an interesting question as to how Gibbs samplers relate to particular choices of conditioning variable ($Z$ in my example). For example, you can imagine drawing $X_1$ from its marginal distribution on your rectangle and then just setting $X_2$ to whatever value you have drawn. But this is obviously different than drawing $X_2$ from its marginal distribution and setting $X_1$ to whatever value you have drawn. I don't have a good sense for how these two methods relate to the formal definition of a conditioning sigma algebra. It would be useful to know about because it is often easier to think in terms of Gibbs samplers than it is to think about formal limits.
