Un estimador de densidad de kernel (KDE), que produce una distribución que es una ubicación de la mezcla de la distribución del núcleo, por lo que para dibujar un valor de la estimación de densidad de kernel todo lo que necesitas hacer es (1) la elaboración de un valor de la densidad de kernel y luego (2) de forma independiente seleccione uno de los puntos de datos al azar y añadir su valor para el resultado de (1).
Aquí está el resultado de este procedimiento se aplica a un conjunto de datos como el de la pregunta.
![Figure]()
El histograma de la izquierda representa la muestra. Para su referencia, la curva negra parcelas de la densidad a partir de la cual la muestra fue tomada. La curva roja representa los KDE de la muestra (con un ancho de banda estrecho). (No es un problema, o incluso inesperadas, que los picos rojos son más cortos que los picos de negro: el KDE diferenciales de las cosas, de modo que los picos recibirá inferior para compensar.)
El histograma de la derecha muestra un ejemplo (del mismo tamaño) de KDE. El negro y el rojo curvas son las mismas que antes.
Evidentemente, el procedimiento utilizado para la muestra de la densidad de las obras. También es extremadamente rápido: el R de ejecución a continuación genera millones de valores por segundo desde cualquier KDE. Me han comentado mucho para ayudar en la migración a Python u otros idiomas. El algoritmo de muestreo mismo es implementado en la función rdens con las líneas
rkernel <- function(n) rnorm(n, sd=width)
sample(x, n, replace=TRUE) + rkernel(n)
rkernel dibuja n iid muestras de que el kernel funcione mientras sample dibuja n de muestras con reemplazo a partir de los datos x. El operador"+", agrega las dos matrices de muestras componente por componente.
Para aquellos que quieran una demostración formal de la corrección, me ofrezco aquí. Deje $K$ representan el núcleo de distribución con CDF $F_K$ y dejar que los datos se $\mathbf{x}=(x_1, x_2, \ldots, x_n)$. Por la definición de un núcleo de estimación, la CDF de KDE es
$$F_{\mathbf{\hat{x}};\, K}(x) = \frac{1}{n}\sum_{i=1}^n F_K(x-x_i).$$
The preceding recipe says to draw $X$ from the empirical distribution of the data (that is, it attains the value $x_i$ with probability $1/n$ for each $i$), independently draw a random variable $S$ from the kernel distribution, and sum them. I owe you a proof that the distribution function of $X+Y$ is that of the KDE. Let's start with the definition and see where it leads. Let $x$ be any real number. Conditioning on $X$ gives
$$\eqalign{
F_{X+Y}(x) &= \Pr(X+Y \le x) \\
&= \sum_{i=1}^n \Pr(X+Y \le x \a mediados de X=x_i) \Pr(X=x_i) \\
&= \sum_{i=1}^n \Pr(x_i + Y \le x) \frac{1}{n} \\
&= \frac{1}{n}\sum_{i=1}^n \Pr(Y \le x-x_i) \\
&= \frac{1}{n}\sum_{i=1}^n F_K(x-x_i) \\
&= F_{\mathbf{\hat{x}};\, K}(x),
}$$
como se reivindica.
#
# 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))
