¿Por qué este algoritmo MCMC para estimar los parámetros de una ecuación lineal no converge a la distribución posterior?

Question

Como una especie de prueba de principio estoy tratando de estimar los parámetros de una ecuación lineal (antes de pasar a las EDOs) utilizando el muestreo de Monte Carlo de cadenas de Markov. El post que estoy siguiendo se puede encontrar en este enlace https://sciencehouse.wordpress.com/2010/06/23/mcmc-and-fitting-models-to-data/

Esencialmente, lo que estoy haciendo no funciona y me pregunto si alguien puede detectar el motivo.

El modelo $$y=mx+c$$ evaluado en $m=2$ y $c=4$ como un modelo lineal simple:

El algoritmo

1. Inicializar $m0$ y $c0$ para que sean valores aleatorios (elegí arbitrariamente 25 para ambos) y evalúo una función que mide la discrepancia entre los datos experimentales y los simulados (donde, por supuesto, en este caso los datos experimentales también se simulan con los parámetros correctos). Utilizo la suma de cuadrados $SSR$ :

$$ SSR=\frac{(SimulatedData-ExperimentalData)^2}{2*sigma^2}$$

Dado que no hay error experimental en este modelo, se eligió arbitrariamente que sigma fuera 0,0001 en todos los casos.

2. A continuación, comienzo un bucle y muestreo de anteriores $_m$ y $_c$ , produciendo $m_prop$ y $c_prop$ respectivamente, que se distribuyen uniformemente entre 0 y 100.
3. Calcule el $SSR$ con los parámetros muestreados
4. Calcular el cociente de probabilidades $LR$ : $$ LR=exp[\frac{(-ln(SSR proposed)+ln(SSR current)}{2*sigma^2}]$$
5. Aceptar aleatoriamente los nuevos parámetros mediante un nuevo muestreo de la distribución uniforme entre 0 y 1 (lo llamamos $r$ ). Si $LR>r$ se aceptan los nuevos parámetros y $m_prop$ se convierte en $m0$ y $c_prop$ se convierte en $c0$

Como pregunta al margen, ¿por qué algunos dicen que la probabilidad de aceptación es $min(1,LR)$ mientras que el puesto que estoy siguiendo utiliza la aceptación aleatoria?

6. Iterar hasta la convergencia

Resultados

---------edit--------------

Por cortesía de la respuesta de Joriki, he eliminado el $2*sigma^2$ de la función de verosimilitud, de modo que el paso 4 queda así:

4. Calcular el cociente de probabilidades $LR$ : $$ LR=exp[(-ln(SSR proposed)+ln(SSR current)]$$

Esto ha mejorado considerablemente los resultados:

Código Matlab

Para completar y dar los detalles del algoritmo, aquí está el código que estoy usando para generar los gráficos:

interval=1:100;
y=mx_plus_c((2),interval,(4));
n=1e5; %iterations
results=zeros(n,8);
m0=25;
c0=25;
y_sim=mx_plus_c(m0,interval,c0);
SSR_cur=SSR(y,y_sim,0.0001);
mlb=log(1);   %lb lower bound, ub upper bound respoectively 
mub=log(100);
clb=log(1);
cub=log(100);
for i=1:n,
    m_prop   =   (mlb+(mub-mlb).*rand(1,1)   )   ;
    c_prop   =    (clb+(cub-clb).*rand(1,1)   )  ;
    SSR_prop=SSR(y,mx_plus_c(m_prop,interval,c_prop),0.0001);
    likelihood_ratio=exp((-log(SSR_prop)+log(SSR_cur)));
    r   =    ( (0+(1-0).*rand(1,1)   )  );
    if likelihood_ratio>r,
        m0=m_prop;
        c0=c_prop;
        SSR_cur=SSR_prop;
    end
    results(i,1)=m0;
    results(i,2)=c0;
    results(i,3)=SSR_cur;
    results(i,4)=m_prop;
    results(i,5)=c_prop;
    results(i,6)=SSR_prop;
    results(i,7)=likelihood_ratio;
    results(i,8)=r;

end
% results
figure(2)
hist(results(:,1),1000)
title('Distribution of m parameters','fontsize',16) 
xlabel('Parameter Values','fontsize',16)
ylabel('Frequency','fontsize',16)

figure(3)
hist(results(:,2),1000)
title('Distribution of c parameters','fontsize',16) 
xlabel('Parameter Values','fontsize',16)
ylabel('Frequency','fontsize',16)

Où : mx_plus_c es una función de Matlab:

function y=mx_plus_c(m,x,c)
y=m*[x]+c;
end

Y SSR también es una función de Matlab

function chi2 = SSR(sim,exp,sigma)
chi2 =sum(  (  ([sim]-[exp]).^2)./ ( 2*sigma^2) ) ;

Answer 1

Hay varios problemas; sin ver los detalles de la simulación, es difícil saber cuál de ellos, si es que hay alguno, está causando el fallo.

• Cada vez eliges los parámetros de nuevo, no aprovechas los avances anteriores para mejorar los parámetros. La idea es proponer nuevos parámetros en la vecindad de los parámetros actuales, no tomar muestras de la misma distribución cada vez. (Lo que es exactamente la "vecindad" depende de la aplicación específica, y se puede cambiar con el tiempo para dar grandes pasos al principio y luego pasos cada vez más pequeños).
• No creo que debas tomar el logaritmo y dividir por $2\sigma^2$ de nuevo en el cálculo de la relación de probabilidad - sólo debe exponer la diferencia entre los errores al cuadrado. Dado que su $\sigma$ es pequeño, esto puede estar causando problemas.
• No dices cómo determinas la convergencia. Con ese pequeño valor de $\sigma$ y esa amplia previa que nunca cambias, si terminas después de un cierto número de intentos fallidos de cambiar los parámetros, esto puede estar dejándote en un lugar aleatorio porque es poco probable que progreses en cada paso.

Por lo general, para depurar este tipo de cosas, hay que observar los valores a medida que evolucionan y tratar de comprobar en cada etapa si lo que está ocurriendo se corresponde con lo que se esperaba que ocurriera, en lugar de limitarse a hacer un gran gráfico al final que no muestre lo que ha ido mal.

En cuanto a su pregunta complementaria, "¿por qué algunos dicen que la probabilidad de aceptación es $\min(1,LR)$ mientras que el post que estoy siguiendo utiliza la aceptación aleatoria", no entiendo la pregunta - utilizando una probabilidad de aceptación de $\min(1,LR)$ es aceptación aleatoria (es decir, con una cierta probabilidad), y es exactamente lo que estás haciendo: Si $LR\gt1$ siempre aceptas, y si no, aceptas con probabilidad $LR$ .