function M = sample_discrete(prob, r, c)% SAMPLE_DISCRETE Like the built in 'rand', except we draw from a non-uniform discrete distrib.% M = sample_discrete(prob, r, c)%% Example: sample_discrete([0.8 0.2], 1, 10) generates a row vector of 10 random integers from {1,2},% where the prob. of being 1 is 0.8 and the prob of being 2 is 0.2.n = length(prob);if nargin == 1  r = 1; c = 1;elseif nargin == 2  c == r;endR = rand(r, c);M = ones(r, c);cumprob = cumsum(prob(:));if n < r*c  for i = 1:n-1    M = M + (R > cumprob(i));  endelse  % loop over the smaller index - can be much faster if length(prob) >> r*c  cumprob2 = cumprob(1:end-1);  for i=1:r    for j=1:c      M(i,j) = sum(R(i,j) > cumprob2)+1;    end  endend% Slower, even though vectorized%cumprob = reshape(cumsum([0 prob(1:end-1)]), [1 1 n]);%M = sum(R(:,:,ones(n,1)) > cumprob(ones(r,1),ones(c,1),:), 3);% convert using a binning algorithm%M=bindex(R,cumprob);