function rand_indices=rdisc(n,p)



% given vector of probabilities p

% rdiscrete(n,p) simulates n random

% indices



q=cumsum(p);

r=rand(1,n);

rand_indices=ones(1,n);

for i=1:length(q)

  rand_indices=rand_indices+(r>(q(i)*ones(1,n)));

end

