Étant donné une chaîne 10MC MCMC, comment puis-je déterminer son ou ses modes postérieurs dans R?

10

Question: Avec une chaîne MCMC à 10 dimensions, disons que je suis prêt à vous remettre une matrice des tirages: 100 000 itérations (lignes) par 10 paramètres (colonnes), comment identifier au mieux les modes postérieurs? Je suis particulièrement préoccupé par plusieurs modes.

Contexte:Je me considère comme un statisticien averti en calcul, mais quand un collègue m'a posé cette question, j'avais honte de ne pas pouvoir trouver de réponse raisonnable. La principale préoccupation est que plusieurs modes peuvent apparaître, mais seulement si au moins huit des dix dimensions sont prises en compte. Ma première pensée serait d'utiliser une estimation de la densité du noyau, mais une recherche dans R n'a rien révélé de prometteur pour des problèmes de plus de trois dimensions. Le collègue a proposé une stratégie de binning ad hoc en dix dimensions et en recherchant un maximum, mais je crains que la bande passante ne conduise à des problèmes de rareté importants ou à un manque de résolution pour discerner plusieurs modes. Cela dit, j'accepterais volontiers des suggestions de suggestions de bande passante automatisée, des liens vers un estimateur de densité de 10 noyaux ou tout autre élément que vous connaissez.

Préoccupations:

  1. Nous pensons que la distribution peut être assez biaisée; par conséquent, nous souhaitons identifier le ou les modes postérieurs et non les moyens postérieurs.

  2. Nous craignons qu'il puisse exister plusieurs modes postérieurs.

  3. Si possible, nous préférons une suggestion basée sur R. Mais tout algorithme fera l'affaire tant qu'il n'est pas incroyablement difficile à mettre en œuvre. Je suppose que je préférerais ne pas implémenter un estimateur de densité de noyau Nd avec une sélection automatisée de la bande passante à partir de zéro.

M. Tibbits
la source
Veuillez consulter le thème sur les méthodes d'estimation du mode rapide stats.stackexchange.com/questions/33625
Pavel Ruzankin

Réponses:

9

Avez-vous envisagé d'utiliser une approche de voisin le plus proche?

Par exemple, construire une liste des kvoisins les plus proches pour chacun des 100'000 points, puis considérer le point de données avec la plus petite distance du kthvoisin un mode. En d'autres termes: trouvez le point avec la «plus petite bulle» contenant d' kautres points autour de ce point.

Je ne sais pas à quel point cela est robuste et le choix pour kinfluencer évidemment les résultats.

Andre Holzner
la source
Parfois, je veux juste me renverser la tête. Excellente suggestion.
M. Tibbits
1
J'ai aussi pensé à utiliser la kmeansfonction dans R. Je ne devrais vraiment pas poser de questions entre minuit et 4 heures du matin.
M. Tibbits
4

Ce n'est qu'une réponse partielle.

J'ai récemment utilisé figtree pour des estimations multidimensionnelles de la densité du noyau. C'est un paquet C et je l'ai fait fonctionner assez facilement. Cependant, je ne l'ai utilisé que pour estimer la densité à des points particuliers, pas pour calculer des statistiques sommaires.

csgillespie
la source
3

Si vous conservez les probabilités du journal, vous pouvez simplement sélectionner celle avec la valeur la plus élevée. De plus, si votre intérêt est principalement le mode, il suffit de faire une optimisation pour trouver le point avec la probabilité de journal la plus élevée.

John Salvatier
la source
C'est la réponse la plus pertinente, au moins la première partie! Dans de nombreuses simulations MCMC, les probabilités (log) sont calculées pour toutes les propositions et peuvent donc être stockées. Ou la valeur la plus élevée jusqu'à présent et son argument peut être stocké. À condition que l'algorithme MCMC ait convergé sur le nombre de simulations que vous avez exécutées, il s'agit d'une approche valide.
Xi'an
2

Avez-vous pensé à «PRIM / bump chasse»? (voir par exemple la section 9.3 de «Les éléments de l'apprentissage statistique» de Tibshirani et al. ou demandez à votre moteur de recherche préféré). Je ne sais pas si c'est implémenté dans R cependant.

[D'après ce que j'ai compris, essayez-vous de trouver le mode de densité de probabilité à partir duquel vos 100 000 lignes sont tirées. Votre problème serait donc partiellement résolu en trouvant une density estimationméthode appropriée ].

Andre Holzner
la source
Oui, il y a un paquet prim , avec une vignette R: Utilisation de prim pour la chasse aux bosses . Cependant, je ne sais pas comment cela fonctionnera dans ce cas.
chl