Confiabilidade do modo de uma amostra do MCMC

Em seu livro Doing Bayesian Data Analysis, John Kruschke afirma que, ao usar JAGS de R

... a estimativa do modo de uma amostra do MCMC pode ser bastante instável porque a estimativa é baseada em um algoritmo de suavização que pode ser sensível a solavancos e ondulações aleatórios na amostra do MCMC. ( Fazendo a análise de dados bayesiana , página 205, seção 8.2.5.1)

Embora eu tenha uma noção do algoritmo Metropolis e formas exatas como a amostragem de Gibbs, não estou familiarizado com o algoritmo de suavização que também é mencionado e por que isso significaria que a estimativa do modo da amostra do MCMC é instável. Alguém é capaz de fornecer uma visão intuitiva do que o algoritmo de suavização está fazendo e por que torna instável a estimativa do modo?

Eu não tenho o livro em mãos, então não tenho certeza do método de suavização que Kruschke usa, mas, para intuição, considere esse gráfico de 100 amostras de um normal padrão, juntamente com as estimativas de densidade de kernel gaussiana usando várias larguras de banda de 0,1 a 1,0. (Resumidamente, os KDEs gaussianos são uma espécie de histograma suavizado: eles estimam a densidade adicionando um gaussiano para cada ponto de dados, com média no valor observado.)

Você pode ver que, mesmo quando a suavização cria uma distribuição unimodal, o modo geralmente fica abaixo do valor conhecido de 0.

Além disso, aqui está um gráfico do modo estimado (eixo y) por largura de banda do kernel usado para estimar a densidade, usando a mesma amostra. Esperemos que isso dê alguma intuição sobre como a estimativa varia com os parâmetros de suavização.

#!/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 deu uma boa resposta; aqui está como é realmente feito pelos scripts R que acompanham o livro de Kruschke. A plotPost()função é definida no script R nomeado DBDA2E-utilities.R. Ele exibe o modo estimado. Dentro da definição da função, existem estas duas linhas:

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

A density()função vem com o pacote de estatísticas básicas de R e implementa um filtro de densidade de kernel do tipo descrito por Sean Easter. Possui argumentos opcionais para a largura de banda do kernel de suavização e para o tipo de kernel a ser usado. O padrão é um kernel gaussiano e possui alguma mágica interna para encontrar uma boa largura de banda. A density()função retorna um objeto com um componente nomeado yque possui as densidades suavizadas em vários valores x. A segunda linha de código, acima, apenas encontra o xvalor onde yé máximo.