MCMC Sampler not Converging for Particular Function

Question

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:

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:

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
clc

%% 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))...
        .*(1/sqrt((2*pi*sigma_b^2))*exp(-1/(2*sigma_b^2)*(x2-mu_b).^2));    



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



%% INITIALISE VARIABLES
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




%% RUN MCMC SAMPLER
while t < nSamples
    t = t + 1;

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


    % CALCULATE THE M-H ACCEPTANCE PROBABILITY
    alpha(t) = min([1, p(xStar(1),xStar(2))/p(x(idx-1,1),x(idx-1,2))]);

    % ACCEPT OR REJECT?
    u = rand;
    if u < alpha(t)
        x(idx,:) = xStar;
        idx = idx + 1;
    else
        x(idx,:) = x(idx-1,:);
    end

    if(mod(t,10000)==0)
        fprintf('%d / %d\n',t,nSamples);
    end
end

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.

Answer 1

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?