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