Generating a random number from array based on requirements

I am using the code below to allow me to select a single "random" element from an array with each element having its own weighting.

A=[1,2,3,4];
p=[10 20 30 40];
c=cumsum(p);
[~,r]=histc(rand(1,1),[0 c/c(end)]);
r=A(r);

Suppose my requirements now change, and only elements 3 and 4 can be selected from A. Without having to redo the percentages (i.e. changing p array), how can I modify the inputs to the histc function to achieve this?

Thank you

Answer 1

[Rik]

A=[1,2,3,4];
p=[10 20 30 40];
c=cumsum(p);
isAllowed=[3,4];
while true
    [~,r]=histc(rand(1,1),[0 c/c(end)]);
    if ismember(r,isAllowed),break;end
end
disp(r)

     4

[Matlab Beginner]

Thanks for your answer.

I did implement it with a while loop but as I am dealing with a large amount of data this takes very long so was looking for something more optimal.

[Rik]

I don't know another way without recalculating c or skewing the distribution.

[Matlab Beginner]

How would you go about skewing the distribution?

[Walter Roberson]

A = [1, 2, 3, 4];
p = [10 20 30 40];
isAllowed=[3,4];

pt = zeros(size(p));
pt(isAllowed) = p(isAllowed);
c = cumsum(pt);

[~,r]=histc(rand(1,100),[0 c/c(end)]);
r = A(r);

stem(r)

However this breaks the fundamental premise about not having to change the p array.

[Rik]

You could use something like the mod function to wrap the results so they will always index isAllowed, however, that will completely break the distribution. So instead of changing p you would skew the results so much that you're effectively ignoring it.

[Walter Roberson]

A = [1, 2, 3, 4];
p = [10 20 30 40];
isAllowed=[3,4];
c = cumsum(p);

N = 100;

[~,r]=histc(rand(1,N),[0 c/c(end)]);
r = A(r(ismember(r,isAllowed)));

stem(r)

This is the rejection process. Notice that the size of the returned data is only 70% of nominal. If you need a fixed number of samples output then you have to go back and ask for more. This can be pretty expensive -- for example you might be asking to accept only items with a 5% probability, and then to get 100 outputs on average you would need 100/0.05 = 2000 tries.