Wiarygodność trybu z próbki MCMC

W swojej książce Doing Bayesian Data Analysis John Kruschke stwierdza, że ​​używając JAGS z R.

... oszacowanie trybu z próbki MCMC może być raczej niestabilne, ponieważ oszacowanie opiera się na algorytmie wygładzania, który może być wrażliwy na przypadkowe nierówności i pomarszczenia w próbce MCMC. ( Doing Bayesian Data Analysis , strona 205, sekcja 8.2.5.1)

Chociaż dobrze rozumiem algorytm Metropolis i dokładne formy, takie jak próbkowanie Gibbsa, nie znam również algorytmu wygładzania, o którym również wspomniano, i dlaczego oznaczałoby to, że oszacowanie trybu z próbki MCMC jest niestabilne. Czy ktokolwiek jest w stanie dać intuicyjny wgląd w to, co robi algorytm wygładzania i dlaczego powoduje, że oszacowanie trybu jest niestabilne?

Nie mam pod ręką tej książki, więc nie jestem pewien, jakiej metody wygładzania używa Kruschke, ale dla intuicji rozważ ten wykres 100 próbek ze standardowej normy, wraz z szacunkami gęstości jądra Gaussa przy użyciu różnych szerokości pasma od 0,1 do 1,0. (W skrócie, Gaussowskie KDE są rodzajem wygładzonego histogramu: szacują gęstość poprzez dodanie Gaussa dla każdego punktu danych, ze średnią dla obserwowanej wartości.)

Widać, że nawet gdy wygładzenie tworzy rozkład nieimodalny, tryb jest generalnie poniżej znanej wartości 0.

Więcej, oto wykres trybu szacunkowego (oś y) według przepustowości jądra użytej do oszacowania gęstości, przy użyciu tej samej próbki. Miejmy nadzieję, że pozwala to na intuicyjną zmianę szacunku w zależności od parametrów wygładzania.

#!/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 podał miłą odpowiedź; oto, jak to faktycznie robi skrypty R dołączone do książki Kruschke. plotPost()Funkcja jest zdefiniowana w skrypcie R nazwie DBDA2E-utilities.R. Wyświetla tryb szacunkowy. W definicji funkcji znajdują się dwie linie:

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

Ta density()funkcja jest dostarczana z pakietem podstawowych statystyk R i implementuje filtr gęstości jądra opisany przez Sean Easter. Ma opcjonalne argumenty za przepustowością jądra wygładzającego i za typem jądra do użycia. Domyślnie jest to jądro Gaussa i ma trochę wewnętrznej magii do znalezienia dobrej przepustowości. density()Zwraca obiekt o nazwie komponentu y, który ma wygładzonej gęstości przy różnych wartościach x. Drugi wiersz kodu powyżej znajduje tylko xwartość, gdzie yjest maksymalna.