simulando amostras aleatórias com um determinado MLE

Essa pergunta de validação cruzada que pergunta sobre simular uma amostra condicional a uma soma fixa me lembrava de um problema que George Casella me colocou .

$f(x|\theta)$$(X_1,\ldots,X_n)$$\theta$θ(X1,...,Xn) θ (X1,...,Xn

$$\hat{\theta}(x_1,\ldots,x_n)=\arg\min \sum_{i=1}^n \log f(x_i|\theta)$$ $\theta$$(X_1,\ldots,X_n)$$\hat{\theta}(X_1,\ldots,X_n)$

Por exemplo, considere uma distribuição , com o parâmetro de localização , cuja densidade é Se como podemos simular condicional em ? Neste exemplo , a distribuição de não possui uma expressão de formulário fechado. μf(x | μ)= Γ ( 3 $\mathfrak{T}_5$$\mu$

$$f(x|\mu)=\dfrac{\Gamma(3)}{\Gamma(1/2)\Gamma(5/2)}\,\left[1+(x-\mu)^2/5\right]^{-3}$$ $$(X_1,\ldots,X_n)\stackrel{\text{iid}}{\sim} f(x|\mu)$$ $(X_1,\ldots,X_n)$$\hat{\mu}(X_1,\ldots,X_n)=\mu_0$$\mathfrak{T}_5$$\hat{\mu}(X_1,\ldots,X_n)$

Uma opção seria usar uma variante HMC restrita, conforme descrito em Uma família de métodos MCMC em coletores definidos implicitamente por Brubaker et al (1). Isso requer que possamos expressar a condição de que a estimativa de probabilidade máxima do parâmetro de localização é igual a alguns fixos, como alguma restrição holonômica implicitamente definida (e diferenciável) c ( { x i } N i = 1 ) = 0 . Podemos então simular uma dinâmica hamiltoniana restrita sujeita a essa restrição e aceitar / rejeitar dentro de uma etapa de Metropolis-Hastings, como no HMC padrão.$\mu_0$$c\left(\lbrace x_i \rbrace_{i=1}^N\right) = 0$

A probabilidade logarítmica negativa é que possui derivadas parciais de primeira e segunda ordem em relação ao parâmetro de localizaçãoμ∂

$$\mathcal{L} = -\sum_{i=1}^N \left[ \log f(x_i \,|\, \mu) \right] = 3 \sum_{i=1}^N \left[ \log\left(1 + \frac{(x_i - \mu)^2}{5}\right)\right] + \text{constant}$$ $\mu$ Uma estimativa de probabilidade máxima deμ0é então implicitamente definida como uma solução para c=N∑i=1[2(μ0-xi $$\frac{\partial \mathcal{L}}{\partial \mu} = 3 \sum_{i=1}^N \left[ \frac{2(\mu - x_i)}{5 + (\mu - x_i)^2}\right] \quad\text{and}\quad \frac{\partial^2 \mathcal{L}}{\partial \mu^2} = 6 \sum_{i=1}^N \left[\frac{5 - (\mu - x_i)^2}{\left(5 + (\mu - x_i)^2\right)^2}\right].$$ $\mu_0$ $$c = \sum_{i=1}^N \left[ \frac{2(\mu_0 - x_i)}{5 + (\mu_0 - x_i)^2}\right] = 0 \quad\text{subject to}\quad \sum_{i=1}^N \left[\frac{5 - (\mu_0 - x_i)^2}{\left(5 + (\mu_0 - x_i)^2\right)^2}\right] > 0.$$

Não tenho certeza se existem resultados sugerindo que haverá um MLE exclusivo para para determinado { x i } N i = 1 - a densidade não é côncava em µ em μ, portanto, não parece trivial garantir isso. Se houver uma única solução única, o acima definido implicitamente define um coletor dimensional N - 1 conectado embutido em R N correspondente ao conjunto de { x i } N i = 1 com MLE para μ igual a μ $\mu$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$N - 1$$\mathbb{R}^N$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$\mu_0$. Se houver várias soluções, o coletor pode consistir em vários componentes não conectados, alguns dos quais podem corresponder a mínimos na função de probabilidade. Nesse caso, precisaríamos ter algum mecanismo adicional para mover-se entre os componentes não conectados (como a dinâmica simulada geralmente permanecerá confinada a um único componente) e verificar a condição de segunda ordem e rejeitar uma movimentação, se corresponder à mudança para um mínimo na probabilidade.

Se usarmos para denotar o vetor [ x 1 … x N ] T e introduzir um estado de momento conjugado p com matriz de massa M e um multiplicador de Lagrange λ para a restrição escalar c ( x ), então a solução para o sistema de EDOs d$\boldsymbol{x}$$\left[ x_1 \dots x_N\right]^{\rm T}$$\boldsymbol{p}$$\mathbf{M}$$\lambda$$c(\boldsymbol{x})$ dada condição inicialx(0)=x0,P(0)=p0comC(x0)=0e ∂

$$\frac{{\rm d}\boldsymbol{x}}{{\rm d}t} = \mathbf{M}^{-1}\boldsymbol{p}, \quad \frac{{\rm d}\boldsymbol{p}}{{\rm d}t} = -\frac{\partial \mathcal{L}}{\partial \mathbf{x}} - \lambda \frac{\partial c}{\partial \boldsymbol{x}} \quad\text{subject to}\quad c(\boldsymbol{x}) = 0 \quad\text{and}\quad \frac{\partial c}{\partial \boldsymbol{x}}\mathbf{M}^{-1}\boldsymbol{p} = 0$$ $\boldsymbol{x}(0) = \boldsymbol{x}_0,~\boldsymbol{p}(0) = \boldsymbol{p}_0$$c(\boldsymbol{x}_0) = 0$, define uma dinâmica hamiltoniana restrita que permanece confinada ao coletor de restrições, é reversível no tempo e conserva exatamente o elemento de volume Hamiltoniano e do coletor. Se usarmos um integrador simplético para sistemas hamiltonianos restritos, como SHAKE (2) ou RATTLE (3), que mantêm exatamente a restrição em cada passo de tempo resolvendo para o multiplicador Lagrange, podemos simular o avanço dinâmico exatoLtimesteps discretosδtde alguma restrição inicial satisfazendox$\left.\frac{\partial c}{\partial \boldsymbol{x}}\right|_{\boldsymbol{x}_0}\,\mathbf{M}^{-1}\boldsymbol{p}_0 = 0$$L$$\delta t$ aceitar o novo par de estados proposto x ′$\boldsymbol{x},\,\boldsymbol{p}$ com probabilidade mínima { 1 $\boldsymbol{x}',\,\boldsymbol{p}'$ Se intercalarmos essas atualizações dinâmicas com reamostragem parcial / total do momento da sua margem gaussiana (restrita ao subespaço linear definido por∂ $$\min\left\lbrace 1, \,\exp\left[ \mathcal{L}(\boldsymbol{x}) - \mathcal{L}(\boldsymbol{x}') + \frac{1}{2}\boldsymbol{p}^{\rm T}\mathbf{M}^{-1}\boldsymbol{p} - \frac{1}{2}\boldsymbol{p}'^{\rm T}\mathbf{M}^{-1}\boldsymbol{p}'\right] \right\rbrace.$$ ), então modulo a possibilidade da existência de múltiplos componentes múltiplos de restrição não-ligado, a dinâmica global MCMC deve ser ergodic e as amostras de estado configuraçãoxvai coverge na distribuição à densidade alvo restringido ao colector de restrição .$\frac{\partial c}{\partial \boldsymbol{x}}\mathbf{M}^{-1}\boldsymbol{p} = 0$$\boldsymbol{x}$

Para ver como o HMC restrito foi executado no caso aqui, executei a implementação do HMC restrito baseado em integrador geodésico descrito em (4) e disponível no Github aqui (divulgação completa: sou autor de (4) e proprietário do repositório do Github), que usa uma variação do esquema integrador 'geodésico-BAOAB' proposto em (5) sem a etapa estocástica de Ornstein-Uhlenbeck. Na minha experiência, esse esquema de integração geodésica é geralmente um pouco mais fácil de ajustar do que o esquema RATTLE usado em (1) devido à flexibilidade extra do uso de várias etapas internas menores para o movimento geodésico no coletor de restrições. Um notebook IPython que gera os resultados está disponível aqui .

Eu usei , μ = 1 e μ 0 = 2 . Um x inicial correspondente a um MLE de µ 0 foi encontrado pelo método de Newton (com a derivada de segunda ordem verificada para garantir que um máximo da probabilidade fosse encontrado). Corri uma dinâmica restrita com δ t = 0,5 , L = 5 intercalada com atualizações de momento completo para 1000 atualizações. O gráfico abaixo mostra os traços resultantes nos três componentes $N=3$$\mu=1$$\mu_0=2$$\boldsymbol{x}$$\mu_0$$\delta t = 0.5$$L=5$$\boldsymbol{x}$

e os valores correspondentes das derivadas de primeira e segunda ordem da probabilidade logarítmica negativa são mostrados abaixo

$\boldsymbol{x}$$\boldsymbol{x}$$\mathbb{R}^3$

$\epsilon \to 0$$\mathbb{R}^N$$|c(\boldsymbol{x})| < \epsilon$$\sqrt{\frac{\partial c}{\partial \boldsymbol{x}}^{\rm \scriptscriptstyle T}\frac{\partial c}{\partial \boldsymbol{x}}}$

Referências

1. MA Brubaker, M. Salzmann e R. Urtasun. Uma família de métodos MCMC em variedades definidas implicitamente. Em Anais da 15ª Conferência Internacional sobre Inteligência Artificial e Estatística , 2012.
   http://www.cs.toronto.edu/~mbrubake/projects/AISTATS12.pdf

2. J.-P. Ryckaert, G. Ciccotti e HJ Berendsen. Integração numérica das equações cartesianas de movimento de um sistema com restrições: dinâmica molecular de n-alcanos. Jornal de Física Computacional , 1977.
   http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.399.6868

3. HC Andersen. CHOCOLATE: Uma versão em "velocidade" do algoritmo SHAKE para cálculos de dinâmica molecular. Jornal de Física Computacional , 1983.
   http://www.sciencedirect.com/science/article/pii/0021999183900141

4. MM Graham e AJ Storkey. Inferência assintoticamente exata em modelos sem probabilidade. pré-impressão do arXiv arXiv: 1605.07826v3 , 2016.
   https://arxiv.org/abs/1605.07826

5. B. Leimkuhler e C. Matthews. Dinâmica molecular eficiente usando integração geodésica e divisão solvente-soluto. Proc. R. Soc. A. Vol. 472. No. 2189. The Royal Society , 2016.
   http://rspa.royalsocietypublishing.org/content/472/2189/20160138.abstract