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

Hello,

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

Thanks

Answer 1

[John D'Errico]

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.

Answer 2

[Star Strider]

Not difficult:

n = 5;
x = rand(n,3);
x = x./sum(x,2);
Check = sum(x,2);                       % Check To Be Certain They Sum To ‘1’

[John D'Errico]

Note that the results as sampled using this scheme are, while "random", not uniformly distributed. To get a truly uniform sampling, you need to use scheme like that used by Roger Stafford in randfixedsum or my own randFixedLinearCombination tool, both found on the file exchange.

[JR]

Hi John, Your randFixedLinearCombination tool is awesome! However I don't know how I can implement it with my actual problem with which I have a sum of squares, ie, x1^2 + x2^2 + ... + x9^2 = 0.4 and I wanted to generate n different combinations for x1 to x9 such that the sum of their squares is always 0.4.

[John D'Errico]

You did not say you needed a sum of squares! :) I'll add that answer separately.

[Star Strider]

Incorporating your new requirements, a slight variation on my code will do what you want:

n = 5;
x = rand(n,9);
xv = x.^2;
x = sqrt(xv./(sum(xv,2)/0.4));
Check1 = sum(x,2);                      % Check The Sum Of ‘x’
Check2 = sum(x.^2,2);                   % Check To Be Certain The ‘x^2’ Sum To ‘0.4’

[John D'Errico]

Be careful with rand here. Rand is NOT a rotationally symmetric distribution. Rescaling those values will fail to be uniform.

A simple counter-example is gained from just two dimensions.

n = 100000;
x = rand(n,2);
xv = x.^2;
x = sqrt(xv./(sum(xv,2)/0.4));
theta = atan2(x(:,2),x(:,1));
hist(theta,1000)

A rotationally symmetric sample would have produced a flat histogram here. As generated from rand, points at some angles around the quarter circle arise roughly twice as often as others. The highest frequency happens at pi/4 (45 degrees). If you think about it, you can see why that happens. The same thing happens in higher numbers of dimensions.

Answer 3

[John D'Errico]

It always works out this way. Someone asks question A, when they really wanted question B answered. :)

You REALLY want to sample from a set in 9 dimensions, such that the sum of squares of the 9 numbers is 0.4. You could have said that in the first place! Grumble. :)

This is the one case where a sample and normalize scheme can work and be efficient, as long as you sample properly. What you need to recognize is that the sum of squares is equivalent to requiring the points lie on the surface of a 9 dimensional hypersphere.

Why does that help? Because you can use a rotationally symmetric distribution to solve the problem, and then normalize the result.

You did not say if the numbers may be positive or negative, or if all must be positive. That actually is irrelevant, as it turns out. We will sample using randn, which produces iid gaussian samples in any number of dimensions. It is completely valid to rescale them.

n = 10000;
targetsum = 0.4;
x = randn(n,9);
x = sqrt(targetsum)*x./sqrt(sum(x.^2,2));

min(sum(x.^2,2))
ans =
          0.4

max(sum(x.^2,2))
ans =
          0.4

The resulting samples will be uniformly distributed on the surface of the desired 9-dimensional hypersphere.

Now, if your problem requires that x be all positive, the symmetry of the unit normal distribution in n-dimensions saves you. Just take the absolute value of x.

It would NOT have been valid to re-scale a sample taken using rand though. That would NOT be at all a uniform sampling, failing for the same essential reason I showed above. Of course, it is difficult to plot in 9 dimensions to prove it. Ok. Just trust me here.

But, please, please, please don't ask me how to generate samples such that the sum of cubes is a constant.

[JR]

Thanks John! You're a legend. I wanted to establish a code initially that worked for a simple case so that I could build it up from there hence me asking question A so opposed to B. :)

[John D'Errico]

Oh, by the way, the code I wrote is valid for R2017 and beyond. Earlier versions would require bsxfun to do the rescaling.

I'd still need to think about a sum of cubes constraint though, so I'm hoping you don't want that. :)