¿Cómo puedo dibujar un valor al azar de una estimación de densidad del núcleo?

Tengo algunas observaciones y quiero imitar el muestreo basado en estas observaciones. Aquí considero un modelo no paramétrico, específicamente, utilizo el suavizado del núcleo para estimar un CDF a partir de las observaciones limitadas. Luego dibujo valores al azar del CDF obtenido. El siguiente es mi código, (la idea es obtener aleatoriamente un resultado acumulativo probabilidad usando distribución uniforme, y tome el inverso de la CDF con respecto al valor de probabilidad)

x = [randn(100, 1); rand(100, 1)+4; rand(100, 1)+8];
[f, xi] = ksdensity(x, 'Function', 'cdf', 'NUmPoints', 300);
cdf = [xi', f'];
nbsamp = 100;
rndval = zeros(nbsamp, 1);
for i = 1:nbsamp
    p = rand;
   [~, idx] = sort(abs(cdf(:, 2) - p));
   rndval(i, 1) = cdf(idx(1), 1);
end
figure(1);
hist(x, 40)
figure(2);
hist(rndval, 40)

Como se muestra en el código, utilicé un ejemplo sintético para probar mi procedimiento, pero el resultado es insatisfactorio, como se ilustra en las dos figuras a continuación (la primera es para las observaciones simuladas y la segunda muestra el histograma extraído del CDF estimado) :

¿Hay alguien que sepa dónde está el problema? Gracias de antemano.

Un estimador de densidad de kernel (KDE) produce una distribución que es una mezcla de ubicación de la distribución de kernel, por lo que para dibujar un valor de la estimación de densidad de kernel, todo lo que necesita hacer es (1) dibujar un valor de la densidad de kernel y luego (2) seleccione independientemente uno de los puntos de datos al azar y agregue su valor al resultado de (1).

Aquí está el resultado de este procedimiento aplicado a un conjunto de datos como el de la pregunta.

El histograma a la izquierda representa la muestra. Como referencia, la curva negra traza la densidad de la que se extrajo la muestra. La curva roja traza el KDE de la muestra (usando un ancho de banda estrecho). (No es un problema, ni siquiera inesperado, que los picos rojos sean más cortos que los negros: el KDE extiende las cosas, por lo que los picos se reducirán para compensar).

El histograma de la derecha muestra una muestra (del mismo tamaño) del KDE. Las curvas negras y rojas son las mismas que antes.

Evidentemente, el procedimiento utilizado para tomar muestras de la densidad funciona. También es extremadamente rápido: la Rimplementación a continuación genera millones de valores por segundo desde cualquier KDE. Lo he comentado mucho para ayudar a portar a Python u otros idiomas. El algoritmo de muestreo en sí se implementa en la función rdenscon las líneas.

rkernel <- function(n) rnorm(n, sd=width) 
sample(x, n, replace=TRUE) + rkernel(n)  

rkerneldibuja nmuestras iid de la función del kernel mientras sampledibuja nmuestras con reemplazo de los datos x. El operador "+" agrega las dos matrices de muestras componente por componente.

Para aquellos que desean una demostración formal de corrección, la ofrezco aquí. Deje que represente la distribución del núcleo con CDF y deje que los datos sean . Según la definición de una estimación del núcleo, el CDF de KDE esF K x = ( x 1 , x 2 , … , x n$K$$F_K$$\mathbf{x}=(x_1, x_2, \ldots, x_n)$

La receta anterior dice que se extrae de la distribución empírica de los datos (es decir, se alcanza el valor con probabilidad para cada ), se extrae independientemente una variable aleatoria de la distribución del núcleo y se suma. Te debo una prueba de que la función de distribución de es la de KDE. Comencemos con la definición y veamos a dónde lleva. Deje ser cualquier número real. Condicionamiento en dax i 1 / n i Y X + Y x $X$$x_i$$1/n$$i$$Y$$X+Y$$x$$X$

según lo reclamado.

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

Primero toma muestras del CDF invirtiéndolo. El CDF inverso se llama función cuantil; es un mapeo desde [0,1] al dominio de la RV. Luego se muestrean RVs aleatorias uniformes como percentiles y se pasan a la función cuantil para obtener una muestra aleatoria de esa distribución.

Aquí, también quiero publicar el código de Matlab siguiendo la idea descrita por whuber, para ayudar a aquellos que están más familiarizados con Matlab que con R.

x = exprnd(3, [300, 1]);
[~, ~, bw] = ksdensity(x, 'kernel', 'normal', 'NUmPoints', 800);

k = 0.25; % control the uncertainty of generated values, the larger the k the greater the uncertainty
mstd = bw*k;
rkernel = mstd*randn(300, 1);
sampleobs = randsample(x, 300, true);
simobs = sampleobs(:) + rkernel(:);

figure(1);
subplot(1,2,1);
hist(x, 50);title('Original sample');
subplot(1,2,2);
hist(simobs, 50);title('Simulated sample');
axis tight;

El siguiente es el resultado:

Por favor, dígame si alguien encuentra algún problema con mi comprensión y el código. Gracias.

Sin mirar demasiado de cerca su implementación, no entiendo completamente su procedimiento de indexación para extraer del ICDF. Creo que sacas del CDF, no es inverso. Aquí está mi implementación:

import sys
sys.path.insert(0, './../../../Python/helpers')
import numpy as np
import scipy.stats as stats
from sklearn.neighbors import KernelDensity

def rugplot(axis,x,color='b',label='draws',shape='+',alpha=1):
    axis.plot(x,np.ones(x.shape)*0,'b'+shape,ms=20,label=label,c=color,alpha=alpha);
    #axis.set_ylim([0,max(axis.get_ylim())])

def PDF(x):
    return 0.5*(stats.norm.pdf(x,loc=6,scale=1)+ stats.norm.pdf(x,loc=18,scale=1));

def CDF(x,PDF):
    temp = np.linspace(-10,x,100)
    pdf = PDF(temp);
    return np.trapz(pdf,temp);

def iCDF(p,x,cdf):
    return np.interp(p,cdf,x);

res = 1000;
X = np.linspace(0,24,res);
P = np.linspace(0,1,res)
pdf  = np.array([PDF(x) for x in X]);#attention dont do [ for x in x] because it overrides original x value
cdf  = np.array([CDF(x,PDF) for x in X]);
icdf = [iCDF(p,X,cdf) for p in P];

#draw pdf and cdf
f,(ax1,ax2) = plt.subplots(1,2,figsize=(18,4.5));
ax1.plot(X,pdf, '.-',label = 'pdf');
ax1.plot(X,cdf, '.-',label = 'cdf');
ax1.legend();
ax1.set_title('PDF & CDF')

#draw inverse cdf
ax2.plot(cdf,X,'.-',label  = 'inverse by swapping axis');
ax2.plot(P,icdf,'.-',label = 'inverse computed');
ax2.legend();
ax2.set_title('inverse CDF');

#draw from custom distribution
N = 100;
p_uniform = np.random.uniform(size=N)
x_data  = np.array([iCDF(p,X,cdf) for p in p_uniform]);

#visualize draws
a = plt.figure(figsize=(20,8)).gca();
rugplot(a,x_data);

#histogram
h = np.histogram(x_data,bins=24);
a.hist(x_data,bins=h[1],alpha=0.5,normed=True);