3
$\begingroup$

I am learning Markov Chain Monte Carlo (MCMC) simulation as of the moment. My background is civil engineering and please excuse my ignorance if some of the questions are quite obvious.

I want to learn MCMC and implement it into my matlab code. To test my understanding, I want to perform an MCMC with an analytical solution based from conjugate pair.

$$f(p|x)=\frac{f(p)f(x|N,p)}{\int_0^1{f(p)f(x|N,p)}dr}$$

Basically, I want to infer the probability $p\in [0,1]$ based from new data $x$ given $N$ trials. I used the beta-binomial conjugate pair as a test example

prior: $f(p) = Beta(1,1)$ likelihood: $f(x|N,p) = Bin(x,N,p)$ Posterior: $f(p|x) = Beta(x+1,N-x+1)$

I provided my code below.

This is how I implemented my metropolis algorithm of MCMC:

  1. Initialize a value for the quantity of interest/parameter $p_0$
  2. I draw a candidate $\xi$ from a proposal distribution of uniform distribution $U(0,1)$
  3. Calculate acceptance rate $\alpha = \min\{1, \frac{\pi(\xi)}{\pi(p_k)}\}$.
  4. Accept candidate if random $r\leq\alpha$ such that $p_{k+1}=\xi$ else reject the sample such that $p_{k+1}=p_k$
  5. Repeat steps 2 to 4 for $N$ iterations

Take note my target distribution is $\pi(p) = Bin(148,150,p)\times Beta(1,1)$

My question is: a) How to choose a correct proposal distribution? In my understanding, the proposal distribution should be dependent from the candidate value $f^{*}(p_k|\xi)$ but the one I used is only uniform distribution from 0 to 1 which does not depend from the previous of value of $p$.

I have seen other codes/scripts that draws candidate from a normal distribution with mean value dependent from the previous value $Normal(p_{k+1},\sigma)$. But this does not work for my case because it will draw samples greater than 1.

How do I choose a correct proposal distrbution?

clc, clear; format compact; format shortG;

%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% MCMC
%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 100000;

seed = 0;
rng(seed);
samples = [];
x0 = 0.5;
sigma = 4;

% Metropolis algorithm
samples = [samples, x0];
for k = 1:N
    k
    % generate candidate value from proposal function of U(0,1)
    candidate = rand;

    % calculate probability of accepting the candidate
    alpha = min(1, target_prob(candidate)/target_prob(samples(end)) );

    % determine whether to accept or reject
    if rand < alpha % accept
        samples = [samples, candidate];
    else
        samples = [samples, samples(end)];
    end
 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Bayesian inference using Beta-Binomial pair
%%%%%%%%%%%%%%%%%%%%%%%%%%%
exp_x = 147.8716;
Nsamples = 150;

alpha = exp_x + 1;
beta = Nsamples - exp_x + 1;

pd = makedist('Beta','A',alpha,'B',beta );

%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = 0.85:0.0001:1;

f2 = figure;
set(f2, 'units','inches','position',[1,1,6,4]);

y = pdf(pd, x );
h = histogram(samples,200,'Normalization','pdf','EdgeColor','none'); hold on;
p = plot(x,y,'LineWidth',1.5,Color='#abdda4',LineStyle='-.', ...
    LineWidth=2); hold on

legend([h,p],{'MCMC samples','Analytical'},location ="northwest")
legend('boxoff')
xlim([0.85,1]); ylim([0,100]);

xlabel('Parameter '); ylabel('Probabilty Density');

function [posterior] = target_prob(r)
%%% function to evaluate the target distribution
    Nsamples = 150;
    x = round(147.8716);

    prior = betapdf(r,1,1);
    likelihood = binopdf(x, Nsamples, r);
    
    posterior = prior * likelihood;
end

enter image description here

$\endgroup$

1 Answer 1

5
$\begingroup$

You've specified a class of MCMC algorithms referred to as Metropolis-Hastings.

With these algorithms, choosing the proposal distribution is the key question when specializing the algorithm. Under very mild conditions for selecting the proposal distribution, the algorithm will cover posterior properly and so your results should be asymptotically correct, as long as you don't make a big mistake like choosing a proposal distribution that doesn't properly cover all possible parameter values.

But as they say, asymptotically we're all dead. You can pick a proposal distribution that asymptotically will properly cover the posterior yet even after running your algorithm for 10 million steps does not have a good representation of the posterior due to not moving quickly through the posterior. So making an efficient algorithm mean making one that covers the posterior quickly. This is a bit tricky.

TLDR is that the closer the proposal distribution is to the target, the more efficient we will be in standard Metropolis Hastings algorithms. But note that this is a chicken and egg problem: if we knew the shape of the posterior then we'd already have our answer!

As such, there can be a bit of trial and error in our problem. For a very simple MH algorithm, we fiddle with the $\sigma$ parameter until we move relatively quickly across the posterior. More advanced algorithms will actually learn the $\sigma$ parameter during the MCMC run. Even more advanced will take advantage of the derivatives of the posterior to move more quickly without proposing values that get rejected too often.

If you're just trying to learn about MH algorithms, I'd suggest playing with a few different $\sigma$'s and get a feeling for how they work (too small: moves too slowly. Too big: rejects proposal too often). If you really want to make an efficient MCMC algorithm, I'd suggest using implementations of the more advanced algorithms.

Also one note: you mention that your proposal is not good since it can propose impossible values (i.e. $p > 1$). Theoretically, this is not an issue: impossible values will automatically be rejected and you will still have valid inference. However, it may not be computationally efficient since you are wasting so many samples on impossible outcomes, so you might reparameterize your problem ($\theta = \log(p/(1-p))$) if that reduces your rejected sample rate.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.