Wie kann ich einen Wert zufällig aus einer Kernel-Dichteschätzung ziehen?

Ich habe einige Beobachtungen und möchte anhand dieser Beobachtungen Stichproben nachahmen. Hier betrachte ich ein nicht parametrisches Modell, insbesondere verwende ich die Kernel-Glättung, um eine CDF aus den begrenzten Beobachtungen zu schätzen. Dann ziehe ich zufällige Werte aus der erhaltenen CDF. Das Folgende ist mein Code (die Idee ist, zufällig eine kumulative zu erhalten Wahrscheinlichkeit unter Verwendung einer gleichmäßigen Verteilung und nehmen die Umkehrung der CDF in Bezug auf den Wahrscheinlichkeitswert)

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)

Wie im Code gezeigt, habe ich ein synthetisches Beispiel verwendet, um mein Verfahren zu testen, aber das Ergebnis ist unbefriedigend, wie die beiden folgenden Abbildungen zeigen (die erste ist für die simulierten Beobachtungen und die zweite Abbildung zeigt das Histogramm, das aus der geschätzten CDF gezogen wurde). ::

Gibt es jemanden, der weiß, wo das Problem liegt? Vielen Dank im Voraus.

Ein Kernel-Dichteschätzer (KDE) erzeugt eine Verteilung, die eine Ortsmischung der Kernel-Verteilung ist. Um also einen Wert aus der Kernel-Dichteschätzung zu ziehen, müssen Sie lediglich (1) einen Wert aus der Kernel-Dichte ziehen und dann (2) Wählen Sie einen der Datenpunkte unabhängig voneinander nach dem Zufallsprinzip aus und addieren Sie seinen Wert zum Ergebnis von (1).

Hier ist das Ergebnis dieser Prozedur, die auf einen Datensatz wie den in der Frage angewendeten angewendet wird.

Das Histogramm links zeigt die Probe. Als Referenz zeigt die schwarze Kurve die Dichte, aus der die Probe gezogen wurde. Die rote Kurve zeigt die KDE der Probe (unter Verwendung einer schmalen Bandbreite). (Es ist kein Problem oder sogar unerwartet, dass die roten Peaks kürzer als die schwarzen Peaks sind: Der KDE verteilt die Dinge, sodass die Peaks zum Ausgleich niedriger werden.)

Das Histogramm rechts zeigt eine Probe (gleicher Größe) aus dem KDE. Die schwarzen und roten Kurven sind die gleichen wie zuvor.

Offensichtlich funktioniert das Verfahren zur Probenahme aus der Dichte. Es ist auch extrem schnell: Die folgende RImplementierung generiert Millionen von Werten pro Sekunde aus jedem KDE. Ich habe es stark kommentiert, um die Portierung auf Python oder andere Sprachen zu unterstützen. Der Abtastalgorithmus selbst ist in der Funktion rdensmit den Linien implementiert

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

rkernelzieht nIId Proben aus der Kernel - Funktion während sampleziehen nProben mit Ersatz aus den Daten x. Der Operator "+" fügt die beiden Arrays von Samples Komponente für Komponente hinzu.

Für diejenigen, die eine formelle Demonstration der Korrektheit wünschen, biete ich sie hier an. Es sei die Kernelverteilung mit CDF und die Daten seien . Nach der Definition einer Kernelschätzung ist die CDF der KDEF K x = ( x 1 , x 2 , … , x n$K$$F_K$$\mathbf{x}=(x_1, x_2, \ldots, x_n)$

Das vorhergehende Rezept besagt, dass aus der empirischen Verteilung der Daten gezogen werden soll (dh es wird der Wert mit einer Wahrscheinlichkeit von für jedes ), unabhängig eine Zufallsvariable aus der Kernverteilung gezogen und summiert werden soll. Ich schulde Ihnen einen Beweis dafür, dass die Verteilungsfunktion von die des KDE ist. Beginnen wir mit der Definition und sehen, wohin sie führt. Sei eine beliebige reelle Zahl. Konditionierung auf gibtx i 1 / n i Y X + Y x $X$$x_i$$1/n$$i$$Y$$X+Y$$x$$X$

wie behauptet.

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

Sie probieren zuerst aus der CDF, indem Sie sie invertieren. Die inverse CDF wird als Quantilfunktion bezeichnet; Es ist eine Zuordnung von [0,1] zur Domäne des Wohnmobils. Anschließend werden zufällige einheitliche RVs als Perzentile abgetastet und an die Quantilfunktion übergeben, um eine zufällige Stichprobe aus dieser Verteilung zu erhalten.

Hier möchte ich auch den Matlab-Code veröffentlichen, der der von whuber beschriebenen Idee folgt, um denjenigen zu helfen, die mit Matlab besser vertraut sind als 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;

Folgendes ist das Ergebnis:

Bitte sagen Sie mir, wenn jemand ein Problem mit meinem Verständnis und dem Code findet. Vielen Dank.

Ohne Ihre Implementierung zu genau zu betrachten, kann ich Ihr Indizierungsverfahren nicht vollständig aus dem ICDF ziehen. Ich denke, Sie ziehen aus der CDF, nicht es ist umgekehrt. Hier ist meine Implementierung:

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