function [ v_variates_vec ] = base_attained_Q4( p_vec, n )

% see nutshell "Random Variables" Section 3.2 

rvarray = cumsum(p_vec); % calculates sum of first probability, then sum of first two probilities, then...

unif = rand(1,n); % create uniform sample realization of "u"
numbases = zeros(1,n); % numbases will be our sample realization of v = g(u)
% at end rename to v_variates_vec


inds = find(unif<=rvarray(1)); % find indices for  u <= p_1
numbases(inds) = zeros(1,length(inds)); % set v to zero at these indices
% since v = g(u) = 0 for u <= p_1

% in fact, numbases already initialized to zero: above is intended only to 
% demonstrate procedure which must be generalized below to the other
% outcomes of v

for j = 1:length(p_vec)-1; % need to loop through other outcomes j for v
    inds = find( (?) ? (?)); % find indices of those values of u which map to outcome v=j (i.e., j = g(u))
    numbases(inds) = ?;   % assign proper value of v at these indices
end


v_variates_vec = numbases; % rename to output in signature


end

