Estymator gęstości jądra (KDE) tworzy rozkład będący mieszanką lokalizacji rozkładu jądra, więc aby wyciągnąć wartość z oszacowania gęstości jądra, wystarczy (1) narysować wartość z gęstości jądra, a następnie (2) niezależnie wybierz jeden z punktów danych losowo i dodaj jego wartość do wyniku (1).
Oto wynik tej procedury zastosowanej do zestawu danych takiego jak ten w pytaniu.
Histogram po lewej stronie przedstawia próbkę. Dla porównania, czarna krzywa przedstawia gęstość, z której narysowano próbkę. Czerwona krzywa przedstawia KDE próbki (przy użyciu wąskiej przepustowości). (Nie jest problemem, ani nawet nieoczekiwanym, że czerwone piki są krótsze niż czarne piki: KDE rozkłada różne rzeczy, więc piki będą niższe, aby to zrekompensować).
Histogram po prawej przedstawia próbkę (tego samego rozmiaru) z KDE. Czarne i czerwone krzywe są takie same jak poprzednio.
Oczywiście procedura stosowana do pobierania próbek z prac gęstościowych. Jest również niezwykle szybki: R
poniższa implementacja generuje miliony wartości na sekundę z dowolnego KDE. Skomentowałem to mocno, aby pomóc w przenoszeniu do Pythona lub innych języków. Sam algorytm próbkowania jest zaimplementowany w funkcji rdens
z liniami
rkernel <- function(n) rnorm(n, sd=width)
sample(x, n, replace=TRUE) + rkernel(n)
rkernel
rysuje n
próbki IDID z funkcji jądra, podczas gdy sample
rysuje n
próbki z zamianą z danych x
. Operator „+” dodaje dwie tablice próbek element po elemencie.
K.faK.x =( x1, x2), … , Xn)
fax^;K.( x ) = 1n∑i = 1nfaK.( x - xja) .
Xxja1 / njaYX+ YxX
faX+ Y( x )= Pr ( X+ Y≤ x )= ∑i = 1nPr ( X+ Y≤ x ∣ X= xja) Pr ( X= xja)= ∑i = 1nPr ( xja+ Y≤ x ) 1n= 1n∑i = 1nPr ( Y≤ x - xja)= 1n∑i = 1nfaK.( x - xja)= F.x^;K.( X ) ,
jak twierdzono.
#
# Define a function to sample from the density.
# This one implements only a Gaussian kernel.
#
rdens <- function(n, density=z, data=x, kernel="gaussian") {
width <- z$bw # Kernel width
rkernel <- function(n) rnorm(n, sd=width) # Kernel sampler
sample(x, n, replace=TRUE) + rkernel(n) # Here's the entire algorithm
}
#
# Create data.
# `dx` is the density function, used later for plotting.
#
n <- 100
set.seed(17)
x <- c(rnorm(n), rnorm(n, 4, 1/4), rnorm(n, 8, 1/4))
dx <- function(x) (dnorm(x) + dnorm(x, 4, 1/4) + dnorm(x, 8, 1/4))/3
#
# Compute a kernel density estimate.
# It returns a kernel width in $bw as well as $x and $y vectors for plotting.
#
z <- density(x, bw=0.15, kernel="gaussian")
#
# Sample from the KDE.
#
system.time(y <- rdens(3*n, z, x)) # Millions per second
#
# Plot the sample.
#
h.density <- hist(y, breaks=60, plot=FALSE)
#
# Plot the KDE for comparison.
#
h.sample <- hist(x, breaks=h.density$breaks, plot=FALSE)
#
# Display the plots side by side.
#
histograms <- list(Sample=h.sample, Density=h.density)
y.max <- max(h.density$density) * 1.25
par(mfrow=c(1,2))
for (s in names(histograms)) {
h <- histograms[[s]]
plot(h, freq=FALSE, ylim=c(0, y.max), col="#f0f0f0", border="Gray",
main=paste("Histogram of", s))
curve(dx(x), add=TRUE, col="Black", lwd=2, n=501) # Underlying distribution
lines(z$x, z$y, col="Red", lwd=2) # KDE of data
}
par(mfrow=c(1,1))