Keandalan Mode dari sampel MCMC

Dalam bukunya Doing Bayesian Data Analysis, John Kruschke menyatakan bahwa dalam menggunakan JAGS dari R

... estimasi mode dari sampel MCMC bisa agak tidak stabil karena estimasi tersebut didasarkan pada algoritma perataan yang bisa peka terhadap benjolan dan riak acak dalam sampel MCMC. ( Melakukan Analisis Data Bayesian , halaman 205, bagian 8.2.5.1)

Meskipun saya memiliki pemahaman tentang algoritma Metropolis dan bentuk-bentuk yang tepat seperti Gibbs sampling, saya tidak akrab dengan algoritma smoothing yang disinggung juga dan mengapa itu berarti estimasi mode dari sampel MCMC tidak stabil. Adakah yang bisa memberikan wawasan intuitif tentang apa yang dilakukan algoritma penghalusan dan mengapa hal itu membuat estimasi mode tidak stabil?

Saya tidak memiliki buku yang ada di tangan jadi saya tidak yakin apa metode penghalusan yang digunakan Kruschke, tetapi untuk intuisi, pertimbangkan plot 100 sampel ini dari standar normal, bersama dengan perkiraan kepadatan kernel Gaussian menggunakan berbagai bandwidth dari 0,1 ke 1,0. (Secara singkat, KDE Gaussian adalah semacam histogram yang dihaluskan: Mereka memperkirakan kepadatan dengan menambahkan Gaussian untuk setiap titik data, dengan rata-rata pada nilai yang diamati.)

Anda dapat melihat bahwa bahkan setelah penghalusan menghasilkan distribusi unimodal, mode umumnya di bawah nilai 0 yang diketahui.

Terlebih lagi, inilah plot dari mode yang diperkirakan (sumbu-y) oleh bandwidth kernel yang digunakan untuk memperkirakan kepadatan, menggunakan sampel yang sama. Semoga ini memberikan beberapa intuisi tentang bagaimana estimasi bervariasi dengan parameter smoothing.

#!/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 memberikan jawaban yang bagus; inilah bagaimana ini sebenarnya dilakukan oleh skrip R yang datang dengan buku Kruschke. The plotPost()Fungsi didefinisikan dalam naskah R bernama DBDA2E-utilities.R. Ini menampilkan mode perkiraan. Di dalam definisi fungsi, ada dua baris ini:

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

The density()Fungsi dilengkapi dengan paket dasar statistik dari R, dan alat kernel filter density dari jenis Sean Paskah dijelaskan. Ini memiliki argumen opsional untuk bandwidth kernel smoothing, dan untuk jenis kernel yang digunakan. Ini default ke kernel Gaussian, dan memiliki beberapa keajaiban internal untuk menemukan bandwidth yang bagus. The density()fungsi mengembalikan sebuah objek dengan komponen bernama yyang memiliki kepadatan merapikan di berbagai nilai x. Baris kode kedua, di atas, hanya menemukan xnilai di mana ymaksimum.