Calculando a razão de hastings, g (x | x ') / g (x' | x) para distribuições de propostas assimétricas no algoritmo MH?

Aug 18 2020

Eu entendo o algoritmo Metropolis. Onde eu fico confuso é o algoritmo MH, onde as distribuições de propostas assimétricas podem ser usadas.

Eu entendo que P (x) e P (x ') representam a probabilidade / densidade de probabilidade de x e x' de acordo com a distribuição alvo. Da mesma forma, eu entendo que g (x | x ') / g (x' | x) é um termo usado para corrigir uma distribuição de proposta assimétrica. Não estou confuso com seu propósito; Não entendo sua execução.

Como um problema de brinquedo, desenvolvi um amostrador de distribuição exponencial. Existem duas variações, uma que usa uma distribuição de proposta simétrica, a distância uniforme. E aquele que não: Ou seja Beta(a=3,b=2) - 0.5,. Escolhi esta distribuição porque (A) é assimétrica e principalmente positiva (mas ocasionalmente negativa, devido ao termo -0,5).

Não tenho ideia de como encontrar 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 distribuição da proposta assimétrica:

Gráfico de distribuição de proposta simétrica (com linhas comentadas invertidas):

Para responder a esta pergunta, edite o código para refletir uma g(x|x')/g(x'|x)distribuição válida para a proposta da Beta(a=3,b=2) - 0.5qual as perturbações são extraídas.

Respostas

2 jbuddy_13 Aug 18 2020 at 02:19

Dê uma olhada no código atualizado e gráfico abaixo. Observe que g (x | x ') / g (x' | x) é essencialmente uma medida de quão provável essas perturbações são vistas sob a distribuição proposta, que foi definida como Beta(a=3,b=2) -0.5.

Primeiro, encontre a diferença entre os eventos atuais e propostos. Em segundo lugar, ajuste para -0,5; chamaremos essas perturbações imparciais (onde -0,5 é uma tendência). Em terceiro lugar, encontre a probabilidade de cada perturbação (curr-> prop & prop-> curr). Por último, retorne a proporção como correction.

Usaremos essa correção e multiplicaremos por outros termos na acceptancedefinição da variável. É basicamente isso!

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