I've written some MCMC code which I thought worked, but for more complicated functions it breaks down. The MCMC algorithm I am using, uses a simple Metropolis algorithm.

In the code which I will attach below, when I use:

f = @(x1,x2) [1,2];   % A simple function which only spits out [1,2] 

Everything converges (i.e. my random walks converge to a mean). This is shown in my image below:

When the simple function is used

However when I use the more complicated function instead:

f = @(x1,x2) x1.^2 + x2.^2 + 20;  % A nonlinearity (when this is used MCMC can't converge)

my random walks go nowhere. To clarify these are the diagrams I am getting: When the complicated function is used.

This is my MATLAB code which I tried to make as easy to follow as I could. Let me know if anything doesn't make sense.

clear all

%% DEFINE THE FIRST FUNCTION (kind of like a likelihood function)
N = 1;
sigma_ML = 0.03;   
cov_ML = eye(2)*sigma_ML;

f = @(x1,x2) x1.^2 + x2.^2 + 20;  % A nonlinearity (when this is used MCMC can't converge)
% f = @(x1,x2) [1,2];          % A simple function which only spits out [1,2] 
% My f(x1,x2) is used in p2 below:
p2 = @(x1,x2) 1./(2*pi*det(cov_ML))^(N/2) * ...
              exp( -1/2*(f(x1,x2) - [x1,x2])*inv(cov_ML)*(f(x1,x2) - [x1,x2])' );

%% DEFINE ANOTHER FUNCTION (basically like a prior function)
sigma_a = 1;
sigma_b = 1;
mu_a = 10;
mu_b = -20;
p1 = @(x1,x2) (1/(sqrt(2*pi*sigma_a^2))*exp(-1/(2*sigma_a^2)*(x1-mu_a).^2))...

p = @(x1,x2) p1(x1,x2).*p2(x1,x2);  % This is the function I will be using in my MCMC

nSamples = 500000;
t = 1;      % To keep track of how many total MCMC steps have been taken
idx = 2;    % To keep track of how many successful MCMC steps have been taken
x(1,:) = randn(1,2)+10;    % To start the algorithm

while t < nSamples
    t = t + 1;

    % SAMPLE FROM PROPOSAL (2D multivariate normal)
    xStar = mvnrnd(x(idx-1,:),eye(2));

    alpha(t) = min([1, p(xStar(1),xStar(2))/p(x(idx-1,1),x(idx-1,2))]);

    u = rand;
    if u < alpha(t)
        x(idx,:) = xStar;
        idx = idx + 1;
        x(idx,:) = x(idx-1,:);

        fprintf('%d / %d\n',t,nSamples);

Something I have noticed is that when I use the more complicated f(x1,x2) my MCMC algorithm accepts basically everything (my alpha is almost always unity). However, with my simpler f(x1,x2) = [1,2] the alpha does very (so some cases are accepted, some other ones are rejected) - which makes sense to me.

Thanks for your help!!!

P.S. You can copy-paste my code directly into MATLAB it is perfectly runnable as is.

EDIT/UPDATE: The same behviour happens even without the prior, 'p1(x1,x2)' function. So if you just let p = @(x1,x2) p2(x1,x2) I still get a non-convergence issue, so fundamentally p2(x1,x2) is causing issues, and I'm not sure why.


It seems like you've got an underflow issue.

The expression exp( -1/2*(f(x1,x2) - [x1,x2])*inv(cov_ML)*(f(x1,x2) - [x1,x2])' ) will evaluate to zero as the quantity inside exp() tends towards -Inf. When you divide by zero Matlab returns a NaN and min([1,NaN]) will return a 1. That's why alpha consistently takes a value of 1.

Do you need to exponentiate that expression? One potential fix is just to leave the interior expression as is, without the exp(). Unless this breaks what you're trying to do (tbh, I don't really understand what all you have going on in your likelihood function).

The other potential issue is that your function f = @(x1,x2) x1.^2 + x2.^2 + 20 only returns one value. Did you mean to have f = @(x1,x2) [x1.^2, x2.^2] + 20?

| cite | improve this answer | |
  • $\begingroup$ Yup my Likelihood function started out as a Likelihood function but I played around with it a bit so it's not 100% one any more. But I disagree with the f = @(x1,x2) x1.^2 + x2.^2 + 20 expression only producing one value. When I just put that into the matlab command window, this function doesn't only return one value. i.e. f(1,1) returns something different to f(5,8). Is this what you mean? $\endgroup$ – tisPrimeTime Sep 21 '16 at 5:24
  • $\begingroup$ Sorry, I meant that it returns a scalar, rather than an array. Your simple function that works returns a two-element array. Maybe this doesn't matter, though. $\endgroup$ – khonegger Sep 21 '16 at 5:33
  • $\begingroup$ Oh that is a good pick up. Lucky for me Matlab handles matrices in a way which didn't really effect which option I chose. $\endgroup$ – tisPrimeTime Sep 21 '16 at 5:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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