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)
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.