MATLAB: How to efficiently generate a random integer within a range from an arbitrary probability distribution

Question

I need to generate a random integer within a range from an arbitrary probability distribution, within a loop of 100000 iterations. My implementation works, but I am not sure it is mathematically clean, and it takes forever:

pdf = [ 0.9 0.3 0.003 0.1 0.07 0.0005 0.003 0.15 0.009 0.08 ]; % discrete prob distrib function

cdf = cumsum(pdf); % cumulative distribution function

cdf = cdf / max(cdf); cdf(1) = 0; % normalization

index = ceil(interp1(cdf, [1:numel(pdf)], rand(1)))

Notice that the pdf above is just an example: my actual case is a vector of about 500 numbers.

Here is a different solution, which seems mathematically cleaner, but does not work for my overall problem, and is just as slow as above:

pdf = [ 0.9 0.3 0.003 0.1 0.07 0.0005 0.003 0.15 0.009 0.08 ]; % discrete prob distrib function
cdf = cumsum(pdf); % cumulative distribution function
cdf = cdf - min(cdf); cdf = cdf / max(cdf); % normalization
index = round(interp1(cdf, [1:numel(pdf)], rand(1)))

Is there a more efficient/correct way to do this?

Answer 1

I think that

index = sum(rand()>cdf)+1;

will be much faster than using interp1 as you do, and will give the same result.