Надежность режима из образца MCMC

В своей книге «Анализ байесовских данных» Джон Крушке утверждает, что при использовании JAGS из R

... оценка режима из выборки MCMC может быть довольно нестабильной, поскольку оценка основана на алгоритме сглаживания, который может быть чувствителен к случайным ударам и пульсациям в выборке MCMC. ( Выполнение байесовского анализа данных , стр. 205, раздел 8.2.5.1)

Хотя у меня есть представление об алгоритме Метрополиса и точных формах, таких как выборка Гиббса, я не знаком с алгоритмом сглаживания, который также упоминается, и почему он будет означать, что оценка режима из выборки MCMC является нестабильной. Может ли кто-нибудь дать интуитивное представление о том, что делает алгоритм сглаживания и почему он делает оценку режима нестабильной?

У меня нет книги под рукой, поэтому я не уверен, какой метод сглаживания использует Крушке, но для интуиции рассмотрим этот график из 100 выборок из стандартной нормали вместе с оценками плотности ядра по Гауссу с использованием различных полос пропускания от 0,1 до 1,0. (Вкратце, гауссовы KDE - это своего рода сглаженная гистограмма: они оценивают плотность, добавляя гауссиан для каждой точки данных со средним значением при наблюдаемом значении.)

Вы можете видеть, что даже когда сглаживание создает унимодальное распределение, режим обычно ниже известного значения 0.

Далее приведен график предполагаемой моды (ось Y) по ширине полосы ядра, используемой для оценки плотности с использованием той же выборки. Надеемся, что это дает некоторую интуицию к тому, как оценка изменяется с параметрами сглаживания.

#!/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))

Шон Пасхальный дал хороший ответ; вот как это на самом деле делается с помощью сценариев R, которые поставляются вместе с книгой Крушке. plotPost()Функция определена в R сценария с именем DBDA2E-utilities.R. Отображает приблизительный режим. Внутри определения функции есть две строки:

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

density()Функция поставляется с базовой STATs пакетом R, и реализует фильтровального ядро плотности рода Шон Easter описал. Он имеет необязательные аргументы для пропускной способности сглаживающего ядра и типа используемого ядра. По умолчанию используется ядро ​​Гаусса, и есть некоторая внутренняя магия для нахождения хорошей пропускной способности. density()Функция возвращает объект с компонентом с именем , yкоторый имеет сглаженные плотности при различных значениях x. Вторая строка кода, приведенная выше, просто находит xзначение, где yэто максимум.