MATLAB: Find n-number of solutions that satisfy x+y+z=1

equation

Hello,

How could one generate n random points (x,y,z) that satisfy the equation x+y+z=1 ?

Thanks

Best Answer

With only one exception that I can think of off the top of my head, the sample & then normalize solutions that people usually tell you to use are not able to produce a uniform sampling. So, while the result may be "random" it won't uniformly sample the possible sets of numbers that satisfy the constraint. In order to get true uniformity, you need to use a tool designed explicitly to solve that problem.

n = 1e4;
x = rand(n,3);
x = x./sum(x,2);

Now, lets plot the points. Are they truly uniform over the region of interest?

plot3(x(:,1),x(:,2),x(:,3),'.')
view(-45,0)

I've rotated the plot so the point set is in the image plane. Does it look random? NO. In fact, the points are NOT sampled uniformly from that region.

The corners of the triangle are much more sparsely sampled than the rest of the triangular region.

Instead, see how randfixedsum does, on the same problem:

Y = randfixedsum(3,n,1,0,1)';
plot3(Y(:,1),Y(:,2),Y(:,3),'.')
view(-45,0)

As you should see, this point set is quite uniformly sampled over the domain of interest.

Alternatively,

Y = randFixedLinearCombination(n,1,[1 1 1],[0 0 0],[1 1 1]);
plot3(Y(:,1),Y(:,2),Y(:,3),'.')
view(-45,0)

As you can see, both tools are able to produce a nicely uniform sampling.

So use randfixedsum , or my randFixedLinearCombination , which is able to solve a slightly more complex problem, but in a very high number of dimensions, won't be as efficient as the one by Roger. 3 dimensions is not high though, so take your pick. They will both be efficient here, and both nicely uniform.

Related Solutions

MATLAB: Generating random numbers from 0 – 1 with limit on the sum

Ok, lets talk about the problem as defined. It is possible to solve the problem, but things are not always simple in high dimensions. In a high number of dimensions, the solution may be highly computationally intensive.

The set of numbers typically generated by rand are uniform, in 3 dimensions or in n dimensions. In n-dimensions, we can think of the points as living in a hypercube, an n-cube.

But if you require that the weighted sum

x1*a1 + x2*a2 + x3*a3 <= D

be limited to some constant D, then you actually want a UNIFORM random sample that lives in a truncated cube, but with one part of the cube sliced away, as if eaten by rats. :)

So, consider the problem in just 3 dimensions. I'll pick some arbitrary values for a1, a2, a3, D to see what happens here. For example, suppose that we chose the a's (and D) as:

a = [2;3;4];
D = 4;

I'm using a vector here, because numbered variables are just SO bad a programming solution. I won't do it. LEARN TO USE VECTORS. Were this a 10-dimensional problem, this solution would be far more readable and writable. In fact, as you see below, to change the solution to 10 dimensions from 3 dimensions, the change was trivial.

I've also used a column vector for a. So if X is a vector, of length 3 that comes from rand:

X = rand(1,3);

Well, actually, lets generate a large sample, say 1000 such samples, in 3 dimensions.

N =1000;
X = rand(N,3);

So now each row of X may be interpreted as one point in the unit 3-cube. We could plot them in 3-d, but many of them are outside of the linear constraint plane. In fact, only 353 of the 1000 samples we generated were within the constraint.

sum(X*a < D)
ans =
   353

Time for a picture.

ind = X*a < D;
plot3(X(ind,1),X(ind,2),X(ind,3),'ro')
hold on
plot3(X(~ind,1),X(~ind,2),X(~ind,3),'g+')
box on
view(42.5,19.6)

Ok, so it looks like roughly 35% of the 1000 points lie underneath the constraint plane. In fact, if D were a little smaller (in this case with these values for a, no larger then 2) then the region of interest would be described by a tetrahedron. But I chose D=4 above. In fact, it is possible to create the region of interest as a solid in 3 dimensions.

So it is a somewhat complex solid. Not just a simple tetrahedron. In fact, the code that I used to generate the picture created 14 different tetrahedrae in the simplicial complex that describes the volume. You can see the edges of them in that plot, where they touch the surface. Ok, better, using transparency...

So a bit messy.

Regardless ... the simple solution to this problem, to generate perhaps 1000 points that lie in the required region, is to just generate too many random samples, then throw away those that lie beyond the constraint.

So, lets suppose that we desire to sample N points, randomly and uniformly from the above volume, so that no point exceeds the required limit D. Simple (sort of) is a rejection sampling scheme:

a = [2;3;4];
D = 4;
N = 1000;
Xfinal = zeros(N,3);
Nfound = 0;
Nneed = N - Nfound;
while Nfound < N
  X = rand(N,3);
  % rejection step
  ind = find(X*a <= D);
  if ~isempty(ind)
    if numel(ind) > Nneed
      ind = ind(1:Nneed);
    end
    Xfinal(Nfound + (1:numel(ind)),:) = X(ind,:);
    Nfound = Nfound + numel(ind);
    Nneed = N - Nfound;
  end
end

Here, I've plotted the final set of 1000 points, with the required volume overlaid in transparent form to see that indeed, they do satisfy the constraint.

The nice thing is, the above solution is rather easily converted to work in a higher number of dimensions. The downside is that in a higher number of dimensions, the rejection step may be come slow. The curse of dimensionality nails you here. There are other, more sophisticated solutions that can be employed, but even there, they get nasty in a high number of dimensions.

MATLAB: Random positive numbers with constant sum

First of all, there is no such thing as a random positive number. If you don't define the distribution, then it is meaningless to talk about a random number.

And you can't have a uniform random number on the interval (0,inf].

You CAN define a uniform random number on A BOUNDED SET. So, the interval[0,1]. But if you want to constrain the sum to be constant, then be careful. This question gets asked over and over again, and people answer it thinking you can solve it trivially. They generate a vector of numbers, and then scale them so the sum is the desired value. WRONG. This generates a distribution that is not in fact uniform.

So use a tool like randfixedsum , as found on the file exchange, or my own randFixedLinearCombination , also on the file exchange.

The latter tool is more capable for some problems, but for long vectors that sum to a constant, it will be inefficient. It probably won't work well in more than about 6 dimensions. But randfixedsum is fine beyond that point.

So to generate a set of 10 numbers that sum to 4, where each element in the vector is in the interval [.2, .6], is as simple as:

S = randfixedsum(10,1,4,.2,.6)
S =
    0.4354
    0.5535
    0.2221
    0.3235
    0.4106
    0.5368
    0.3963
    0.2636
    0.4705
    0.3877
sum(S)
ans =
    4.0000

Again, DON'T use a rescaling scheme. Randfixedsum is perfect for the problem.

Best Answer

Related Solutions

MATLAB: Generating random numbers from 0 – 1 with limit on the sum

MATLAB: Random positive numbers with constant sum

Related Question