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))