MCMC para lidar com problemas de probabilidade plana

9

Eu tenho uma probabilidade bastante plana de levar o amostrador Metropolis-Hastings a se mover pelo espaço de parâmetros de maneira muito irregular, ou seja, nenhuma convergência pode ser alcançada, independentemente dos parâmetros de distribuição da proposta (no meu caso, é gaussiano). Não há alta complexidade no meu modelo - apenas 2 parâmetros, mas parece que o MH não pode lidar com essa tarefa. Então, existe algum truque para solucionar esse problema? Existe um amostrador que não produziria cadeias de Markov se movendo muito longe para as caudas posteriores?

Atualização do problema:
tentarei reformular minha pergunta dando mais detalhes. Antes de tudo, descreverei o modelo.
Eu tenho um modelo gráfico com dois nós. Cada nó é governado por um modelo de auto-Poisson (Besag, 1974) da seguinte maneira: Ou, uma vez que existem apenas dois nós e assumindo intensidades globais iguais :

p(Xj|Xk=xk,kj,Θ)Poisson(eθj+jkθkjxk)
p ( X 2 | X 1 = x 1 , θ , α ) ~ P o i s s o n ( e θ + α x 1 )
p(X1|X2=x2,θ,α)Poisson(eθ+αx2)
p(X2|X1=x1,θ,α)Poisson(eθ+αx1)

Como é um campo de Markov, a distribuição conjunta (ou probabilidade de realização ) é a seguinte: Como eu assumi as anteriores planas para e , a posterior é então proporcional a DesdeX=[x1,x2]

p(X)=exp(θ(x1+x2)+2x1x2α)Z(θ,α)=exp(E(θ,α,X))Z(θ,α)
αθ
π(θ,α|X)exp(E(θ,α,X))Z(θ,α)
Z(θ,α)em geral, é muito difícil de avaliar (muitos somatórios). Estou usando o método de variável auxiliar devido a J. Moller (2006). De acordo com esse método, primeiro desenhei uma amostra de dados do amostrador Gibbs (como condicionais são apenas distribuições de poisson), depois desenhei uma proposta da distribuição gaussiana e calculei adequadamente os critérios de aceitação . E aqui eu recebo uma cadeia selvagem de Markov. Quando imponho alguns limites dentro dos quais a cadeia pode se mover, o amostrador parece convergir para alguma distribuição, mas uma vez que movo pelo menos um limite, a distribuição resultante também se move e sempre mostra transe. Eu acho que @ Xi'an é errado - o posterior pode ser impróprio.XH(X,α,θ|X,α,θ)
Tomas
fonte
1
Uma possibilidade é usar um parâmetro de escala maior para obter etapas maiores. Você pode estar interessado no pacote R mcmce no comando metroptambém. Você provavelmente precisará de um amostrador adaptável. Este amostrador (o twalk) pode ser usado nesse tipo de caso, uma vez que é adaptável (talvez apenas como uma "segunda opinião"). É implementado em R, C e Python. Os códigos podem ser baixados de uma das páginas do autor .
@Procrastinator Você pode elaborar mais sobre o que você quer dizer com "parâmetro de escala maior"? Significa usar parâmetros de variação maiores para propostas?
Tomas
2
Apenas deixe-me esclarecer primeiro que, se a probabilidade for baixa, você realmente não quer que o seu amostrador "não se mova muito para a cauda posterior". O que se deseja é amostrar adequadamente a partir da distribuição (ambas, caudas e centro). Ao usar um algoritmo MH com propostas gaussianas, você precisa escolher parâmetros de escala / matriz de covariância que determinam o comprimento das etapas. Estes devem ser escolhidos para: 1. Amostrar adequadamente da distribuição e 2. Obter uma taxa de aceitação razoável.
se você só tem dois parâmetros, em seguida, integração numérica é provavelmente uma alternativa melhor
probabilityislogic
há algo errado com a expressão de probabilidade conjunta. Se você tentar somar obtém . portanto, a probabilidade é imprópria conforme está escrito atualmente. p ( x 2 | α θ ) = g ( x 2 ) Σ x 1 = 0 exp ( x 1 [ θ + 2 α x 2 ] ) = x1p(x2|αθ)=g(x2)x1=0exp(x1[θ+2αx2])=
probabilityislogic

Respostas:

8

Acho surpreendente que uma probabilidade plana produza problemas de convergência: geralmente é o caso oposto que causa problemas! A primeira verificação usual para essas situações é garantir que o posterior seja adequado : caso contrário, isso explicaria excursões intermináveis ​​nas "caudas". Se a posterior for realmente adequada, você pode usar propostas de cauda mais gorda, como uma distribuição de Cauchy ... E um algoritmo adaptável à Roberts e Rosenthal.

Se isso ainda "não funcionar", sugiro considerar uma reparameterização do modelo, usando, por exemplo (ou seja, se não houver outra parametrização natural), uma transformação logística, (com um possível parâmetro de escala), que traz o parâmetro para o quadrado da unidade.

φ(x)=exp(x)/{1+exp(x)}

Em relação às respostas anteriores, a amostragem de Gibbs parece uma solução mais provável do que a aceitação-rejeição, que exige encontrar um limite e escalar a distribuição t para a posterior, o que não parecia viável para o amostrador mais robusto de Metropolis-Hastings ...

Xi'an
fonte
@Xian, obrigado pelo feedback sobre o voto negativo. Existe realmente alguma situação em que você preferiria aceitar-rejeitar sobre MH?
gui11aume
@ gui11aume: se você pode produzir um algoritmo de aceitação / rejeição com um limite suficientemente pequeno para garantir uma taxa de aceitação razoável, a aceitação / rejeição é, sem dúvida, preferível a Metropolis-Hastings. No entanto, é improvável que isso aconteça com (a) grandes dimensões e / ou (b) complexos, possivelmente multimodais, metas ...
Xi'an
2

Você pode anotar a distribuição do seu primeiro parâmetro condicional no seu segundo parâmetro e vice-versa? Nesse caso, a amostragem de Gibbs seria uma opção viável. São apenas algumas linhas de código e podem se misturar quase instantaneamente em muitos casos.

David J. Harris
fonte
1

EDIT: Veja a resposta de @ Xi'an e a discussão depois para ver os problemas com a seguinte abordagem.

Se Metropolis-Hastings falhar e seu modelo for relativamente simples, você pode pensar em usar o algoritmo de aceitação / rejeição com a distribuição de Student com um baixo grau de liberdade (1-6) para as propostas.t

Se você usar R, você pode facilmente simular um Student com . Se você não tem uma maneira fácil de gerar variáveis ​​com seu software, mas pode simular a , então desenhar a variação de um Gaussiano de em cada etapa e simular um Gaussiano com essa variação é equivalente.t Γ Γtrt()tΓΓ

gui11aume
fonte