So far we’ve had a look at rejection sampling and importance sampling. Here we take a quick look at slice sampling, although rather than implementing it ourselves, we will use the built in Matlab slicesample function.

Using our parameter estimation example, we will use slice sampling to estimate the mean and sigma of some samples from a normal distribution. Here is the code.

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 = @(params)
prod(normpdf(repmat(x,[1 numel(params(1))]),
repmat(params(1), [1 numel(x)])',
repmat(params(2),[1 numel(x)])' ), 1);
%% generate data
global x
N=20;
x = normrnd(true_mean, true_sigma, [N 1]);
%% DO SLICE SAMPLING
initial_samples = [0,1];
nSamples = 10^4;
trace = slicesample(initial_samples, nSamples,
'burnin', 100,
'pdf',likelihood_func);
mean_samples = trace(:,1);
sigma_samples = trace(:,2);
subplot(1,2,1)
plot(mean_samples,sigma_samples,'.')
xlabel('mean')
ylabel('sigma')
hold on
plot(true_mean, true_sigma, 'r.','MarkerSize',5^2)
axis square
subplot(2,2,2)
plot(mean_samples)
ylabel('mean')
subplot(2,2,4)
plot(sigma_samples)
ylabel('sigma')

view raw
sliceSamplingDemo.m
hosted with ❤ by GitHub

Which results in this plot of samples estimating the posterior P(params|data)

slice-sample-result

Leave a comment

Leave a Reply