Sum of a 4D matrix

조회 수: 10 (최근 30일)
Susan
Susan 2019년 4월 17일
댓글: Susan 2019년 4월 18일
Hey Guys!
I want to generate bunch of matrices in which each element takes value betwwn 0 and 1 and the summation of all elements in each matrix be less than or equal to one. I wrote the following code but it seems it takes forever.
Could someone please take a look at this code and let me know: 1) if the code is written correctly or not 2) Is it possible to rewrite the code in such a way that I have a faster one. Thanks in advance
K = 1;
nbrOfSubCarriers = 2*ones(1, K);
nt_L = 2;
nr_L = 2;
ulim=1; %max value of each element
llim=0; %min value of each element
sumlim=1; %max sum for each matrix
c = zeros(nr_L,nt_L, K, max(nbrOfSubCarriers(:)));
for k = 1 : K
for m = 1 : max(nbrOfSubCarriers(:))
RMat=rand(nr_L,nt_L)*(ulim-llim)+llim;
Rsum=sum(RMat(:));
Rcheck=Rsum>sumlim;
while any(Rcheck)
I=find(Rcheck);
RMat(I,:)=rand(length(I),nt_L)*(ulim-llim)+llim;
Rsum=sum(RMat(:));
Rcheck=Rsum>sumlim;
end
end
c(:,:,k,m) = RMat;
end
  댓글 수: 5
Walter Roberson
Walter Roberson 2019년 4월 17일
My tests suggest that the sum of N random variables is less than or equal to 1, 1/factorial(N) of the time. With you summing nr_L by nt_L all into one sum, then only 1/factorial(nr_L * nt_L) of the generated values will pass the test -- about 1/24 with the values or nr_L and nt_L that you use.
Perhaps it would be acceptable for your purpose to use
if Rsum <= sumlim
accept RMat as-is
else
Rmat = Rmat ./ Rsum; %makes the sum exactly 1
end
Susan
Susan 2019년 4월 17일
Wow. Good to know. Thanks. Let me try that then.

댓글을 달려면 로그인하십시오.

채택된 답변

John D'Errico
John D'Errico 2019년 4월 17일
편집: John D'Errico 2019년 4월 17일
The problem is a relatively easy one, as long as the total sum does not exceed 1. It gets nasty above that.
I'll look at the problem in TWO dimensions, that is, imagine I want to find random numbers that lie within the unit square [0,1]X[0,1], such that the sum, x1 + x2 does not exceed 1.
The set of viable solutions lies within a triangle, the triangle with vertices [0, 0], [0,1], and [1,0]. So an isosceles triangle, with legs of length 1.
This is a different problem than what randfixedsum solves, in a subtle way. It solves the problem where the sum is constrained to be exactly 1. So you cannot use randfixedsum directly. (Though you CAN use it, with a subtle tweak, as long as the maximum on the sum does not exceed 1. That upper limit is important.)
Your admissable set is shown in the figure:
fill([0 0 1 0],[0 1 0 0],'r')
We want to see that the 2-d case is just a simple variation of the n-d case. In n-dimensions, 3 for example, the admissable set becomes an isosceles tetrahedron, so a 3-simplex with legs connecting the four points [0,0,0], [0,0,1], [0,1,0], and [1,0,0]. The number of dimensions is just the number of total elements in your vectors. So, if your vector is of length 100, then the points live in a 100 dimensional space. And then the admissable set will live in a 100-simplex, composed of 101 vertices. Hard to visualize, but a basic extension of the red triangle.
So that is your goal. Again, randfixedsum does not solve that specific problem. But it is close, in a sense.
How can we solve your problem? The idea is we can start with a set that sums to 1. Then take a nonlinear combination of those points with the point at the origin. And since the origin is at zero, it is simple. In two-dimensions, this reduces to:
ndim = 2;
nsample = 10000;
S = randfixedsum(ndim,nsample,1,0,1);
p = rand(1,nsample).^(1/ndim);
S = S.*p;
plot(S(1,:),S(2,:),'.')
Note that this will fail miserably if the maximum sum is greater than 1.
Be careful now. Because the first thing you might do is look at the distribution of the sum of those numbers. Even in the 2-dimensional case, the expected median value for the sum will be
1/sqrt(2)
ans =
0.70711
But in 100 dimensions, the median of the sums will be quite close to 1.
ndim = 100;
nsample = 10000;
S = randfixedsum(ndim,nsample,1,0,1);
p = rand(1,nsample).^(1/ndim);
S = S.*p;
median(sum(S,1))
ans =
0.99305
These points are still uniformly distributed in theory in the domain in question. But you need to understand that the probability of finding a point inside the isosceles 100-hypersimplex with leg length 0.5 is VERY small. That would be:
1/2^100
ans =
7.8886e-31
In fact, what is the the probability of finding a point in that uniformly distributed simplex, such that NONE of the variables x1,...,x100 exceed 0.99?
(0.99)^100
ans =
0.36603
That happens only about 37% of the time.
So don't compute the sum of the samples, and be surprised or disappointed. It will get even worse in a much higher number of dimensions.
Oh, and again, things go to hell if that sum exceeds 1, because then the domain that contains the points is not a single hyper simplex, but much more complex.
  댓글 수: 1
Susan
Susan 2019년 4월 18일
John, thank you very much for your detailed and clear explanation. it was very helpful and understandable. I appreciate your time.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Walter Roberson
Walter Roberson 2019년 4월 17일
  댓글 수: 4
John D'Errico
John D'Errico 2019년 4월 17일
Let me see what I can write up.
Susan
Susan 2019년 4월 17일
Thanks John. Looking forward to hearing from you then :)

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Linear Algebra에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by