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