Using the Statistics Toolbox
The randsample
function can do that directly:
result = randsample(supp_epsilon, n, true, pr_mass_epsilon);
Without using toolboxes
Manual approach:
- Generate
n
samples of a uniform random variable in the interval (0,1). - Compare each sample with the distribution function (cumulative sum of mass function).
- See in which interval of the distribution function each uniform sample lies.
- Index into the array of possible values
result = supp_epsilon(sum(rand(1,n)>cumsum(pr_mass_epsilon(:)), 1)+1);
For your example, with n=1e6
either of the two approaches gives a histogram similar to this:
histogram(result, 'normalization', 'probability')