Fiabilidad del modo de una muestra MCMC

En su libro Doing Bayesian Data Analysis, John Kruschke afirma que al usar JAGS de R

... la estimación del modo de una muestra de MCMC puede ser bastante inestable porque la estimación se basa en un algoritmo de suavizado que puede ser sensible a golpes y ondulaciones aleatorias en la muestra de MCMC. ( Análisis de datos bayesianos , página 205, sección 8.2.5.1)

Si bien entiendo el algoritmo Metropolis y las formas exactas como el muestreo de Gibbs, no estoy familiarizado con el algoritmo de suavizado que también se alude y por qué significaría que la estimación del modo de la muestra MCMC es inestable. ¿Alguien puede dar una idea intuitiva de lo que está haciendo el algoritmo de suavizado y por qué hace que la estimación del modo sea inestable?

No tengo el libro a la mano, así que no estoy seguro de qué método de suavizado utiliza Kruschke, pero por intuición considere este gráfico de 100 muestras de una normal estándar, junto con las estimaciones de densidad de kernel gaussianas que utilizan varios anchos de banda de 0.1 a 1.0. (Brevemente, los KDE gaussianos son una especie de histograma suavizado: estiman la densidad agregando un gaussiano para cada punto de datos, con la media en el valor observado).

Puede ver que incluso una vez que el suavizado crea una distribución unimodal, el modo generalmente está por debajo del valor conocido de 0.

Más, aquí hay una gráfica del modo estimado (eje y) por ancho de banda del kernel utilizado para estimar la densidad, usando la misma muestra. Esperemos que esto preste algo de intuición sobre cómo varía la estimación con los parámetros de suavizado.

#!/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 dio una buena respuesta; Así es como lo hacen los guiones R que vienen con el libro de Kruschke. La plotPost()función se define en el script R llamado DBDA2E-utilities.R. Muestra el modo estimado. Dentro de la definición de la función, hay estas dos líneas:

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

La density()función viene con el paquete de estadísticas base de R e implementa un filtro de densidad del núcleo del tipo que describió Sean Easter. Tiene argumentos opcionales para el ancho de banda del kernel de suavizado y para el tipo de kernel a utilizar. Su valor predeterminado es un núcleo gaussiano y tiene algo de magia interna para encontrar un buen ancho de banda. La density()función devuelve un objeto con un componente llamado yque tiene las densidades suavizadas en varios valores x. La segunda línea de código, arriba, solo encuentra el xvalor donde yes máximo.