¿Calcular la relación de hastings, g (x | x ') / g (x' | x) para distribuciones de propuesta asimétricas en el algoritmo MH?
Entiendo el algoritmo de Metropolis. Donde me confundo es el algoritmo MH donde se pueden usar distribuciones de propuestas asimétricas.

Entiendo que P (x) y P (x ') representan la probabilidad / densidad de probabilidad de x y x' de acuerdo con la distribución objetivo. Asimismo, entiendo que g (x | x ') / g (x' | x) es un término utilizado para corregir una distribución de propuesta asimétrica. No estoy confundido por su propósito; No entiendo su ejecución.
Como problema de juguete, he desarrollado un muestreador de distribución exponencial. Hay dos variaciones, una que usa una distribución de propuesta simétrica, la dist uniforme. Y uno que no lo hace, a saber: Beta(a=3,b=2) - 0.5
. Elegí esta distribución porque (A) es asimétrica y en su mayoría positiva (pero ocasionalmente negativa, debido al término -0,5).
No tengo ni idea de cómo encontrarlo g(x|x')/g(x'|x)
.
Código:
def target(x,lam):
return int(x>0) * lam * np.exp(-x * lam)
def exponential_MH(hops,lam=3):
states = []
burn_in = int(hops*0.2)
current = lam
for i in range(hops):
states.append(current)
# movement = current + random.uniform(-1,1) # does not require asymmetric correction
movement = current + np.random.beta(a=3,b=2)-0.5 # requires asymmetric correction
acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)
event = random.uniform(0,1)
if acceptance > event:
current = movement
return states[burn_in:]
lam = 1
exp_samples = exponential_MH(hops=10_000,lam=lam)
lines = np.linspace(0,5,10_000)
exp_curve = [lam*np.exp(-l*lam) for l in lines]
plt.hist(exp_samples,normed=1,bins=20)
plt.plot(lines,exp_curve)
Gráfico de distribución de propuesta asimétrica:

Gráfico de distribución de propuesta simétrica (con líneas comentadas invertidas):

Para responder a esta pregunta, edite el código para reflejar una g(x|x')/g(x'|x)
distribución válida para la propuesta de la Beta(a=3,b=2) - 0.5
que se derivan las perturbaciones.
Respuestas
Eche un vistazo al código actualizado y grafique a continuación. Tenga en cuenta que g (x | x ') / g (x' | x) es esencialmente una medida de la probabilidad de que estas perturbaciones se vean bajo la distribución propuesta, que se ha definido como Beta(a=3,b=2) -0.5
.
Primero, encuentre la diferencia entre los eventos actuales y propuestos. En segundo lugar, ajuste el -0,5; llamaremos a estas perturbaciones insesgadas (donde -0.5 es un sesgo). En tercer lugar, encuentre la probabilidad de cada perturbación (curr-> prop & prop-> curr). Por último, devuelva la relación como correction
.
Usaremos esta corrección y la multiplicaremos por otros términos en la acceptance
definición de variable. ¡Eso es practicamente todo!
def target(x,lam):
return int(x>0) * lam * np.exp(-x * lam)
def correct(prop,curr,a=3,b=2):
x0 = curr - prop + 0.5
x1 = prop - curr + 0.5
b0 = beta.pdf(x=x0, a=a, b=b)
b1 = beta.pdf(x=x1, a=a, b=b)
return b0/b1
def exponential_MH(hops,lam=3):
states = []
burn_in = int(hops*0.2)
current = 1
for i in range(hops):
states.append(current)
movement = current + np.random.beta(a=3,b=2)-0.5 # requires assymetric correction
correction = correct(curr=current,prop=movement)
acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)*correction
event = random.uniform(0,1)
if acceptance > event:
current = movement
return states[burn_in:]
