Cálculo da probabilidade marginal de amostras MCMC

Esta é uma pergunta recorrente (veja este post , este post e este post ), mas eu tenho uma opinião diferente.

Suponha que eu tenha várias amostras de um amostrador genérico do MCMC. Para cada amostra , eu sei o valor da probabilidade do log e do log anterior . Se ajudar, também sei o valor da probabilidade do log por ponto de dados, (essas informações ajudam em certos métodos, como WAIC e PSIS-LOO).log f ( x | θ ) log f ( θ ) log f ( x i | θ $\theta$$\log f(\textbf{x} | \theta)$$\log f(\theta)$$\log f(x_i | \theta)$

Quero obter uma estimativa (bruta) da probabilidade marginal, apenas com as amostras que tenho e, possivelmente, algumas outras avaliações de função (mas sem executar novamente um MCMC ad hoc ).

Primeiro de tudo, vamos limpar a mesa. Todos sabemos que o estimador harmônico é o pior estimador de todos os tempos . Vamos continuar. Se você estiver fazendo uma amostragem de Gibbs com anteriores e posteriores na forma fechada, poderá usar o método de Chib ; mas não sei como generalizar fora desses casos. Também existem métodos que exigem que você modifique o procedimento de amostragem (como via posterior temperado ), mas não estou interessado nisso aqui.

A abordagem em que estou pensando consiste em aproximar a distribuição subjacente com uma forma paramétrica (ou não paramétrica) e depois descobrir a constante de normalização como um problema de otimização 1-D (ou seja, o que minimiza algum erro entre e , avaliadas sobre as amostras). No caso mais simples, suponha que o posterior seja aproximadamente multivariado normal, eu posso ajustar como um normal multivariado e obter algo semelhante a uma aproximação de Laplace (eu poderia querer usar algumas avaliações de funções adicionais para refinar a posição de o modo). No entanto, eu poderia usar comoZ Z Z g ( θ ) f ( x | θ ) f ( θ ) g ( θ ) g ( θ $g(\theta)$$Z$$Z$$Z g(\theta)$$f(\textbf{x}|\theta) f(\theta)$$g(\theta)$$g(\theta)$uma família mais flexível, como uma mistura variacional de distribuições multivariadas de $t$ .

Compreendo que esse método funcione apenas se $Z g(\theta)$ for uma aproximação razoável de $f(\textbf{x}|\theta) f(\theta)$ , mas qualquer motivo ou conto preventivo sobre por que seria muito imprudente faça? Alguma leitura que você recomendaria?

A abordagem totalmente não paramétrica usa alguma família não paramétrica, como um processo Gaussiano (GP), para aproximar (ou alguma outra transformação não linear do mesmo, como como raiz quadrada) e quadratura bayesiana para integrar implicitamente sobre o alvo subjacente (veja aqui e aqui ). Essa parece ser uma abordagem alternativa interessante, mas de espírito análogo (observe também que os GPs seriam difíceis de manejar no meu caso).$\log f(\textbf{x}|\theta) + \log f(\theta)$

Infelizmente, a extensão de Chib e Jeliazkov (2001) fica rapidamente cara ou altamente variável, razão pela qual não é muito usada fora dos casos de amostragem de Gibbs.

Embora existam muitas maneiras e abordagens para o problema de estimativa de constante normalização (como ilustrado pelas palestras bastante diversas no workshop Estimando Constante que realizamos na semana passada na Universidade de Warwick, slides disponíveis lá ), algumas soluções exploram diretamente a saída do MCMC .$\mathfrak{Z}$

1. Como você mencionou, o estimador de média harmônica de Newton e Raftery (1994) é quase sempre invariavelmente pobre por ter uma variação infinita. No entanto, existem maneiras de evitar a maldição de variação infinita usando um alvo de suporte finito na identidade média harmônica escolhendoαcomo o indicador de uma região HPD para a parte posterior. Isso garante uma variação finita, removendo as caudas na média harmônica. (Detalhes podem ser encontrados emum artigo que escrevi com Darren Wraithe em umcapítulo sobre normalização de constantesescritas com Jean-Michel Marin.) Em resumo, o método recicla a saída do MCMCθ1,…,θMidentificando oβ( 20% dizem que os maiores valores do alvoπ(θ)f(x|θ)e criando

   $$\int \dfrac{\alpha(\theta)}{\pi(\theta)f(x|\theta)}\text{d}\pi(\theta|x)=\frac{1}{\mathfrak{Z}}$$ $\alpha$$\theta_1,\ldots,\theta_M$$\beta$$\pi(\theta)f(x|\theta)$$\alpha$como um uniforme através da união das bolas centrado no aqueles maior densidade (HPD) simulações e com raio ρ , o que significa que a estimativa da constante de normalização Z é dada por Z - 1 = $\theta^0_i$$\rho$$\mathfrak{Z}$ sedé a dimensão deθ(as correções se aplicam a bolas que se cruzam) e seρé pequeno o suficiente para que as bolas nunca se cruzem (significando que, na melhor das hipóteses, apenas um indicador nas bolas é diferente de zero). A explicação para odenominadorαM2é que esta é uma soma dupla determosβM2: $$\hat{\mathfrak{Z}}^{-1}=\underbrace{\frac{1}{\beta M^2}\sum_{m=1}^M}_{\text{double sum over}\\\beta M\text{ ball centres }\theta_i^0\\\text{and $M$ simulations } \theta_m} \underbrace{\mathbb{I}_{(0,\rho)}(\min_i||\theta_m-\theta^0_i||)\{\pi(\theta_m)f(x|\theta_m)\}^{-1}\big/\overbrace{\pi^{d/2}\rho^d\Gamma(d/2+1)^{-1}}^{\text{volume of ball with radius $\rho$}}}_{\dfrac{\beta M\alpha(\theta_m)}{\pi(\theta_m)f(x|\theta_m)}}$$ $d$$\theta$$\rho$$\alpha M^2$$\beta M^2$ com cada termo emθmintegrado aZ-1. $$\frac{1}{\beta M}\sum_{i=1}^{\beta M} \underbrace{\frac{1}{M}\sum_{m=1}^M {\cal U}(\theta_i^0,\rho)(\theta_m)}_{\text{same as with $\min$}} \times \frac{1}{\pi(\theta_m)f(x|\theta_m)}$$ $\theta_m$${\mathfrak{Z}}^{-1}$

2. Outra abordagem é transformar a constante de normalização em um parâmetro. Isso soa como uma heresia estatística, mas o artigo de Guttmann e Hyvärinen (2012) me convenceu do contrário. Sem entrar muito em detalhes, a idéia pura é transformar a probabilidade logarítmica observada n ∑ i = 1 f ( x i | θ ) - n log ∫ exp f ( x | θ ) d x em uma probabilidade logarítmica conjunta n ∑ i = 1 [ $\mathfrak{Z}$

   $$\sum_{i=1}^n f(x_i|\theta) - n \log \int \exp f(x|\theta) \text{d}x$$ que é a probabilidade logarítmica de um processo de ponto de Poisson com função de intensidade exp { f ( x | θ ) + ν + log n $$\sum_{i=1}^n[f(x_i|\theta)+\nu]-n\int\exp[f(x|\theta)+\nu]\text{d}x$$ $$\exp\{ f(x|\theta) + \nu +\log n\}$$ Este é um modelo alternativo, pois a probabilidade original não aparece como marginal do que foi mencionado acima. Somente os modos coincidem, com o modo condicional em ν fornecendo a constante de normalização. Na prática, a probabilidade do processo de Poisson acima está indisponível e Guttmann e Hyvärinen (2012) oferecem uma aproximação por meio de uma regressão logística. Para se conectar ainda melhor com sua pergunta, a estimativa de Geyer é um MLE, portanto, solução para um problema de maximização.
3. $\pi(\theta|x)$$\pi(\theta|x)$$g(\theta)$$\pi(\theta|x)$$g(\theta)$) Com os regressores sendo os valores de ambas as densidades, normalizados ou não. Isso está diretamente relacionado à amostragem de ponte de Gelman e Meng (1997), que também recicla amostras de diferentes alvos. E versões posteriores, como o MLE de Meng.
4. Uma abordagem diferente que obriga a executar um amostrador MCMC específico é a amostragem aninhada de Skilling . Embora eu [e outros] tenhamos algumas reservas quanto à eficiência do método, ele é bastante popular em astrostatística e cosmologia, com softwares disponíveis como multinacionais .
5. $H_0: \theta=\theta_0$$\xi$$\pi_1(\theta)\pi_2(\xi)$$H_0$ $$\mathfrak{B}_{01}(x)=\dfrac{\pi^\theta(\theta_0|x)}{\pi_1(\theta_0)}$$ $\pi^\theta(\theta_0|x)$$\theta$$\theta_0$$H_0: \theta=\theta_0$ $$m_0(x)=\int_\Xi f(x|\theta_0,\xi)\pi_2(\xi)\text{d}\xi$$ $$m_a(x)=\int_{\Theta\times\Xi} f(x|\theta,\xi)\pi_1(\theta)\pi_2(\xi)\text{d}\theta\text{d}\xi$$

[Aqui está um conjunto de slides que escrevi sobre a estimativa de constantes de normalização para um workshop do NIPS em dezembro passado.]