Multiple solutions for linear equations 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:
x1+2.x2+6.x3+4.x4+2.x5=10
3.x1+x2+4.x3+7.x4+4.x5=16
Can I for example generate 100 instances of [x1, x2, x3, x4, x5]? And how?
 A: 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:

*

*Start from an arbitrary $x^{(0)}\in\mathfrak A$

*For iteration $t$, elect at random three different integers in $\{1,2,3,4,5\}$, $i_1,i_2,i_3$

*Generate $y_{i_1},y_{i_2},y_{i_3}$ uniformly in $[0,10]^3$

*If there exist a solution to $Ay=b$ with the three above components, switch $x^{(t)}$ to this solution $y^*$, else repeat $x^{(t)}$

*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!
m=rbind(m<-1:5,o<-m)
A=rbind(c(1,2,6,4,2),c(3,1,4,7,4))
b=c(10,16)
u=c(6,5,2,3,4)
for(t in 1:1e4){
  o[i<-sample(1:5,3)]=y=runif(3,0,u[i])
  p=ifelse(rep(!prod((o[-i]<-solve(B<-A[,-i],
           b-B%*%y))),5),m[length(m[,1]),],o)
  m=rbind(m,p)}

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).
A: 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.
