Search code examples
matlabrandomdistributionprobability-theory

Draw random numbers from pre-specified probability mass function in Matlab


I have a support (supp_epsilon) and a probability mass function (pr_mass_epsilon) in Matlab, constructed as follows.

supp_epsilon=[0.005 0.01 0.015 0.02]; 

suppsize_epsilon=size(supp_epsilon,2);

pr_mass_epsilon=zeros(suppsize_epsilon,1);

alpha=1;
beta=4;

for j=1:suppsize_epsilon
    pr_mass_epsilon(j)=betacdf(supp_epsilon(j),alpha,beta)/sum(betacdf(supp_epsilon,alpha,beta));
end

Note that the components of pr_mass_epsilon sum up to 1. Now, I want to draw n random numbers from pr_mass_epsilon. How can I do this? I would like a code that works for any suppsize_epsilon.

In other words: I want to randomly draw elements from supp_epsilon, each element with a probability given by pr_mass_epsilon.


Solution

  • 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:

    1. Generate n samples of a uniform random variable in the interval (0,1).
    2. Compare each sample with the distribution function (cumulative sum of mass function).
    3. See in which interval of the distribution function each uniform sample lies.
    4. 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')
    

    enter image description here