% Takes a vector of weigths (non-negative numbers) and returns a 
% realization of a discrete rv with pmf given by the normalized weights.
% w(i) is the probability that i is selected (after normalizing so 
% that w sums to 1)
% j=discrete_rv(w)


function j=discrete_rv(w)

% Check that w is non-negative
if min(w)<0 
    error('Weights must be non-negative.');
end

% Check that some weight is positive
if sum(w)==0
    error('Some weight must be positive');
end


% Let p be the cumulative distribution function
p=cumsum(w)/sum(w);

% Pick a random number uniformly from [0,1] and store 
% it to temp
temp=rand;

% Choose the last integer j for which p(j)<temp
j=find(p>=temp,1,'first');