Wiarygodność trybu z próbki MCMC


12

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?


2
Mówię o Johnie Kruschke mówiącym o algorytmie szacowania trybu opartym na szacowaniu gęstości jądra.
Andrey Kolyadin

2
Ten link może być pomocny.
Andrey Kolyadin

O ile nie mylę się, że jestem nowy w tym obszarze statystyki, JAGS generuje zestaw próbek z rozkładu tylnego zamiast funkcji gęstości prawdopodobieństwa, więc nie jestem pewien, czy w to wchodzi oszacowanie gęstości jądra. Dzięki za link.
Morgan Ball

Myślę, że jest to prawdopodobnie bardziej związane z tym, jak uzyskać tryb z dużej próbki zmiennej ciągłej, w której może nie być więcej niż jedna z jakiejkolwiek konkretnej wartości, więc musisz zgrupować (lub wygładzić) próbkę.
Morgan Ball

1
Możesz uzyskać tryb jako wartość z maksymalną gęstością przy szacowaniu gęstości jądra. (Przynajmniej tak robię i jeśli się nie mylę, J. Kruschke stosuje takie samo podejście w swoich przykładach)
Andrey Kolyadin

Odpowiedzi:


8

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.

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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

5

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.

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.