El muestreador MCMC no converge para una función concreta

Question

He escrito algún código MCMC que creía que funcionaba, pero para funciones más complicadas se rompe. El algoritmo MCMC que estoy utilizando, utiliza un simple algoritmo de Metrópolis.

En el código que adjuntaré a continuación, cuando uso:

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

Todo converge (es decir, mis paseos aleatorios convergen a una media). Esto se muestra en mi imagen de abajo:

Sin embargo, cuando uso la función más complicada en su lugar:

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

mis paseos al azar no van a ninguna parte. Para aclarar estos son los diagramas que estoy recibiendo:

Este es mi código de MATLAB, que he intentado hacer tan fácil de seguir como he podido. Hazme saber si algo no tiene sentido.

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

Algo que he notado es que cuando uso la f(x1,x2) más complicada mi algoritmo MCMC acepta básicamente todo (mi alfa es casi siempre la unidad). Sin embargo, con mi más simple f(x1,x2) = [1,2] el alfa hace muy (por lo que algunos casos son aceptados, algunos otros son rechazados) - que tiene sentido para mí.

¡¡¡Gracias por su ayuda!!!

P.D. Puedes copiar y pegar mi código directamente en MATLAB, es perfectamente ejecutable tal cual.

EDITAR/ACTUALIZAR: El mismo comportamiento ocurre incluso sin la función anterior, 'p1(x1,x2)'. Así que si usted acaba de dejar p = @(x1,x2) p2(x1,x2) Todavía tengo un problema de no convergencia, por lo que fundamentalmente p2(x1,x2) está causando problemas, y no estoy seguro de por qué.

Answer 1

Parece que tienes un problema de desbordamiento.

La expresión exp( -1/2*(f(x1,x2) - [x1,x2])*inv(cov_ML)*(f(x1,x2) - [x1,x2])' ) se evaluará a cero como la cantidad dentro de exp() tiende a -Inf . Cuando se divide por cero Matlab devuelve un NaN y min([1,NaN]) devolverá un 1 . Por eso alpha toma sistemáticamente un valor de 1 .

¿Es necesario exponer esa expresión? Una posible solución es dejar la expresión interior tal como está, sin el exp() . A menos que esto rompa lo que estás tratando de hacer (tbh, no entiendo muy bien todo lo que tienes en tu función de probabilidad).

El otro problema potencial es que su función f = @(x1,x2) x1.^2 + x2.^2 + 20 sólo devuelve un valor. ¿Quiere decir que f = @(x1,x2) [x1.^2, x2.^2] + 20 ?