Comment puis-je tirer une valeur au hasard à partir d'une estimation de la densité du noyau?

10

J'ai quelques observations, et je veux imiter l'échantillonnage basé sur ces observations. Ici, je considère un modèle non paramétrique, en particulier, j'utilise le lissage du noyau pour estimer un CDF à partir des observations limitées.Puis je tire des valeurs au hasard à partir du CDF obtenu.Le code suivant est (l'idée est d'obtenir au hasard un cumulatif probabilité en utilisant une distribution uniforme, et prendre l'inverse du CDF par rapport à la valeur de probabilité)

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)

Comme indiqué dans le code, j'ai utilisé un exemple synthétique pour tester ma procédure, mais le résultat n'est pas satisfaisant, comme illustré par les deux figures ci-dessous (la première est pour les observations simulées, et la deuxième figure montre l'histogramme tiré de CDF estimé) :

Figure 1 Figure 2

Y a-t-il quelqu'un qui sait où est le problème? Merci d'avance.

emberbillow
la source
L' échantillonnage par transformée inverse repose sur l'utilisation du CDF inverse . en.wikipedia.org/wiki/Inverse_transform_sampling
Sycorax dit Reinstate Monica
1
Votre estimateur de densité de noyau produit une distribution qui est un mélange d'emplacements de la distribution de noyau, de sorte que tout ce dont vous avez besoin pour tirer une valeur de l'estimation de densité de noyau est (1) dessiner une valeur à partir de la densité de noyau puis (2) sélectionner indépendamment l'un des la donnée pointe au hasard et ajoute sa valeur au résultat de (1). Tenter d'inverser directement KDE sera beaucoup moins efficace.
whuber
@Sycorax Mais je suis en effet la procédure d'échantillonnage à transformée inverse décrite dans Wiki. Veuillez voir le code: p = rand; [~, idx] = sort (abs (cdf (:, 2) - p)); rndval (i, 1) = cdf (idx (1), 1);
emberbillow
@whuber Je ne sais pas si ma compréhension de votre idée est correcte ou non. Veuillez aider à vérifier: commencez par rééchantillonner une valeur à partir des observations; puis tirer une valeur du noyau, disons distribution normale standard; enfin, ajoutez-les ensemble?
emberbillow

Réponses:

12

Un estimateur de densité de noyau (KDE) produit une distribution qui est un mélange d'emplacements de la distribution de noyau, donc pour tirer une valeur de l'estimation de densité de noyau, tout ce que vous devez faire est (1) de tirer une valeur de la densité de noyau, puis (2) sélectionnez indépendamment l'un des points de données au hasard et ajoutez sa valeur au résultat de (1).

Voici le résultat de cette procédure appliquée à un ensemble de données comme celui de la question.

Figure

L'histogramme à gauche représente l'échantillon. Pour référence, la courbe noire représente la densité à partir de laquelle l'échantillon a été prélevé. La courbe rouge trace le KDE de l'échantillon (en utilisant une bande passante étroite). (Ce n'est pas un problème, ni même inattendu, que les pics rouges sont plus courts que les pics noirs: le KDE répartit les choses, donc les pics s'abaissent pour compenser.)

L'histogramme à droite représente un échantillon (de la même taille) du KDE. Les courbes noires et rouges sont les mêmes qu'avant.

De toute évidence, la procédure utilisée pour échantillonner à partir de la densité fonctionne. C'est aussi extrêmement rapide: l' Rimplémentation ci-dessous génère des millions de valeurs par seconde à partir de n'importe quel KDE. Je l'ai fortement commenté pour aider au portage vers Python ou d'autres langages. L'algorithme d'échantillonnage lui-même est implémenté dans la fonction rdensavec les lignes

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

rkerneldessine des néchantillons iid à partir de la fonction noyau tandis que sampledessine des néchantillons avec remplacement à partir des données x. L'opérateur "+" ajoute les deux tableaux d'échantillons composant par composant.


KFKX=(X1,X2,,Xn)

FX^;K(X)=1nje=1nFK(X-Xje).

XXje1/njeOuiX+OuiXX

FX+Oui(X)=Pr(X+OuiX)=je=1nPr(X+OuiXX=Xje)Pr(X=Xje)=je=1nPr(Xje+OuiX)1n=1nje=1nPr(OuiX-Xje)=1nje=1nFK(X-Xje)=FX^;K(X),

comme revendiqué.


#
# 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))
whuber
la source
Salut @whuber, je veux citer cette idée dans mon article. Avez-vous des articles publiés à ce sujet? Je vous remercie.
emberbillow
2

Vous échantillonnez à partir du CDF d'abord en l'inversant. Le CDF inverse est appelé la fonction quantile; c'est une mise en correspondance de [0,1] avec le domaine du RV. Vous échantillonnez ensuite des RV uniformes aléatoires en tant que centiles et les passez à la fonction quantile pour obtenir un échantillon aléatoire à partir de cette distribution.

AdamO
la source
2
C'est la voie difficile: voir mon commentaire sur la question.
whuber
2
@whuber bon point. Sans être trop empêtré dans les aspects programmatiques, je supposais que nous devons travailler avec un CDF dans ce cas. Nul doute que les éléments internes d'une telle fonction prennent une densité lissée du noyau et l' intègrent ensuite pour obtenir un CDF. À ce stade, il est probablement préférable et plus rapide d'utiliser l'échantillonnage par transformée inverse. Cependant, votre suggestion d'utiliser simplement la densité et l'échantillon directement à partir du mélange est meilleure.
AdamO
@AdamO Merci pour votre réponse. Mais mon code suit en effet la même idée que vous avez dit ici. Je ne sais pas pourquoi les motifs tri-modaux ne peuvent pas être reproduits.
emberbillow
@AdamO Voici si le mot "internes" dans votre commentaire doit être "intervalles"? Je vous remercie.
emberbillow
Ember, "internes" est parfaitement logique pour moi. Une telle fonction doit intégrer la densité du mélange et construire un inverse: c'est un processus compliqué et numériquement compliqué comme le laisse entendre AdamO, et serait donc enfoui dans la fonction - ses «internes».
whuber
1

Ici, je veux également publier le code Matlab suivant l'idée décrite par whuber, pour aider ceux qui connaissent mieux Matlab que 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;

Le résultat est le suivant: résultats

Veuillez me dire si quelqu'un trouve un problème avec ma compréhension et le code. Je vous remercie.

emberbillow
la source
1
De plus, j'ai trouvé que mon code dans la question est correct. L'observation selon laquelle le motif ne peut pas être reproduit est largement due au choix de la bande passante.
emberbillow
0

Sans regarder de trop près dans votre implémentation, je n'obtiens pas entièrement votre procédure d'indexation à tirer de l'ICDF. Je pense que vous tirez du CDF, pas son inverse. Voici ma mise en œuvre:

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);
Jan
la source
2
Si vous avez un cdf F, il est vrai que F (X) est uniforme. Vous obtenez donc X en prenant l'inverse cdf d'un nombre aléatoire à partir d'une distribution uniforme. Le problème, je pense, est de savoir comment déterminer l'inverse lorsque vous produisez une densité de noyau.
Michael R. Chernick
Merci pour votre réponse. Je n'ai pas échantillonné directement à partir du CDF. Le code montre que j'ai en effet fait la même chose que l'échantillonnage par transformée inverse. p = rand; % cette ligne obtient un nombre aléatoire uniforme comme probabilité cumulative. [~, idx] = sort (abs (cdf (:, 2) - p)); rndval (i, 1) = cdf (idx (1), 1);% ces deux lignes doivent déterminer le quantile correspondant à la probabilité cumulative
emberbillow