Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
rejection sampling demo - Matlab - parameter estimation
true_mean = 0;
true_sigma = 1;
% likelihood_func = @(x, mean, sigma) normpdf(x, mean, sigma);
% the above function to calcalate in matrix form, for speed
likelihood_func = @(x, mean, sigma)...
prod(normpdf(repmat(x,[1 numel(mean)]),...
repmat(mean, [1 numel(x)])',...
repmat(sigma,[1 numel(x)])' ), 1);
%% generate data
N=20;
observed_data = normrnd(true_mean, true_sigma, [N 1]);
%% Do rejection sampling
% create many samples for mean and sigma
n_samples = 10^6;
mean_samples = (rand(n_samples,1)-0.5)*5;
sigma_samples = rand(n_samples, 1) * 10;
% evaluate likelihood for each (mean, sigma) sample
sample_value = likelihood_func(observed_data, mean_samples, sigma_samples);
% accept in proportion to highest
max_value = max(sample_value);
accepted = rand(1,n_samples) < (sample_value./max_value);
mean_samples_accepted = mean_samples(accepted);
sigma_samples_accepted = sigma_samples(accepted);
%% plot
subplot(1,2,1)
plot(mean_samples,sigma_samples,'.')
xlabel('mean')
ylabel('sigma')
title('initial samples')
subplot(1,2,2)
plot(mean_samples_accepted,sigma_samples_accepted,'.')
hold on
plot(true_mean, true_sigma, 'r.','MarkerSize',5^2)
xlabel('mean')
ylabel('sigma')
title(['posterior samples, ' num2str(sum(accepted==0)/n_samples*100)...
'% rejection rate'])
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment