Como posso desenhar um valor aleatoriamente a partir de uma estimativa de densidade do kernel?

10

Tenho algumas observações e quero imitar a amostragem com base nessas observações. A partir do momento em que o CDF é capturado, o CDF é transferido para o CDF e o CDF é transferido para o CDF, o que significa que o CDF pode ser usado como um arquivo não-paramétrico, para facilitar a estimativa de um CDF a partir de observações limitadas. probabilidade usando distribuição uniforme e tomar o inverso do CDF com relação ao valor da probabilidade)

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 mostrado no código, usei um exemplo sintético para testar meu procedimento, mas o resultado é insatisfatório, conforme ilustrado pelas duas figuras abaixo (a primeira é para as observações simuladas e a segunda figura mostra o histograma extraído do CDF estimado) :

figura 1 Figura 2

Existe alguém que sabe onde está o problema? Agradeço antecipadamente.

emberbillow
fonte
A amostragem de transformação inversa depende do uso do CDF inverso . pt.wikipedia.org/wiki/Inverse_transform_sampling
Sycorax diz Reinstate Monica
11
Seu estimador de densidade do kernel produz uma distribuição que é uma mistura local da distribuição do kernel; portanto, tudo o que você precisa para extrair um valor da estimativa de densidade do kernel é (1) extrair um valor da densidade do kernel e (2) selecionar independentemente um dos os dados apontam aleatoriamente e adicionam seu valor ao resultado de (1). Tentar inverter o KDE diretamente será muito menos eficiente.
whuber
@ Sycorax Mas eu realmente sigo o procedimento de amostragem por transformação inversa, conforme descrito no Wiki. Por favor, veja o código: p = rand; [~, idx] = classificação (abs (cdf (:, 2) - p)); rndval (i, 1) = cdf (idx (1), 1);
precisa saber é o seguinte
@whuber Não sei se meu entendimento de sua ideia está correto ou não. Por favor, ajude a verificar: primeiro redefina um valor das observações; e então desenhe um valor do kernel, digamos distribuição normal padrão; finalmente, adicione-os juntos?
emberbillow

Respostas:

12

Um estimador de densidade do kernel (KDE) produz uma distribuição que é uma mistura local da distribuição do kernel, portanto, para extrair um valor da estimativa de densidade do kernel, tudo o que você precisa fazer é: (1) extrair um valor da densidade do kernel e, em seguida, (2) selecione independentemente um dos pontos de dados aleatoriamente e adicione seu valor ao resultado de (1).

Aqui está o resultado desse procedimento aplicado a um conjunto de dados como o da pergunta.

Figura

O histograma à esquerda mostra a amostra. Para referência, a curva preta representa a densidade a partir da qual a amostra foi retirada. A curva vermelha representa o KDE da amostra (usando uma largura de banda estreita). (Não é um problema, ou mesmo inesperado, que os picos vermelhos sejam mais curtos que os pretos: o KDE espalha as coisas para que os picos fiquem mais baixos para compensar.)

O histograma à direita mostra uma amostra (do mesmo tamanho) do KDE. As curvas em preto e vermelho são as mesmas de antes.

Evidentemente, o procedimento usado para amostrar a partir da densidade funciona. Também é extremamente rápido: a Rimplementação abaixo gera milhões de valores por segundo a partir de qualquer KDE. Eu comentei bastante para ajudar na portabilidade para Python ou outras linguagens. O próprio algoritmo de amostragem é implementado na função rdenscom as linhas

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

rkerneldesenha namostras de iid da função kernel enquanto sampledesenha namostras com substituição dos dados x. O operador "+" adiciona as duas matrizes de amostras componente por componente.


KFKx=(x1 1,x2,,xn)

Fx^;K(x)=1 1nEu=1 1nFK(x-xEu).

XxEu1 1/nEuYX+YxX

FX+Y(x)=Pr(X+Yx)=Eu=1 1nPr(X+YxX=xEu)Pr(X=xEu)=Eu=1 1nPr(xEu+Yx)1 1n=1 1nEu=1 1nPr(Yx-xEu)=1 1nEu=1 1nFK(x-xEu)=Fx^;K(x),

como reivindicado.


#
# 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
fonte
Olá @whuber, quero citar essa ideia no meu trabalho. Você tem alguns artigos que foram publicados para isso? Obrigado.
emberbillow
2

Você colhe amostras do CDF primeiro, invertendo-o. O CDF inverso é chamado de função quantil; é um mapeamento de [0,1] para o domínio do VD. Em seguida, você amostra RVs uniformes aleatórios como percentis e os passa para a função quantil para obter uma amostra aleatória dessa distribuição.

AdamO
fonte
2
Esta é a maneira mais difícil: veja meu comentário à pergunta.
whuber
2
@whuber bom ponto. Sem estar muito envolvido nos aspectos programáticos, eu estava assumindo que deveríamos estar trabalhando com um CDF nesse caso. Sem dúvida, os internos para tal função dar uma densidade de kernel alisado e , em seguida, integrá-lo a obter uma CDF. Nesse ponto, provavelmente é melhor e mais rápido usar a amostragem por transformação inversa. No entanto, sua sugestão de usar apenas a densidade e a amostra diretamente da mistura é melhor.
AdamO 4/0118
@AdamO Obrigado pela sua resposta. Mas meu código realmente segue a mesma idéia que você disse aqui. Não sei por que os padrões tri-modais não podem ser reproduzidos.
precisa saber é o seguinte
@AdamO Aqui, se a palavra "internos" no seu comentário deve ser "intervalos"? Obrigado.
emberbillow
Ember, "internos" faz todo sentido para mim. Essa função precisa integrar a densidade da mistura e construir um inverso: esse é um processo confuso e numericamente complicado, como sugere o AdamO, e assim seria enterrado dentro da função - seus "internos".
whuber
1

Aqui, também quero postar o código do Matlab seguindo a idéia descrita pelo whuber, para ajudar aqueles que estão mais familiarizados com o Matlab do que com o 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;

O seguinte é o resultado: resultados

Diga-me se alguém encontrar algum problema com o meu entendimento e o código. Obrigado.

emberbillow
fonte
11
Além disso, descobri que meu código na pergunta está correto. A observação de que o padrão não pode ser reproduzido se deve em grande parte à escolha da largura de banda.
emberbillow
0

Sem olhar muito de perto para sua implementação, não obtenho completamente seu procedimento de indexação no ICDF. Eu acho que você tira do CDF, não é inverso. Aqui está a minha implementação:

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
fonte
2
Se você tem um cdf F, é verdade que F (X) é uniforme. Então você obtém o X pegando o cdf inverso de um número aleatório de uma distribuição uniforme. O problema que penso é como determinar o inverso quando você está produzindo uma densidade de kernel.
Michael R. Chernick
Obrigado pela sua resposta. Não colhi amostras diretamente do CDF. O código mostra que eu realmente fiz a mesma coisa que a amostragem por transformação inversa. p = margem; % esta linha obtém um número aleatório uniforme como probabilidade cumulativa. [~, idx] = classificação (abs (cdf (:, 2) - p)); rndval (i, 1) = CDF (IDX (1), 1);% estas duas linhas são para determinar o quantil correspondente à probabilidade cumulativa
emberbillow