How can I select a random element from each column of an array, ignoring zeros

I have an array of m*n elements, some of them zeros, for example:

A = [1 0 2 3; 0 1 5 6; 0 0 0 4]

I would like to draw randomly, one single non- zero element from each column. one such vector could be:

[1 1 2 6] for example.

I have tried replacing zeros with nan's and using the datasample function to do this, but it does not ignore nan's.

is there a simple way to do this without a loop?

thanks

Answer 1

[Bruno Luong]

B = A;
B(B==0) = NaN;
B = sort(B,1);
[~,rlast] = max(B,[],1);
rrand = ceil(rand(size(rlast)).*rlast);
rval = B(sub2ind(size(B),rrand,1:size(B,2)))

[Walter Roberson]

Interesting approach, Bruno.

I do notice that the generation process has a small statistical bias, which could potentially grow to be serious for large enough arrays.

There is some discussion of bias in random generation at https://www.mathworks.com/matlabcentral/answers/246999-probability-1-2-with-matlab-random-number-generators, at the end of which I wondered how they take care of bias in randi(); John Chilleri proposed one solution.

[Bruno Luong]

For continuous real number (infinite precision) and real uniform distribution, the question is not applicable, and should be take care by people who work on finite-precision computer (so us).

The issue is not even the difference between ">" or ">=" to generate {0,1}; but I'm much more concern about the gaps between discrete numbers when the value is approaching 0. So we talk about might be some quantification effect when scale RAND up to 1E53, that can make some cross effect with CEIL.

The origin Mersenne Twister generates integer(s) up to some fix number of bits, so the gaps are constant. So RANDI could uses the integer as based engine, and not related to RAND.

IMHO, all that effect, could be cause of concern but still negligible on the intrinsic drawback usage of a pseudo random generator. But I'm no expert on this field.

[Walter Roberson]

randi uses the current random stream, but it is not defined whether it calls rand() or uses bits more directly.

The default stream for rand() generates integers and divides by 2^53 to get the floating point values. That is also the case for all of the other random generators except one, which generates floating point more directly. Most have equal gaps between possible outcomes; I would need to have another look at the one I discussed as being able to generate exact halves as I am not sure that the gaps for it are exactly equal. And the one very different one does have uneven gaps because it generates from all representable from 0 to 1.

[Bruno Luong]

Something I don't get it.

According to this wiki page , MT algo generates numbers MT in [0, 2^w−1], so the discrete sample are even (w=53 according to you).

If

rand = MT/2^53

then rand() generates number in [0,1), this is contradict with the doc and rand, which states the number are in open interval (0,1).

But let assumes RAND just simple divide MT, then to have exactly half probability you should call

r01 = rand < 0.5

If however RAND discard 0 (by rejection) then the number of discrete samples is odd, then to have non bias you have to do

r = rand;
while r==0.5
  r = rand;
end
r01 = r < 0.5

[Walter Roberson]

rand tests for 0 and generates another in that case.

Answer 2

[Walter Roberson]

is there a simple way to do this without a loop?

No. There are ways to do it with disguised loops, but I would not call them simple.

Index = @(V,K) V(K)
OneOfOneCol = @(V) Index(V(find(V)), randi(nnz(V)) );
OneOfEachCol = @(M) arrayfun(@(COL) OneOfOneCol(M(:,COL)), 1:size(M,2));

then

OneOfEachCol(A)

This code will fail if there is a column that has only zeros.

[KALYAN ACHARJYA]

@Walter Sir

No. There are ways to do it with disguised loops, but I would not call them simple.

If I use single for loop for the problem, Is it a major cause?

[Bruno Luong]

Test my answer there is no for-loop there.

[Walter Roberson]

Is it a major cause?

I am not sure what you are asking?

If you are asking whether a for loop is always slower then vectorization, the answer is NO, there are situations where for loops can be faster than vectorization. Even when vectorization can be used to solve a situation, it is common that it involves scanning an array multiple times, or involves temporary variables as large as the original variables.