¿Calcular la relación de hastings, g (x | x ') / g (x' | x) para distribuciones de propuesta asimétricas en el algoritmo MH?

Aug 18 2020

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.5que se derivan las perturbaciones.

Respuestas

2 jbuddy_13 Aug 18 2020 at 02:19

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 acceptancedefinició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:]