Generating random numbers from 0 - 1 with limit on the sum

22 ビュー (過去 30 日間)
Alexandra
Alexandra 2016 年 7 月 14 日
コメント済み: Alexandra 2016 年 7 月 18 日
Hi,
I am trying to generate values from 0 to 1 for several variables (x1 x2 x3).
These results must assure that x1*a1 + x2*a2 + x3*a3 <= D. a1 a2 a3 and D are other variables I define. How can I do this?
Thank you very much,
  2 件のコメント
Alexandra
Alexandra 2016 年 7 月 15 日
The numbers can't be biased. The x variables must be uniform.

サインインしてコメントする。

採用された回答

John D'Errico
John D'Errico 2016 年 7 月 15 日
編集済み: John D'Errico 2016 年 7 月 15 日
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.
  5 件のコメント
Alexandra
Alexandra 2016 年 7 月 18 日
Hi John, That makes sense, thanks. It was what I meant by the data in "a" generating those results, because of the combinations. Just wanted to make sure I'm seing everything. Statistics is not the hurdle =)

サインインしてコメントする。

その他の回答 (2 件)

Azzi Abdelmalek
Azzi Abdelmalek 2016 年 7 月 14 日
編集済み: Azzi Abdelmalek 2016 年 7 月 14 日
a1=2
a2=3
a3=4
D=5
v=rand(1,2)
x1=v(1)
x2=v(2)
x3=rand*(D-a1*x1-a2*x2)/a3
Check
x1*a1 + x2*a2 + x3*a3<D
  5 件のコメント
Alexandra
Alexandra 2016 年 7 月 15 日
Hi John, thank you so much. The exact number of random variables x is 10.

サインインしてコメントする。


Walter Roberson
Walter Roberson 2016 年 7 月 15 日
If you were looking for exact equality with D, then you could probably adapt the code at https://www.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum to add in a weight vector, a
  1 件のコメント
John D'Errico
John D'Errico 2016 年 7 月 15 日
I recall looking at Roger's code, wanting to modify it for a general linear constraint. That might take some serious thought.

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by