rejection sampling demo - Matlab - parameter estimation

rejectionSamplingDemo2.m

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'])