# Sampling point uniformly at random satisfying equality constraints

First of all, I apologize in advance if the question has already been asked in some way on this site and/or if there is a widely known solution to this problem.

The description of my problem is the following:

Sample a point $x$ uniformly at random from the solutions satisfying linear equalities $Ax = b$.

Specifically, I'm interested in sampling a $n$ x $m$ joint probability distribution of two discrete random variables that satisfies given marginal probabilities. This is a special case of the more general described above.

Also, this problem is similar in sampling contingency tables with fixed margins, but it is not restricted to integral solutions.

What I have found / tried so far:

Uniformly Sampling from Convex Polytopes (not exactly what I want, as I'm only interested about on the boundaries and not in the polytope)

Random Sampling a linearly constrained region in n-dimensions... (This is closer to what I want, but can (?) only be applied to the case of $n$ x $2$ dimensional probability distributions)

http://www.mzlabs.com/JMPubs/Sampling%20Contingency%20Tables-Dyer.pdf (Paper about sampling contingency tables with fixed margins)

http://www.mathworks.com/matlabcentral/fileexchange/34208-uniform-distribution-over-a-convex-polytope/content/cprnd.m (Implementation of a method to sample a point that satisfies a set of linear inequalities. I tried enforcing equality constraints but the method does not work)

## 1 Answer

The equality constraints determine an affine subspace of $\mathbb R^n$. After suitable change of variables, this becomes $\mathbb R^m$ for suitable dimension $m$, and the nonnegativity of the entries of your probability distribution gives you a convex polytope. So it is exactly the case of uniformly sampling from convex polytopes.

• Thanks for the reply. It turns out that only (n-1)*(m-1) values are needed to fully determine the probability table. So, it suffices to constrain the first (n-1) and (m-1) values of each row and column respectively to be less than the given marginals and also constrain all values to be greater than 0. The remaining values can then be easily calculated. – George B. Jan 22 '15 at 1:59