Fiabilité du mode à partir d'un échantillon MCMC

12

Dans son livre Doing Bayesian Data Analysis, John Kruschke déclare qu'en utilisant JAGS de R

... l'estimation du mode à partir d'un échantillon MCMC peut être plutôt instable car l'estimation est basée sur un algorithme de lissage qui peut être sensible aux bosses et ondulations aléatoires dans l'échantillon MCMC. ( Faire une analyse des données bayésiennes , page 205, section 8.2.5.1)

Bien que je connaisse l'algorithme de Metropolis et les formes exactes comme l'échantillonnage de Gibbs, je ne connais pas l'algorithme de lissage qui y est également mentionné et pourquoi cela signifierait que l'estimation du mode à partir de l'échantillon MCMC est instable. Quelqu'un peut-il donner un aperçu intuitif de ce que fait l'algorithme de lissage et pourquoi il rend l'instabilité du mode instable?

Morgan Ball
la source
2
Je pense que John Kruschke parle de l'algorithme d'estimation de mode basé sur l'estimation de la densité du noyau.
Andrey Kolyadin
2
Ce lien peut être utile.
Andrey Kolyadin
À moins que je ne me trompe étant nouveau dans ce domaine de la statistique, JAGS génère un ensemble d'échantillons à partir de la distribution postérieure plutôt qu'une fonction de densité de probabilité, donc je ne suis pas sûr que l'estimation de la densité kernal entre en jeu. Merci pour le lien cependant.
Morgan Ball
Je pense que cela est probablement plus lié à la façon dont vous obtenez un mode à partir d'un grand échantillon d'une variable continue où il ne peut y avoir plus d'une valeur particulière et vous devez donc grouper (ou lisser) l'échantillon.
Morgan Ball
1
Vous pouvez obtenir le mode en tant que valeur avec une densité maximale lors de l'estimation de la densité du noyau. (Au moins, c'est ce que je fais, et si je ne me trompe pas, J. Kruschke utilise la même approche dans ses exemples)
Andrey Kolyadin

Réponses:

8

Je n'ai pas le livre à portée de main, donc je ne sais pas quelle méthode de lissage Kruschke utilise, mais pour l'intuition, considérez ce tracé de 100 échantillons à partir d'une normale standard, ainsi que les estimations de la densité du noyau gaussien en utilisant différentes bandes passantes de 0,1 à 1,0. (En bref, les KDE gaussiens sont une sorte d'histogramme lissé: ils estiment la densité en ajoutant un gaussien pour chaque point de données, avec une moyenne à la valeur observée.)

Vous pouvez voir que même une fois que le lissage crée une distribution unimodale, le mode est généralement inférieur à la valeur connue de 0.

entrez la description de l'image ici

De plus, voici un tracé du mode estimé (axe y) par la bande passante du noyau utilisée pour estimer la densité, en utilisant le même échantillon. Espérons que cela donne une idée de la façon dont l'estimation varie en fonction des paramètres de lissage.

entrez la description de l'image ici

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Feb  1 09:35:51 2017

@author: seaneaster
"""

import numpy as np
from matplotlib import pylab as plt
from sklearn.neighbors import KernelDensity

REAL_MODE = 0
np.random.seed(123)

def estimate_mode(X, bandwidth = 0.75):
    kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
    u = np.linspace(-3,3,num=1000)[:, np.newaxis]
    log_density = kde.score_samples(u)
    return u[np.argmax(log_density)]

X = np.random.normal(REAL_MODE, size = 100)[:, np.newaxis] # keeping to standard normal

bandwidths = np.linspace(0.1, 1., num = 8)

plt.figure(0)
plt.hist(X, bins = 100, normed = True, alpha = 0.25)

for bandwidth in bandwidths:
    kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
    u = np.linspace(-3,3,num=1000)[:, np.newaxis]
    log_density = kde.score_samples(u)
    plt.plot(u, np.exp(log_density))

bandwidths = np.linspace(0.1, 3., num = 100)
modes = [estimate_mode(X, bandwidth) for bandwidth in bandwidths]
plt.figure(1)
plt.plot(bandwidths, np.array(modes))
Sean Easter
la source
5

Sean Easter a fourni une belle réponse; voici comment les scripts R sont fournis avec le livre de Kruschke. La plotPost()fonction est définie dans le script R nommé DBDA2E-utilities.R. Il affiche le mode estimé. À l'intérieur de la définition de fonction, il y a ces deux lignes:

mcmcDensity = density(paramSampleVec)
mo = mcmcDensity$x[which.max(mcmcDensity$y)]

La density()fonction est livrée avec le package de statistiques de base de R et implémente un filtre de densité du noyau du type décrit par Sean Easter. Il a des arguments facultatifs pour la bande passante du noyau de lissage et pour le type de noyau à utiliser. Il utilise par défaut un noyau gaussien et possède une certaine magie interne pour trouver une belle bande passante. La density()fonction renvoie un objet avec un composant nommé yqui a les densités lissées à différentes valeurs x. La deuxième ligne de code, ci-dessus, trouve juste la xvaleur où yest maximum.

John K. Kruschke
la source