Lors de l'approximation d'un postérieur à l'aide de MCMC, pourquoi ne sauvegardons-nous pas les probabilités postérieures mais utilisons-nous ensuite les fréquences des valeurs des paramètres?

J'évalue actuellement les paramètres d'un modèle défini par plusieurs équations différentielles ordinaires (ODE). J'essaie ceci avec une approche bayésienne en approximant la distribution postérieure des paramètres étant donné certaines données en utilisant la chaîne de Markov Monte Carlo (MCMC).

Un échantillonneur MCMC génère une chaîne de valeurs de paramètres où il utilise la probabilité postérieure (non normalisée) d'une certaine valeur de paramètre pour décider (stochastiquement) s'il ajoutera cette valeur à la chaîne ou ajoutera à nouveau la valeur précédente. Mais, il semble être la pratique que les probabilités postérieures réelles n'ont pas besoin d'être enregistrées, mais plutôt un histogramme à n dimensions des valeurs de paramètres résultantes générées et des statistiques récapitulatives comme les régions de densité la plus élevée (HDR) d'une distribution postérieure des paramètres est calculée à partir de cet histogramme. C'est du moins ce que je pense avoir appris du livre de didacticiel de Kruschkes sur l'inférence bayésienne .

Ma question: ne serait-il pas plus simple de sauvegarder les probabilités postérieures des valeurs des paramètres échantillonnés avec celles-ci et d'approximer la distribution postérieure à partir de ces valeurs et non à partir des fréquences des valeurs des paramètres dans la chaîne MCMC? Le problème de la phase de rodage ne se poserait pas car l'échantillonneur échantillonnerait encore plus souvent des régions à faible probabilité qu'il ne le mériterait par leurs probabilités postérieures, mais ce ne serait plus le problème de leur donner des valeurs de probabilité indûment élevées.

Ceci est une question intéressante, avec différents problèmes:

1. Les algorithmes MCMC ne recyclent pas toujours le calcul de la densité postérieure à toutes les valeurs proposées, mais certaines techniques de réduction de variance comme Rao-Blackwellisation le font. Par exemple, dans un article de Biometrika de 1996 avec George Casella, nous proposons d'utiliser toutes les valeurs simulées,$\theta_i$ $(i=1,\ldots,T)$, acceptés ou non, en introduisant des poids $\omega_i$ qui tournent la moyenne $$\sum_{i=1}^T \omega_ih(\theta_i)\big/\sum_{i=1}^T \omega_i$$ dans un estimateur presque sans biais. (Le presque étant dû à la normalisation par la somme des poids.)
2. MCMC est souvent utilisé sur des problèmes de grande dimension (paramètre). Proposer une approximation de l'ensemble postérieur sur la base des valeurs de densité observées à certaines valeurs de paramètres est tout un défi, y compris la question de la constante de normalisation mentionnée dans la réponse et les commentaires de Tim. On peut imaginer une approche qui est un mélange d'estimation non paramétrique du noyau (comme par exemple le krigging ) et de régression, mais les experts avec lesquels j'ai discuté de cette solution [il y a quelques années] étaient assez sceptiques. Le problème est que l'estimateur résultant reste non paramétrique et "jouit" donc de vitesses de convergence non paramétriques qui sont plus lentes que les vitesses de convergence de Monte-Carlo, plus la dimension est pire, plus la dimension est grande.
3. Une autre utilisation potentielle de la disponibilité des valeurs postérieures est de pondérer chaque valeur simulée par son postérieur associé, comme dans Malheureusement, cela crée un biais car les valeurs simulées sont déjà simulées à partir de la partie postérieure: Même sans problème de normalisation, ces simulations devraient donc cibler et utiliser un poids proportionnel à$\pi(\theta|\mathcal{D})$ $$\frac{1}{T}\sum_{t=1}^T h(\theta_t) \pi(\theta_t|\mathcal{D})$$ $$\mathbb{E}[h(\theta_t) \pi(\theta_t|\mathcal{D})]=\int h(\theta) h(\theta_t) \pi(\theta_t|\mathcal{D})^2\text{d}\theta$$ $\pi(\theta|\mathcal{D})^{1/2}$$\pi(\theta|\mathcal{D})^{1/2}$mais je ne connais pas de résultats prônant ce changement de cible. Comme vous le mentionnez dans les commentaires, cela est lié au revenu en ce sens que toutes les simulations produites dans un cycle de revenu simulé peuvent être recyclées à des fins de Monte Carlo (intégration) de cette façon. Un problème numérique, cependant, est de gérer plusieurs fonctions d'importance de la forme avec des constantes de normalisation manquantes.$\pi(\theta)^{1/T}$

Comme vous l'avez bien remarqué, les probabilités dont nous traitons ne sont pas normalisées . Fondamentalement, nous utilisons MCMC pour calculer le facteur de normalisation dans le théorème de Bayes. Nous ne pouvons pas utiliser les probabilités car elles ne sont pas normalisées. La procédure que vous proposez: enregistrer les probabilités non normalisées puis les diviser par leur somme est incorrecte.

Permettez-moi de vous le montrer par l'exemple. Imaginez que vous utilisiez Monte Carlo pour tirer dix valeurs de la distribution de Bernoulli paramétrée par$p=0.9$, ils sont les suivants:

1 0 1 1 1 1 1 1 1 1

vous avez également des probabilités correspondantes:

0.9 0.1 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9

Dans ce cas, les probabilités sont normalisées, mais les diviser par leur somme (celle des axiomes de probabilité est égale à l'unité) ne devrait rien changer. Unfortunatelly, en utilisant votre procédure , elle ne change les résultats:

> f/sum(f)
 [1] 0.10975610 0.01219512 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610

Pourquoi donc? La réponse est simple, dans votre échantillon, chaque "probabilité" enregistrée fapparaît avec une probabilité f, vous pondérez donc les probabilités par elles-mêmes!