Jak losowo narysować wartość z oszacowania gęstości jądra?


10

Mam pewne spostrzeżenia i chcę naśladować próbkowanie na podstawie tych spostrzeżeń. Rozważam tutaj model nieparametryczny, w szczególności używam wygładzania jądra, aby oszacować CDF na podstawie ograniczonych obserwacji, a następnie losowo dobieram wartości z uzyskanego CDF. Oto mój kod (chodzi o to, aby losowo uzyskać kumulatywny prawdopodobieństwo za pomocą równomiernego rozkładu i weź odwrotność CDF w stosunku do wartości prawdopodobieństwa)

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)

Jak pokazano w kodzie, użyłem syntetycznego przykładu do przetestowania mojej procedury, ale wynik jest niezadowalający, co ilustrują dwie liczby poniżej (pierwsza dotyczy symulacji obserwacji, a druga rysunek pokazuje histogram sporządzony z szacowanego CDF) :

Rycina 1 Rysunek 2

Czy jest ktoś, kto wie, gdzie jest problem? Z góry dziękuję.


Odwrotne próbkowanie transformaty zależy od zastosowania odwrotnego CDF. en.wikipedia.org/wiki/Inverse_transform_sampling
Sycorax mówi Przywróć Monikę

1
Twój estymator gęstości jądra tworzy rozkład będący mieszanką lokalizacji rozkładu jądra, więc wszystko, czego potrzebujesz, aby wyciągnąć wartość z oszacowania gęstości jądra, to (1) narysuj wartość z gęstości jądra, a następnie (2) niezależnie wybierz jedno z dane wskazują losowo i dodają swoją wartość do wyniku (1). Próba bezpośredniego odwrócenia KDE będzie znacznie mniej wydajna.
whuber

@Sycorax Ale rzeczywiście postępuję zgodnie z procedurą próbkowania odwrotnej transformacji opisaną w Wiki. Zobacz kod: p = rand; [~, idx] = sort (abs (cdf (:, 2) - p)); rndval (i, 1) = cdf (idx (1), 1);
emberbillow

@ whuber Nie jestem pewien, czy rozumiem twój pomysł, czy nie. Pomóż sprawdzić: najpierw ponownie próbkuj wartość z obserwacji; a następnie narysuj wartość z jądra, powiedzmy standardowy rozkład normalny; wreszcie dodać je razem?
emberbillow

Odpowiedzi:


12

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.

Postać

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: Rponiż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 rdensz liniami

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

rkernelrysuje npróbki IDID z funkcji jądra, podczas gdy samplerysuje npróbki z zamianą z danych x. Operator „+” dodaje dwie tablice próbek element po elemencie.


K.faK.x=(x1,x2),,xn)

fax^;K.(x)=1nja=1nfaK.(x-xja).

Xxja1/njaYX+YxX

faX+Y(x)=Par(X+Yx)=ja=1nPar(X+YxX=xja)Par(X=xja)=ja=1nPar(xja+Yx)1n=1nja=1nPar(Yx-xja)=1nja=1nfaK.(x-xja)=fax^;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))

Cześć @ Whuber, chcę zacytować ten pomysł w mojej pracy. Czy masz jakieś artykuły na ten temat? Dziękuję Ci.
emberbillow

2

Próbkujesz najpierw z CDF, odwracając go. Odwrotny CDF nazywany jest funkcją kwantylową; jest to mapowanie z [0,1] na domenę RV. Następnie próbkuje się losowe jednolite RV jako percentyle i przekazuje się je do funkcji kwantylu, aby uzyskać losową próbkę z tego rozkładu.


2
To jest trudna droga: patrz mój komentarz do pytania.
whuber

2
@ whuber dobry punkt. Bez zbytniego zaangażowania w aspekty programistyczne, zakładałem, że w tym przypadku musimy współpracować z CDF. Bez wątpienia elementy wewnętrzne takiej funkcji przyjmują wygładzoną gęstość jądra, a następnie integrują ją w celu uzyskania CDF. W tym momencie prawdopodobnie lepiej i szybciej jest używać odwrotnego próbkowania transformacji. Jednak twoja sugestia, aby po prostu użyć gęstości i próbki prosto z mieszaniny jest lepsza.
AdamO

@AdamO Dziękujemy za odpowiedź. Ale mój kod rzeczywiście podąża za tym samym pomysłem, co powiedziałeś tutaj. Nie wiem, dlaczego nie można odtworzyć wzorów trójmodalnych.
emberbillow

@AdamO Tutaj, czy słowo „internals” w twoim komentarzu powinno być „interwałami”? Dziękuję Ci.
emberbillow

Ember, „wewnętrzne” ma dla mnie idealny sens. Taka funkcja musi zintegrować gęstość mieszaniny i zbudować odwrotność: jest to bałaganiarski, skomplikowany numerycznie proces, jak sugeruje AdamO, i dlatego zostałby pochowany w obrębie funkcji - jej „elementów wewnętrznych”.
whuber

1

W tym miejscu chcę również opublikować kod Matlab zgodnie z ideą opisaną przez Whucera, aby pomóc tym, którzy znają Matlaba bardziej niż 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;

Oto wynik: wyniki

Poinformuj mnie, jeśli ktoś znajdzie jakiś problem z moim zrozumieniem i kodem. Dziękuję Ci.


1
Dodatkowo stwierdziłem, że mój kod w pytaniu jest prawidłowy. Obserwacja, że ​​wzorca nie można odtworzyć, wynika w dużej mierze z wyboru szerokości pasma.
emberbillow

0

Nie przyglądając się zbytnio swojej implementacji, nie otrzymuję w pełni twojej procedury indeksowania, aby czerpać z ICDF. Myślę, że czerpiesz z CDF, a nie odwrotnie. Oto moja implementacja:

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

2
Jeśli masz cdf F, prawdą jest, że F (X) jest jednolity. Otrzymujesz X, biorąc odwrotną liczbę cdf liczby losowej z jednolitego rozkładu. Myślę, że problem polega na tym, jak określić odwrotność podczas tworzenia gęstości jądra.
Michael R. Chernick

Dziękuję za Twoją odpowiedź. Nie próbowałem bezpośrednio z CDF. Kod pokazuje, że rzeczywiście zrobiłem to samo, co próbkowanie z transformacją odwrotną. p = rand; % ta linia otrzymuje jednolitą liczbę losową jako skumulowane prawdopodobieństwo. [~, idx] = sort (abs (cdf (:, 2) - p)); rndval (i, 1) = cdf (idx (1), 1);% te dwie linie mają określić kwantyl odpowiadający skumulowanemu prawdopodobieństwu
emberbillow
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.