Как я могу получить значение случайно из оценки плотности ядра?

10

У меня есть некоторые наблюдения, и я хочу повторить выборку на основе этих наблюдений. Здесь я рассматриваю непараметрическую модель, в частности, я использую сглаживание ядра для оценки CDF из ограниченных наблюдений. Затем я рисую значения случайным образом из полученного CDF. Ниже приведен мой код (идея состоит в том, чтобы случайным образом получить кумулятивный вероятность с использованием равномерного распределения, и принять обратный CDF по отношению к значению вероятности)

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)

Как показано в коде, я использовал синтетический пример для тестирования моей процедуры, но результат неудовлетворительный, как показано двумя рисунками ниже (первый - для смоделированных наблюдений, а второй - гистограмма, полученная из оценочного CDF) :

фигура 1 фигура 2

Есть кто-нибудь, кто знает, где проблема? Заранее спасибо.

emberbillow
источник
Выборка обратного преобразования зависит от использования обратного CDF. ru.wikipedia.org/wiki/Inverse_transform_sampling
Sycorax сообщает о восстановлении Monica
1
Ваш оценщик плотности ядра создает распределение, которое представляет собой смесь местоположения распределения ядра, поэтому все, что вам нужно, чтобы извлечь значение из оценки плотности ядра, это (1) извлечь значение из плотности ядра, а затем (2) независимо выбрать один из данные указывают случайным образом и добавляют свое значение к результату (1). Попытка инвертировать KDE напрямую будет гораздо менее эффективной.
whuber
@Sycorax Но я действительно следую процедуре выборки обратного преобразования, как описано в Wiki. Пожалуйста, смотрите код: p = rand; [~, idx] = sort (abs (cdf (:, 2) - p)); rndval (i, 1) = cdf (idx (1), 1);
Эмбербиллоу
@whuber Я не уверен, правильно ли я понимаю твою идею или нет. Пожалуйста, помогите проверить: сначала пересчитайте значение из наблюдений; а затем извлечь значение из ядра, скажем, стандартное нормальное распределение; наконец, добавить их вместе?
emberbillow

Ответы:

12

Оценщик плотности ядра (KDE) создает распределение, представляющее собой смесь локаций распределения ядра, поэтому для получения значения из оценки плотности ядра все, что вам нужно сделать, это (1) извлечь значение из плотности ядра, а затем (2) независимо выберите случайным образом одну из точек данных и добавьте ее значение к результату (1).

Вот результат этой процедуры, примененной к набору данных, подобному тому, который указан в вопросе.

фигура

Гистограмма слева изображает образец. Для справки, черная кривая показывает плотность, из которой был взят образец. Красная кривая отображает KDE образца (с использованием узкой полосы пропускания). (Это не проблема или даже неожиданность, что красные пики короче, чем черные пики: KDE распределяет вещи, поэтому пики будут уменьшаться для компенсации.)

Гистограмма справа показывает образец (того же размера) из KDE. Черные и красные кривые такие же, как и раньше.

Очевидно, процедура, используемая для отбора проб из плотины работает. Это также чрезвычайно быстро: приведенная Rниже реализация генерирует миллионы значений в секунду из любого KDE. Я прокомментировал это сильно, чтобы помочь в портировании на Python или другие языки. Сам алгоритм выборки реализован в функции rdensсо строками

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

rkernelрисует nобразцы из функции ядра, а sampleрисует nобразцы с заменой из данных x. Оператор «+» добавляет два массива выборок компонент за компонентом.


КFКИксзнак равно(Икс1,Икс2,...,ИксN)

FИкс^;К(Икс)знак равно1NΣязнак равно1NFК(Икс-Икся),

ИксИкся1/NяYИкс+YИксИкс

FИкс+Y(Икс)знак равноPr(Икс+YИкс)знак равноΣязнак равно1NPr(Икс+YИкс|Иксзнак равноИкся)Pr(Иксзнак равноИкся)знак равноΣязнак равно1NPr(Икся+YИкс)1Nзнак равно1NΣязнак равно1NPr(YИкс-Икся)знак равно1NΣязнак равно1NFК(Икс-Икся)знак равноFИкс^;К(Икс),

как утверждено.


#
# 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
источник
Привет @ whuber, я хочу процитировать эту идею в своей статье. У вас есть какие-то статьи, которые были опубликованы для этого? Спасибо.
Эмбербиллоу
2

Вы сначала делаете выборку из CDF, переворачивая ее. Обратный CDF называется функцией квантиля; это отображение из [0,1] в область RV. Затем вы выбираете случайные однородные RV в виде процентилей и передаете их в функцию квантиля, чтобы получить случайную выборку из этого распределения.

Adamo
источник
2
Это трудный путь: см. Мой комментарий к вопросу.
whuber
2
@ Хороший вопрос. Не слишком погружаясь в программные аспекты, я предполагал, что в этом случае мы должны работать с CDF. Без сомнения, внутреннее устройство такой функции принимает сглаженную плотность ядра и затем интегрирует ее для получения CDF. В этот момент, вероятно, лучше и быстрее использовать выборку с обратным преобразованием. Тем не менее, ваше предложение просто использовать плотность и образец прямо из смеси лучше.
AdamO
@AdamO Спасибо за ваш ответ. Но мой код действительно следует той же идее, что и здесь. Я не знаю, почему тримодальные модели не могут быть воспроизведены.
Эмбербиллоу
@AdamO Вот должно ли слово «внутренности» в вашем комментарии быть «интервалами»? Спасибо.
emberbillow
Эмбер, "внутренности" имеют для меня смысл. Такая функция должна интегрировать плотность смеси и построить обратное: это сложный, численно сложный процесс, как подсказывает AdamO, и поэтому он будет скрыт внутри функции - ее «внутренних органов».
whuber
1

Здесь я также хочу опубликовать код Matlab, следуя идее, описанной whuber, чтобы помочь тем, кто лучше знаком с Matlab, чем с 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;

Ниже приводится результат: Результаты

Пожалуйста, скажите мне, если кто-нибудь обнаружит какие-либо проблемы с моим пониманием и кодом. Спасибо.

emberbillow
источник
1
Кроме того, я обнаружил, что мой код в вопросе является правильным. Замечание о том, что рисунок не может быть воспроизведен, во многом связано с выбором полосы пропускания.
emberbillow
0

Не смотря слишком близко к вашей реализации, я не в полной мере понимаю вашу процедуру индексации из ICDF. Я думаю, что вы черпаете из CDF, а не наоборот. Вот моя реализация:

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
Если у вас есть cdf F, то это правда, что F (X) равномерно. Таким образом, вы получите X, взяв обратный cdf случайного числа из равномерного распределения. Я думаю, что проблема заключается в том, как определить обратное значение при создании плотности ядра.
Майкл Р. Черник
Спасибо за ваш ответ. Я не пробовал напрямую с CDF. Код показывает, что я действительно сделал то же самое, что и выборка с обратным преобразованием. р = ранд; % эта строка получает равномерное случайное число в качестве совокупной вероятности. [~, idx] = sort (abs (cdf (:, 2) - p)); rndval (i, 1) = cdf (idx (1), 1);% эти две строки должны определить квантиль, соответствующую кумулятивной вероятности
emberbillow