MCMC.m

% Created by: Dima Damen
% Date: December 2008

% This function takes a probability density function, and draws a set of
% samples of size n and plots the histogram of those samples along with the
% original pdf showing convergence.

% examples of some pdf you can experiment with are:
% pdf = @(x) 0.3*normpdf (x,0.3,0.05)+0.7*normpdf(x,0.7,0.2);
% pdf = @(x) normpdf (x, 4, 3);
% pdf = @(x) 0.1*normpdf (x,0.3,0.1)+0.3*normpdf(x,0.6,0.3)+0.6*normpdf(x,0.8,0.2);
% pdf = @(x) exppdf (x, 2);
% pdf = @(x) 0.4*exppdf(x,2) + 0.6*normpdf (x, 4, 2);
% pdf = @(x) betapdf (x, 2, 1);
% pdf = @(x) gampdf (x,1.5,2);

function MCMC (pdf, nsamples)

    %proppdf is the distribution from which you can sample/propose (Q)
    delta = 1;
    proppdf = @(x,y) unifpdf (y-x, -delta, 1);
    proprnd = @(x) x + rand*2*delta - delta;

    % metropolis-hastings algorithm 
    x = mhsample (0.4, nsamples, 'pdf', pdf, 'proppdf', proppdf, 'proprnd', proprnd);

    % d is the domain of the function you want to run
    d = min(x):0.01:max(x);

    % this draws the actual (hidden) distribution
    d_x = pdf (d);

    n = hist (x,50);
    prod = max(n)/max(d_x);
    yVal = d_x * prod;
    figure; hist(x,50);
    hold on;
    plot (d, yVal, 'r-', 'linewidth', 3);
end