# Generating random variables satisfying constraints

I need to generate a list of random variables $\bf{x}$ subject to constraints that can be expressed in the form $\bf{E}x=b$ where $\bf{E}$ is an $m \times n$ matrix if $\bf{x}$ has $n$ entries. In all the cases I'm dealing with, $n >> m$, for example $n$ will be around 14,000 and $m$ will be 50. I'm not sure what method I will use for random sampling, either normal or uniform, it is not clear which is best for the problem I'm trying to solve, but I need all the variables to be sampled from distributions with the same mean and range/variance.

What I've been doing to solve this is reducing $\bf{E}$ to row-echelon form, setting all the variables corresponding to columns to the right of the last pivot to random values, and then solving the remaining square matrix equality.

There is a problem however, to solve the square matrix equality, I subtract the already set values from the right hand side. Unfortunately, the variances add as well, so my last 50 values tend to vary hugely, which is unfortunately unacceptable in this problem.

Is there a better way to do this? I cannot think of a way to fix the current method I am using. I use R.

• Unfortunately, you won't be able to do this, unless you get really lucky with your constraint matrix. For example, consider a constraint matrix with two rows, one nonzero entry in the first row that constrains $x_1 = 0$ and two nonzero entries in the second row that constrain $x_2+x_3 = 1$. Obviously, $x_1$ will have a different mean than at least one of $x_2$ and $x_3$, and unless you set the variances for $x_2$ and $x_3 = 0$, a different variance as well. Jun 16, 2012 at 13:50

This paper and R package completely solved my problem. It uses the Markov Chain Monte Carlo method, which relies on the fact that if you can find an initial solution of the constraint, through linear programming, you can find an arbitrary number of them by using a matrix that when multiplied by $E$, the constraints, gives zero. Read about it here:

http://www.vliz.be/imisdocs/publications/149403.pdf

and here is the package:

http://cran.r-project.org/web/packages/limSolve/index.html

Might seem trivial (and not terribly machine efficient), but consider repeating the process until you get a suitable answer? Preferably only modifying a smaller subset each time.