% probability density function is based on this data
probabilities = [0, 1, 10, 2, 0, 0, 4, 5, 3, 1, 0];
x = 1 : length(probabilities);

% do spline interpolation for smoothing
xq = 1 : 0.05 : length(probabilities);
pdf = interp1(x, probabilities, xq, 'spline');

% remove negative elements due to the spline interpolation
pdf(pdf < 0) = 0;

% normalize the function to have an area of 1.0 under it
pdf = pdf / sum(pdf);

% the integral of PDF is the cumulative distribution function
cdf = cumsum(pdf);
% these two variables holds and describes the CDF
xq;         % x
cdf;        % P(x)

% remove non-unique elements
[cdf, mask] = unique(cdf);
xq = xq(mask);
% create an array of 2500 random numbers
randomValues = rand(1, 2500);

% inverse interpolation to achieve P(x) -> x projection of the random values
projection = interp1(cdf, xq, randomValues);
hist(projection, 50);
