simulation d'échantillons aléatoires avec un MLE donné

Cette question de validation croisée portant sur la simulation d'un échantillon conditionnel à une somme fixe m'a rappelé un problème posé par George Casella .

Étant donné un modèle paramétrique $f(x|\theta)$ et un échantillon iid de ce modèle , le MLE de est donné par Pour une valeur donnée de \ theta , existe-t-il un moyen générique de simuler un échantillon iid (X_1, \ ldots, X_n) conditionnel sur la valeur du MLE \ hat {\ theta} (X_1, \ ldots, X_n) ?$(X_1,\ldots,X_n)$θ ( x 1 , ... , x n ) = arg min n Σ i = 1 log f ( x i | θ ) θ ( X 1 , ... , X n ) θ ( X 1 , … , X n$\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)$

Par exemple, prenez une distribution $\mathfrak{T}_5$ , avec le paramètre d'emplacement $\mu$ , dont la densité est

$$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)$$ comment simuler $(X_1,\ldots,X_n)$ conditionnellement à $\hat{\mu}(X_1,\ldots,X_n)=\mu_0$ ? Dans cet exemple $\mathfrak{T}_5$ , la distribution de $\hat{\mu}(X_1,\ldots,X_n)$ n'a pas d'expression de forme fermée.

Une option consisterait à utiliser une variante HMC contrainte comme décrit dans A Family of MCMC Methods on Implicitly Defined Manifolds par Brubaker et al (1). Cela nécessite que nous puissions exprimer la condition que l'estimation du maximum de vraisemblance du paramètre de localisation est égale à une certaine fixe comme une contrainte holonomique implicitement définie (et différenciable) . Nous pouvons ensuite simuler une dynamique hamiltonienne contrainte soumise à cette contrainte, et accepter / rejeter au sein d'une étape Metropolis-Hastings comme en HMC standard. c ( { x i } N i = 1 ) = $\mu_0$$c\left(\lbrace x_i \rbrace_{i=1}^N\right) = 0$

La log-vraisemblance négative est qui a des dérivées partielles de premier et deuxième ordre par rapport à le paramètre d'emplacement Une estimation du maximum de vraisemblance de est alors implicitement définie comme une solution à μ∂

$$\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.$$

Je ne sais pas s'il y a des résultats suggérant qu'il y aura un MLE unique pour pour un - la densité n'est pas log-concave dans donc il ne semble pas trivial pour garantir cela. S'il n'y a qu'une seule solution unique, ce qui précède définit implicitement une variété dimensionnelle connectée intégrée dans correspondant à l'ensemble de avec MLE pour égal à$\mu$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$N - 1$$\mathbb{R}^N$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$\mu_0$. S'il existe plusieurs solutions, le collecteur peut être composé de plusieurs composants non connectés, dont certains peuvent correspondre à des minima dans la fonction de vraisemblance. Dans ce cas, nous aurions besoin d'un mécanisme supplémentaire pour se déplacer entre les composants non connectés (car la dynamique simulée restera généralement confinée à un seul composant) et vérifier la condition de second ordre et rejeter un mouvement s'il correspond à un déplacement vers un minimum dans la vraisemblance.

Si nous utilisons pour désigner le vecteur et introduisons un état de moment conjugué avec la matrice de masse et une Lagrange multiplicateur pour la contrainte scalaire puis la solution au système des ODE $\boldsymbol{x}$$\left[ x_1 \dots x_N\right]^{\rm T}$$\boldsymbol{p}$$\mathbf{M}$$\lambda$$c(\boldsymbol{x})$

$$\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$$ étant donné la condition initiale avec et , définit une dynamique hamiltonienne contrainte qui reste confinée au collecteur de contraintes, est réversible dans le temps et conserve exactement le hamiltonien et l'élément de volume du collecteur. Si nous utilisons un intégrateur symplectique pour les systèmes hamiltoniens contraints tels que SHAKE (2) ou RATTLE (3), qui maintiennent exactement la contrainte à chaque pas de temps en résolvant pour le multiplicateur de Lagrange, nous pouvons simuler les pas de temps discrets dynamiques vers l'avant exacts$\boldsymbol{x}(0) = \boldsymbol{x}_0,~\boldsymbol{p}(0) = \boldsymbol{p}_0$$c(\boldsymbol{x}_0) = 0$$\left.\frac{\partial c}{\partial \boldsymbol{x}}\right|_{\boldsymbol{x}_0}\,\mathbf{M}^{-1}\boldsymbol{p}_0 = 0$$L$$\delta t$à partir d'une certaine contrainte initiale satisfaisant et accepter la nouvelle paire d'états proposée avec probabilité Si nous entrelacons ces mises à jour dynamiques avec un rééchantillonnage partiel / complet des impulsions de leur marginal gaussien (limité au sous-espace linéaire défini par$\boldsymbol{x},\,\boldsymbol{p}$$\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$) modulo alors la possibilité qu'il y ait plusieurs composants de collecteur de contraintes non connectés, la dynamique MCMC globale devrait être ergodique et les échantillons d'état de configuration couvriront la distribution à la densité cible limitée au collecteur de contraintes.$\boldsymbol{x}$

Pour voir les performances de la console HMC contrainte dans le cas présent, j'ai exécuté l'implémentation de la console HMC contrainte basée sur l'intégrateur géodésique décrite dans (4) et disponible sur Github ici (divulgation complète: je suis l'auteur de (4) et propriétaire du référentiel Github), qui utilise une variante du schéma d'intégrateur «géodésique-BAOAB» proposé en (5) sans l'étape stochastique d'Ornstein-Uhlenbeck. D'après mon expérience, ce schéma d'intégration géodésique est généralement un peu plus facile à régler que le schéma RATTLE utilisé dans (1) en raison de la flexibilité supplémentaire d'utiliser plusieurs étapes internes plus petites pour le mouvement géodésique sur le collecteur de contraintes. Un bloc-notes IPython générant les résultats est disponible ici .

J'ai utilisé , et . Un initial correspondant à un MLE de été trouvé par la méthode de Newton (avec la dérivée du second ordre vérifiée pour s'assurer qu'un maximum de vraisemblance a été trouvé). J'ai exécuté une dynamique contrainte avec , entrelacé avec des rafraîchissements de momentum complet pour 1000 mises à jour. Le graphique ci-dessous montre les traces résultantes sur les trois composants$N=3$$\mu=1$$\mu_0=2$$\boldsymbol{x}$$\mu_0$$\delta t = 0.5$$L=5$$\boldsymbol{x}$

et les valeurs correspondantes des dérivées du premier et du second ordre de la log-vraisemblance négative sont indiquées ci-dessous

à partir de laquelle on peut voir que nous sommes au maximum de la log-vraisemblance pour tous les échantillonnés . Bien qu'il ne soit pas facilement apparent à partir des tracés de trace individuels, le échantillonné se trouve sur une variété non linéaire 2D intégrée dans - l'animation ci-dessous montre les échantillons en 3D$\boldsymbol{x}$$\boldsymbol{x}$$\mathbb{R}^3$

Selon l'interprétation de la contrainte, il peut également être nécessaire d'ajuster la densité cible par un facteur jacobien comme décrit dans (4). En particulier, si nous voulons des résultats cohérents avec la limite d'utiliser une approche de type ABC pour maintenir approximativement la contrainte en proposant des mouvements non contraints dans et en acceptant if , alors nous devons multiplier la densité cible par . Dans l'exemple ci-dessus, je n'ai pas inclus cet ajustement, de sorte que les échantillons proviennent de la densité cible d'origine limitée au collecteur de contraintes.$\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}}}$

Les références

1. MA Brubaker, M. Salzmann et R. Urtasun. Une famille de méthodes MCMC sur des variétés implicitement définies. Dans les actes de la 15e Conférence internationale sur l'intelligence artificielle et les statistiques , 2012.
   http://www.cs.toronto.edu/~mbrubake/projects/AISTATS12.pdf

2. J.-P. Ryckaert, G. Ciccotti et HJ Berendsen. Intégration numérique des équations cartésiennes de mouvement d'un système à contraintes: dynamique moléculaire des n-alcanes. Journal of Computational Physics , 1977.
   http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.399.6868

3. HC Andersen. RATTLE: Une version "vitesse" de l'algorithme SHAKE pour les calculs de dynamique moléculaire. Journal of Computational Physics , 1983.
   http://www.sciencedirect.com/science/article/pii/0021999183900141

4. MM Graham et AJ Storkey. Inférence asymptotiquement exacte dans les modèles sans vraisemblance. arXiv pré-impression arXiv: 1605.07826v3 , 2016.
   https://arxiv.org/abs/1605.07826

5. B. Leimkuhler et C. Matthews. Dynamique moléculaire efficace utilisant l'intégration géodésique et la séparation solvant-soluté. Proc. R. Soc. A. Vol. 472. No 2189. The Royal Society , 2016.
   http://rspa.royalsocietypublishing.org/content/472/2189/20160138.abstract