Zuverlässigkeit des Modus anhand eines MCMC-Beispiels

12

John Kruschke erklärt in seinem Buch Doing Bayesian Data Analysis, dass bei der Verwendung von JAGS von R

... die Schätzung des Modus aus einem MCMC-Sample kann ziemlich instabil sein, da die Schätzung auf einem Glättungsalgorithmus basiert, der für zufällige Unebenheiten und Welligkeiten im MCMC-Sample empfindlich sein kann. ( Bayesianische Datenanalyse durchführen , Seite 205, Abschnitt 8.2.5.1)

Obwohl ich den Metropolis-Algorithmus und genaue Formen wie Gibbs-Sampling kenne, kenne ich den ebenfalls angedeuteten Glättungsalgorithmus nicht und weiß nicht, warum die Schätzung des Modus aus dem MCMC-Sample instabil ist. Kann jemand einen intuitiven Einblick geben, was der Glättungsalgorithmus tut und warum er die Schätzung des Modus instabil macht?

Morgan Ball
quelle
2
John Kruschke spricht über einen Algorithmus zur Modenschätzung, der auf der Schätzung der Kerneldichte basiert.
Andrey Kolyadin
2
Dieser Link kann hilfreich sein.
Andrey Kolyadin
Sofern ich mich in diesem Bereich der Statistik nicht irre, gibt JAGS eine Reihe von Stichproben aus der posterioren Verteilung und nicht aus einer Wahrscheinlichkeitsdichtefunktion aus, sodass nicht sicher ist, ob eine Kerndichteschätzung vorliegt. Vielen Dank für den Link.
Morgan Ball
Ich denke, das hat wahrscheinlich mehr damit zu tun, wie Sie einen Modus aus einem großen Sample einer stetigen Variablen erhalten, bei dem möglicherweise nicht mehr als einer von einem bestimmten Wert vorhanden ist und Sie das Sample gruppieren (oder glätten) müssen.
Morgan Ball
1
Sie können den Modus als Wert mit maximaler Dichte bei der Schätzung der Kerneldichte erhalten. (Zumindest ist es das, was ich tue, und wenn ich mich nicht irre, benutze J. Kruschke in seinen Beispielen den gleichen Ansatz.)
Andrey Kolyadin

Antworten:

8

Ich habe das Buch nicht zur Hand, daher bin ich mir nicht sicher, welche Glättungsmethode Kruschke verwendet, aber für die Intuition betrachten Sie diese Darstellung von 100 Proben aus einer Standardnormalen zusammen mit Schätzungen der Gaußschen Kerneldichte unter Verwendung verschiedener Bandbreiten von 0,1 bis 1,0. (Kurz gesagt, Gaußsche KDEs sind eine Art geglättetes Histogramm: Sie schätzen die Dichte, indem sie für jeden Datenpunkt einen Gaußschen Wert addieren, dessen Mittelwert beim beobachteten Wert liegt.)

Sie können sehen, dass der Modus im Allgemeinen unter dem bekannten Wert von 0 liegt, selbst wenn durch das Glätten eine unimodale Verteilung erstellt wird.

Bildbeschreibung hier eingeben

Weiter ist hier eine grafische Darstellung des geschätzten Modus (y-Achse) nach Kernbandbreite, der zum Schätzen der Dichte unter Verwendung derselben Stichprobe verwendet wird. Hoffentlich gibt dies einen Anhaltspunkt dafür, wie sich die Schätzung mit den Glättungsparametern ändert.

Bildbeschreibung hier eingeben

#!/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
quelle
5

Sean Easter gab eine nette Antwort; Das machen die R-Skripte in Kruschkes Buch so. Die plotPost()Funktion ist im genannten R-Skript definiert DBDA2E-utilities.R. Es zeigt den geschätzten Modus an. Innerhalb der Funktionsdefinition gibt es diese zwei Zeilen:

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

Die density()Funktion wird mit dem Basisstatistikpaket von R ausgeliefert und implementiert einen Kernel-Dichtefilter der beschriebenen Art Sean Easter. Es gibt optionale Argumente für die Bandbreite des Glättungskerns und für die Art des zu verwendenden Kernels. Der Standardwert ist ein Gaußscher Kernel und es gibt eine gewisse Magie, eine gute Bandbreite zu finden. Die density()Funktion gibt ein Objekt mit einer benannten Komponente zurück y, die die geglätteten Dichten bei verschiedenen Werten aufweist x. In der zweiten Codezeile oben wird nur der xWert gefunden, bei dem ydas Maximum erreicht ist.

John K. Kruschke
quelle