MATLAB: Random sampling of elements from an array based on a target condition in MATLAB

randomrandom sampling

I have an array (let's call it ElmInfo) of size Nx2 representing a geometry. In that array the element number and element volume are on the column 1 and column 2 respectively. The volume of elements largely vary. The sum of the volume of all elements leads to a value V which can be obtained in MATLAB as:

V=sum(ElmInfo(:,2));

I want to randomly sample elements from the array ElmInfo in such a way that the volume of sampled elements (with no repetition) will lead to a target volume V1. Note: V1 is less than V. So I don't know the number of elements to be sampled. I am giving an example. For a sampling case number of sampled element can be '10' whereas in other sampling number of sampled element can be '15'.

There is no straightforward MATLAB in-built function to meet the target condition. How can I implement the code in MATLAB?

Best Answer

If your process is just to sample until you are within a tolerance of V_1 and you want things to be truly random, how about the following algorithm:

Randomly shuffle the entries of V
Take a cumulative sum
If there is an entry within a certain tolerance of V_1, keep the corresponding entries and exit
Otherwise, repeat

This is a completely naive brute-force approach, which won't work well at all if your tolerance is small. But if your tolerance is generous, the following might do the job well enough, and will actually sample things at random.

n = 300;
v = rand(n, 1);
V1 = 17;
tol = 1e-2;
sample = [];
maxits = 10000;
for count = 1:maxits
    p = randperm(n);
    s = cumsum(v(p));
    k = find(abs(s - V1) < tol);
    if ~isempty(k)
        sample_indices = p(1:k(1));
        sample = v(sample_indices);
        fprintf('sample found after %d iterations\n', count);
        break
    end
end

Related Solutions

MATLAB: Generating dispersed (non-integer) random matrix/array that sums to a particular value

I think you do not understand what you are asking.

randfixedssum indeed produces results that are uniformly sistributed within the sub-set in question. That is, any point in a 10 dimensional space that satisfies the requirements of a fixed sum is equally likely to arise.

However, that does not mean that it is at all probable you would find something that satisfies your goal, of "dispersion".

For example, suppose you were to choose one element that is greater than 0.5? Then the probability that the other 9 elements were ALL small enough that the sum is 1, is pretty low. In the 9 dimensional space that remains, that event would be actually very uncommon.

Thus, you want to generate 10 numbers, all of which lie between 0.05 and 0.9, such that the sum is 1.

Suppose, just suppose that one of the numbers was say, 0.6? Now what are the odds that you can find 9 other numbers that make the total sum exactly 1, but none of them are less than 0.05? SURPRISE! It can never be done.

In fact, if any simgle element was any larger than 0.55 in this example, your goal will never be doable. So if one element is as large as even 0.55+eps, it is mathematically impossible to find 9 numbers, all of which are between 0.05 and 0.9, such that the sum is 0.45-eps.

Next, suppose one element was even as large as 0.5? Just one element that large?

Now the other 9 elements must all be very close to 0.05. What is the probability of that event? Not surprisingly, it is pretty darn small. I can compute the actual probability of such an event to happen if you need. Being too lazy to think at this time of day...

X = randfixedsum(10,10000000,1,0,0.9);
sum(max(X) >= 0.5)
ans =
      195844

So 1.96e5 such events in 1e7. A little under 2% of the time. As expected, a rare event, and that is EXACTLY as it should be.

You ask for dispersion. But you don't seem to understand what dispersion means or what it implies in this context.

If I look at the distribution of the maximum of all 10 elements, I get something that is actually pretty reasonable.

X = randfixedsum(10,10000,1,0.05,0.9);
   Min     0.1207
  1.0%     0.1342
  5.0%     0.1445
 10.0%     0.1524
 25.0%     0.1674
 50.0%     0.1884
 75.0%     0.2167
 90.0%     0.2503
 95.0%     0.2738
 99.0%     0.3143
   Max     0.4039

Most of the time, we get a maximum value that is pretty small in context. And that is because the sample truly is uniformly distributed around the constraint space. One point in that space is equally as likely to arise as any other point. But that does NOT mean that the maximum is ever likely to be larger than 0.55. In fact, that would be an impossible event.

Suppose instead, that we change the way things were generated? Now, instead of requiring that the min be 0.05. Just make it 0. How do the statistics change?

X = randfixedsum(10,10000,1,0,0.9);
   Min     0.1395
  1.0%     0.1681
  5.0%     0.1902
 10.0%      0.205
 25.0%     0.2353
 50.0%     0.2784
 75.0%     0.3359
 90.0%      0.401
 95.0%     0.4492
 99.0%     0.5479
   Max     0.8123

As you now see, the maximum element is now considerably larger. In the same size sample, I once got something as large as 0.8123. There is now much more room for those "dispersed" events to arise.

MATLAB: Generate random coordinates inside a convex polytope

This is actually easier than you think, IF you understand how to solve the necessary sub-problems. I even have code that does it, but I tend not to give it out, because people seem not to understand that toolbox. Interface is everything. Well, not everything, but important as hell. And since I've never had the ambition to create a complete new interface for those tools, they will languish on my list of things to do.

So, lets see if we can solve the problem. Is this a 2-d problem? 3-d? N-d? I'll show you how to do it in 2-d, because I can make pretty plots there that are easy to look at and verify we were successful. The extension to n-d is trivial.

First, create a convex tessellation.

xy = randn(10,2);
CH = convhulln(xy)
CH =
   1     8
   8    10
  10     6
   6     2
   2     7
   7     1

Those are edges, because I'm working in 2-d. n-d is the same idea.

plot(xy(:,1),xy(:,2),'bo');
hold on
plot([xy(CH(:,1),1),xy(CH(:,2),1)]',[xy(CH(:,1),2),xy(CH(:,2),2)]','r-')

Compute a point in the interior of the polygon. Since it is known to be a convex polygon (polyhedron in n-d), mean will suffice.

xycent = mean(xy,1)
xycent =
   0.026363      0.42644
nxy = size(xy,1)
ncent = nxy+1;
xy(ncent,:) = xycent;

now convert the polygon into a triangulation of that convex domain. Each edge will become a triangle, or (n-1)-simplex in n-d. (A k-simplex is defined by k+1 vertices. It may be full volume, or it may live in a k-manifold, embedded in a higher dimensional space.) These are now triangles, living in the 2-d data plane.

tri = [CH,repmat(ncent,ntri,1)]
tri =
   1     8    11
   8    10    11
  10     6    11
   6     2    11
   2     7    11
   7     1    11
figure
plot(xy(:,1),xy(:,2),'bo');
hold on
plot([xy(tri(:,1),1),xy(tri(:,2),1),xy(tri(:,3),1)]',[xy(tri(:,1),2),xy(tri(:,2),2),xy(tri(:,3),2)]','g-')

You can see one new point in that set, the one we added. It is a common vertex for every triangle.

Next, compute the volume of each simplex. In 2-d, it is just the area. This is mot easily done using determinants, with a loop. Yes, you might know that I hate determinants in general. But sometimes they can be useful, and I'm feeling too lazy.

Note that the loop below uses aR2016b feature. Earlier releases will need to use bsxfun.

V = zeros(1,ntri);
for ii = 1:ntri
  V(ii) = abs(det(xy(tri(ii,1:2),:) - xycent));
end
V
V =
       1.7526       1.2067      0.79587      0.94023       1.0889        2.993

So V is a list of areas (or volumes) for each simplex. Actually I needed to divide by factorial(ndim), where ndim is the dimension of the space we are living in to truly compute volume, but I need only something proportional to volume for the next computation. If that bothers you here, you could divide V by 2, but save the cpu cycles, because the next step is to normalize V to sum to 1.

V = V/sum(V)
V =
    0.19968      0.13747     0.090674      0.10712      0.12406      0.34099

All we need are the relative areas of each simplex.

We are almost done. Lets say we wish to generate M points in total.

M = 1000;

For each point, we first need to know which simplex it lives in. That will be proportional to the relative volume of the corresponding simplexes.

[~,~,simpind] = histcounts(rand(M,1),cumsum([0,V]));

So the bigger simplexes will have proportionally more points to fall in them.

FINALLY, generate the desired random points. Yes, I know this seems like a trip to get here, but most of it was me typing explanation and making pictures.

Again, the code that follows will require R2016b. Easy to fix though.

r1 = rand(M,1);
uv = xy(tri(simpind,1),:).*r1 + xy(tri(simpind,2),:).*(1-r1);
r2 = sqrt(rand(M,1));
uv = uv.*r2 + xy(tri(simpind,3),:).*(1-r2);

uv is the complete set of random points. Trust me? Yeah, I know, trust but verify.

figure
plot(xy(:,1),xy(:,2),'bo');
hold on
plot([xy(tri(:,1),1),xy(tri(:,2),1),xy(tri(:,3),1)]',[xy(tri(:,1),2),xy(tri(:,2),2),xy(tri(:,3),2)]','g-')
plot(uv(:,1),uv(:,2),'m.')

As you can see, the sampling is indeed uniform in that domain.

In 3 or more dimensions, just extend what I did as is appropriate. Ask if you have questions, but the extension really is easy. (Ok, there is one part that requires some thought by the programmer, and that is how to extend the very last sampling step into n-d. The trick is to understand why I used sqrt there.) Again, ask if you need the trick.

I suppose one day I should post code on the FEX that does this, but sampling from a convex n-dimensional polyhedron is not something for which I see people clamoring.

Best Answer

Related Solutions

MATLAB: Generating dispersed (non-integer) random matrix/array that sums to a particular value

MATLAB: Generate random coordinates inside a convex polytope

Related Question