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