simulando muestras aleatorias con un MLE dado

Esta pregunta de validación cruzada sobre la simulación de una muestra condicional a tener una suma fija me recordó un problema que George Casella me planteó .

Dado un modelo paramétrico y una muestra de iid de este modelo , el MLE de viene dado por Para un valor dado de \ theta , ¿hay una forma genérica de simular una muestra de iid (X_1, \ ldots, X_n) condicional al valor de MLE \ hat {\ theta} (X_1, \ ldots, X_n) ?$f(x|\theta)$$(X_1,\ldots,X_n)$θ ( x 1 , ... , x n ) = arg min n Σ i = 1 log f ( x i | θ ) θ ( X 1 , ... , X n ) θ ( X 1 , ... $\theta$

$$\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 ejemplo, tome una distribución $\mathfrak{T}_5$ , con el parámetro de ubicación $\mu$ , cuya densidad es

$$f(x|\mu)=\dfrac{\Gamma(3)}{\Gamma(1/2)\Gamma(5/2)}\,\left[1+(x-\mu)^2/5\right]^{-3}$$ If $$(X_1,\ldots,X_n)\stackrel{\text{iid}}{\sim} f(x|\mu)$$ ¿cómo podemos simular $(X_1,\ldots,X_n)$ condicional en $\hat{\mu}(X_1,\ldots,X_n)=\mu_0$ ? En este ejemplo de $\mathfrak{T}_5$ , la distribución de $\hat{\mu}(X_1,\ldots,X_n)$ no tiene una expresión de forma cerrada.

Una opción sería utilizar una variante HMC restringida como se describe en Una familia de métodos MCMC en manifiestos definidos implícitamente por Brubaker et al (1). Esto requiere que podamos expresar la condición de que la estimación de máxima verosimilitud del parámetro de ubicación sea igual a alguna fija como alguna restricción holonómica implícitamente definida (y diferenciable) . Entonces podemos simular una dinámica hamiltoniana restringida sujeta a esta restricción, y aceptar / rechazar dentro de un paso Metropolis-Hastings como en la HMC estándar. c ( { x i } N i = 1 ) = $\mu_0$$c\left(\lbrace x_i \rbrace_{i=1}^N\right) = 0$

La probabilidad de registro negativa es que tiene derivadas parciales de primer y segundo orden con respecto a el parámetro de ubicación Una estimación de máxima verosimilitud de se define implícitamente como una solución para μ∂

$$\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$ μ0c=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.$$

No estoy seguro de si hay resultados que sugieran que habrá un MLE único para para dado - la densidad no es cóncava en por lo que no parece trivial para garantizar esto. Si hay una única solución única, lo anterior define implícitamente una variedad dimensional conectada incrustada en correspondiente al conjunto de con MLE para igual to{ x i } N i = 1 μ N - 1 R N { x i } N i = 1 μ μ $\mu$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$N - 1$$\mathbb{R}^N$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$\mu_0$. Si hay múltiples soluciones, entonces el múltiple puede consistir en múltiples componentes no conectados, algunos de los cuales pueden corresponder a mínimos en la función de probabilidad. En este caso, necesitaríamos tener algún mecanismo adicional para moverse entre los componentes no conectados (ya que la dinámica simulada generalmente permanecerá confinada a un solo componente) y verificar la condición de segundo orden y rechazar un movimiento si corresponde al movimiento a un mínimo en la probabilidad.

Si usamos para denotar el vector e introducimos un estado de momento conjugado con matriz de masa y un Lagrange multiplicador para la restricción escalar luego la solución al sistema de EDO [ x 1 … x N ] T p M λ c ( x ) d$\boldsymbol{x}$$\left[ x_1 \dots x_N\right]^{\rm T}$$\boldsymbol{p}$$\mathbf{M}$$\lambda$$c(\boldsymbol{x})$x(0)=x0,p(0)=p0c(x0)=0 ∂

$$\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$$ condición inicial dada con y , define una dinámica hamiltoniana restringida que permanece confinada al múltiple de restricción, es reversible en el tiempo y conserva exactamente el elemento de volumen hamiltoniano y múltiple. Si utilizamos un integrador simpléctico para sistemas hamiltonianos restringidos como SHAKE (2) o RATTLE (3), que mantienen exactamente la restricción en cada paso de tiempo resolviendo para el multiplicador de Lagrange, podemos simular el paso de tiempo discreto directo dinámico$\boldsymbol{x}(0) = \boldsymbol{x}_0,~\boldsymbol{p}(0) = \boldsymbol{p}_0$$c(\boldsymbol{x}_0) = 0$Lδtx$\left.\frac{\partial c}{\partial \boldsymbol{x}}\right|_{\boldsymbol{x}_0}\,\mathbf{M}^{-1}\boldsymbol{p}_0 = 0$$L$$\delta t$de alguna restricción inicial que satisfaga y acepte el nuevo par de estados propuesto con probabilidad Si intercalamos estas actualizaciones dinámicas con remuestreo parcial / total de los momentos de su marginal gaussiano (restringido al subespacio lineal definido porx ′$\boldsymbol{x},\,\boldsymbol{p}$ min { 1 $\boldsymbol{x}',\,\boldsymbol{p}'$ $$\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.$$ $\frac{\partial c}{\partial \boldsymbol{x}}\mathbf{M}^{-1}\boldsymbol{p} = 0$) luego module la posibilidad de que existan múltiples componentes del múltiple de restricción no conectados, la dinámica general de MCMC debe ser ergódica y el estado de configuración samples cubrirá en distribución a la densidad objetivo restringida al múltiple de restricción.$\boldsymbol{x}$

Para ver cómo funcionó HMC restringido para el caso aquí, ejecuté la implementación de HMC restringida basada en integrador geodésico descrita en (4) y disponible en Github aquí (divulgación completa: soy autor de (4) y propietario del repositorio de Github), que utiliza una variación del esquema integrador 'geodésico-BAOAB' propuesto en (5) sin el paso estocástico Ornstein-Uhlenbeck. En mi experiencia, este esquema de integración geodésica es generalmente un poco más fácil de ajustar que el esquema RATTLE utilizado en (1) debido a la flexibilidad adicional de usar múltiples pasos internos más pequeños para el movimiento geodésico en el múltiple de restricción. Un cuaderno de IPython que genera los resultados está disponible aquí .

Usé , y . El método de Newton encontró un inicial correspondiente a un MLE de (con la derivada de segundo orden verificada para asegurar que se encontró un máximo de la probabilidad). Ejecuté una dinámica restringida con , intercalada con actualizaciones completas de impulso para 1000 actualizaciones. La siguiente gráfica muestra las trazas resultantes en los tres componentes$N=3$$\mu=1$$\mu_0=2$$\boldsymbol{x}$$\mu_0$$\delta t = 0.5$$L=5$$\boldsymbol{x}$

y los valores correspondientes de las derivadas de primer y segundo orden de la probabilidad de registro negativa se muestran a continuación

de lo cual se puede ver que tenemos un máximo de log-verosimilitud para todos los muestreados . Aunque no es evidente a partir de las trazas de trazos individuales, el muestreado se encuentra en un múltiple no lineal 2D incrustado en - la animación a continuación muestra las muestras en 3D$\boldsymbol{x}$$\boldsymbol{x}$$\mathbb{R}^3$

Dependiendo de la interpretación de la restricción, también puede ser necesario ajustar la densidad objetivo por algún factor jacobiano como se describe en (4). En particular, si queremos resultados consistentes con el límite de de usar un enfoque similar a ABC para mantener aproximadamente la restricción proponiendo movimientos sin restricciones en y aceptando if , entonces necesitamos multiplicar la densidad objetivo por . En el ejemplo anterior, no incluí este ajuste, por lo que las muestras provienen de la densidad objetivo original restringida al múltiple de restricción.$\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}}}$

Referencias

1. MA Brubaker, M. Salzmann y R. Urtasun. Una familia de métodos MCMC en múltiples definidos implícitamente. En Actas de la 15ª Conferencia Internacional sobre Inteligencia Artificial y Estadísticas , 2012.
   http://www.cs.toronto.edu/~mbrubake/projects/AISTATS12.pdf

2. J.-P. Ryckaert, G. Ciccotti y HJ Berendsen. Integración numérica de las ecuaciones cartesianas de movimiento de un sistema con restricciones: dinámica molecular de n-alcanos. Journal of Computational Physics , 1977.
   http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.399.6868

3. HC Andersen. RATTLE: Una versión de "velocidad" del algoritmo SHAKE para cálculos de dinámica molecular. Journal of Computational Physics , 1983.
   http://www.sciencedirect.com/science/article/pii/0021999183900141

4. MM Graham y AJ Storkey. Inferencia asintóticamente exacta en modelos libres de probabilidad. pre-impresión de arXiv arXiv: 1605.07826v3 , 2016.
   https://arxiv.org/abs/1605.07826

5. B. Leimkuhler y C. Matthews. Dinámica molecular eficiente utilizando integración geodésica y división solvente-soluto. Proc. R. Soc. A. vol. 472. No. 2189. The Royal Society , 2016.
   http://rspa.royalsocietypublishing.org/content/472/2189/20160138.abstract