mhsample

Metropolis-Hastings sample

Syntax

smpl = mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf, 'proprnd',proprnd)
smpl = mhsample(...,'symmetric',sym)
smpl = mhsample(...,'burnin',K)
smpl = mhsample(...,'thin',m)
smpl = mhsample(...,'nchain',n)
[smpl,accept] = mhsample(...)


Description

smpl = mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf, 'proprnd',proprnd) draws nsamples random samples from a target stationary distribution pdf using the Metropolis-Hastings algorithm.

start is a row vector containing the start value of the Markov Chain, nsamples is an integer specifying the number of samples to be generated, and pdf, proppdf, and proprnd are function handles created using @. proppdf defines the proposal distribution density, and proprnd defines the random number generator for the proposal distribution. pdf and proprnd take one argument as an input with the same type and size as start. proppdf takes two arguments as inputs with the same type and size as start.

smpl is a column vector or matrix containing the samples. If the log density function is preferred, 'pdf' and 'proppdf' can be replaced with 'logpdf' and 'logproppdf'. The density functions used in Metropolis-Hastings algorithm are not necessarily normalized.

The proposal distribution q(x,y) gives the probability density for choosing x as the next point when y is the current point. It is sometimes written as q(x|y).

If the proppdf or logproppdf satisfies q(x,y) = q(y,x), that is, the proposal distribution is symmetric, mhsample implements Random Walk Metropolis-Hastings sampling. If the proppdf or logproppdf satisfies q(x,y) = q(x), that is, the proposal distribution is independent of current values, mhsample implements Independent Metropolis-Hastings sampling.

smpl = mhsample(...,'symmetric',sym) draws nsamples random samples from a target stationary distribution pdf using the Metropolis-Hastings algorithm. sym is a logical value that indicates whether the proposal distribution is symmetric. The default value is false, which corresponds to the asymmetric proposal distribution. If sym is true, for example, the proposal distribution is symmetric, proppdf and logproppdf are optional.

smpl = mhsample(...,'burnin',K) generates a Markov chain with values between the starting point and the k^th point omitted in the generated sequence. Values beyond the k^th point are kept. k is a nonnegative integer with default value of 0.

smpl = mhsample(...,'thin',m) generates a Markov chain with m-1 out of m values omitted in the generated sequence. m is a positive integer with default value of 1.

smpl = mhsample(...,'nchain',n) generates n Markov chains using the Metropolis-Hastings algorithm. n is a positive integer with a default value of 1. smpl is a matrix containing the samples. The last dimension contains the indices for individual chains.

[smpl,accept] = mhsample(...) also returns accept, the acceptance rate of the proposed distribution. accept is a scalar if a single chain is generated and is a vector if multiple chains are generated.

Examples

Estimate Moments Using Independent Metropolis-Hastings Sampling

Use Independent Metropolis-Hastings sampling to estimate the second order moment of a Gamma distribution.

rng default;  % For reproducibility
alpha = 2.43;
beta = 1;
pdf = @(x)gampdf(x,alpha,beta); % Target distribution
proppdf = @(x,y)gampdf(x,floor(alpha),floor(alpha)/alpha);
proprnd = @(x)sum(...
              exprnd(floor(alpha)/alpha,floor(alpha),1));
nsamples = 5000;
smpl = mhsample(1,nsamples,'pdf',pdf,'proprnd',proprnd,...
                'proppdf',proppdf);

Plot the results.

xxhat = cumsum(smpl.^2)./(1:nsamples)';
figure;
plot(1:nsamples,xxhat)

Random Walk Metropolis-Hastings Sampling

Use Random Walk Metropolis-Hastings sampling to generate sample data from a standard normal distribution.

rng default  % For reproducibility
delta = .5;
pdf = @(x) normpdf(x);
proppdf = @(x,y) unifpdf(y-x,-delta,delta);
proprnd = @(x) x + rand*2*delta - delta;   
nsamples = 15000;
x = mhsample(1,nsamples,'pdf',pdf,'proprnd',proprnd,'symmetric',1);

Plot the sample data.

figure;
h = histfit(x,50);
h(1).FaceColor = [.8 .8 1];

Version History

Introduced in R2006a