Can I generate multiple (random) vectors [x1, x2, x3, x4, x5] where all x's are positive numbers such that they all satisfy the linear constraints:



Can I for example generate 100 instances of [x1, x2, x3, x4, x5]? And how?


2 Answers 2


Call $\mathfrak A$ the set $$\mathfrak A=\left\{x\in\mathbb R_+^5;Ax=b \right\}$$ where $$A=\left(\begin{matrix} 1 &2 &6 &4 &2\\ 3 &1 &4 &7 &4\end{matrix} \right)\qquad b={10 \choose 16}$$ It is the intersection of two hyperplanes of $\mathbb R^5$ restricted to the positive quadrant. One can then generate a Markov chain on $\mathfrak A$ as follows:

  1. Start from an arbitrary $x^{(0)}\in\mathfrak A$
  2. For iteration $t$, elect at random three different integers in $\{1,2,3,4,5\}$, $i_1,i_2,i_3$
  3. Generate $y_{i_1},y_{i_2},y_{i_3}$ uniformly in $[0,10]^3$
  4. If there exist a solution to $Ay=b$ with the three above components, switch $x^{(t)}$ to this solution $y^*$, else repeat $x^{(t)}$
  5. Increment $t$ and go to 2.

This Markov chain should have uniform distribution over $\mathfrak A$.

Here is for instance an R code implementation:

#warning: codegolfed!
for(t in 1:1e4){

An alternative is to consider the uniform distribution over $$\begin{cases} x_1+2x_2+6x_3\le 10\\ 3x_1+x_2+4x_3\le16\end{cases}\tag{1}$$ which can be simulated by drawing points uniformly in $$[0,16/3]\times[0,5]\times[0,5/3]$$ and keeping only those satisfying (1).


There are two equations with five unknowns so that in general you can solve only two unknowns. For example, you can solve $x_4$ and $x_5$ from the system. Now you only need to sample $x_1$, $x_2$, and $x_3$. Since they are positive, you can sample from the range as given by @Xi'an: $[0, 16/3] \times [0,5] \times [0, 5/3]$. After that you calculate $x_4$ and $x_5$ and check whether they are positive. If yes, keep that vector, otherwise, sample and calculate again.

EDIT: I have edit the half-close interval into close interval. As required, we have to check whether the sampled values are 0 or positive.

  • $\begingroup$ Yes, what I mean is to check whether they are positive or not. If 0, then the procedure has to be repeated. $\endgroup$
    – TrungDung
    Dec 10, 2020 at 21:33
  • $\begingroup$ You are correct at least in R as R says "runif will not generate either of the extreme values unless max = min or max-min is small compared to min, and in particular not for the default arguments". But is it correct in general, and how do we know that? I mean, mathematically, (0, 1) is not the same as [0,1]. Do I miss something from Uniform(a, b) or Uniform[a, b]. Could you help to clarify? $\endgroup$
    – TrungDung
    Dec 11, 2020 at 9:38
  • 1
    $\begingroup$ As I stated in the first comment, a Uniform$(a,b)$ is the same thing as a Uniform$([a,b]$ since changing the density over a set of measure zero does not modify the associated distribution. $\endgroup$
    – Xi'an
    Dec 11, 2020 at 10:44
  • $\begingroup$ Yes, you are right mathematically. But what will be wrong if for example R sets to sample 0 or 1 from U(0,1)? I do not see any disapproval for R to sample the end points? In other words, I do not see why R has to restrict in (0, 1). Of course, it is not mathematically necessary, but a machine can do, right? $\endgroup$
    – TrungDung
    Dec 11, 2020 at 22:04
  • $\begingroup$ This discussion is quite tangential to the purpose of the original question: please check this thread on stackoverflow. $\endgroup$
    – Xi'an
    Dec 12, 2020 at 9:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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