rejection sampling demo - Matlab

rejectionSamplingDemo.m

%% true probability distribution
true_func = @(x) betapdf(x,1+1,1+10);

%% Do rejection sampling
% create many samples on interval 0-1
x_samples = rand(10^6,1);
% evaluate for each sample
sample_value = true_func(x_samples);
% accept in proportion to highest
max_value = max(sample_value);
accepted = rand(10^6,1) < (sample_value/max_value);
samples = x_samples(accepted);

%% plot
subplot(1,2,1)
x = linspace(0,1,1000);
plot(x, true_func(x) )

subplot(1,2,2)
hist(samples,1000)
proportion_rejected = sum(accepted==0) / numel(x_samples);
title([num2str(proportion_rejected*100) '% of samples rejected'])

% example inspired by http://www.r-bloggers.com/a-simple-explanation-of-rejection-sampling-in-r/