# Gibbs sampler with linear constraints

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.

• What do you mean by estimation over a polyhedron? – Noah Stein Jun 11 '14 at 16:16
• I want to sample multivariate normal distributions over a polyhedra, for example over a rectangle. This problem is handle with the Gibbs sampler and gives truncated multivariate normal distributions. Hower, i don't know how to sample these distributions when are added general constraints bewteen the random variables. (like in the example.) – user41037 Jun 11 '14 at 16:23
• When you say sample over a polyhedron you mean condition on the sample being in that polyhedron? When you say constrain $X_1 = X_2$ you mean further condition on that event? Off-the-shelf Gibbs sampling will not handle such hard constraints and will perform very poorly (mix slowly) for softened versions of these. But given that you are working with just two univariate normals, it should be possible to work things out analytically. Am I correct in understanding that (unconditionally on the polyhedron and desired equality) $X_1$ and $X_2$ are independent? – Noah Stein Jun 11 '14 at 18:23
• Very interesting, but i was wondering what happened with more random variables and more constraints ? For example : Consider the random variables : $X_1 \sim \mathcal{N_1}(m_1,\sigma_1^{2})$ $X_2 \sim \mathcal{N_2}(m_2,\sigma_2^{2})$ $X_3 \sim \mathcal{N_3}(m_3,\sigma_3^{2})$ $X_4 \sim \mathcal{N_4}(m_4,\sigma_4^{2})$ with the constraints $X_1 = X_2+X_3$ and $X_3=X_4$. – user41037 Jun 11 '14 at 20:37

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 , who proposes a Gibbs sampling procedure, that does not handle equality constraints and no more than N inequality constraints. Rodriguez-Yam et al.  extends the method in  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.