10
$\begingroup$

I have some observations, and I want to mimick sampling based on these observations. Here I consider a non-parametric model, specifically, I use kernel smoothing to estimate a CDF from the limited observations.Then I draw values at random from the obtained CDF.The following is my code,(the idea is to get randomly a cumulative probability using uniform distribution, and take the inverse of the CDF with respect to the probability value)

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)

As shown in the code, I used a synthetic example to test my procedure, but the result is unsatisfactory, as illustrated by the two figures below (the first is for the simulated observations, and the second figure shows the histogram drawn from estimated CDF):

Figure 1 Figure 2

Is there anyone who knows where the problem is? Thank you in advance.

$\endgroup$
  • $\begingroup$ Inverse transform sampling hinges on using the inverse CDF. en.wikipedia.org/wiki/Inverse_transform_sampling $\endgroup$ – Sycorax Jan 4 '18 at 16:03
  • 1
    $\begingroup$ Your kernel density estimator produces a distribution that is a location mixture of the kernel distribution, so all you need to draw a value from the kernel density estimate is (1) draw a value from the kernel density and then (2) independently select one of the data points at random and add its value to the result of (1). Attempting to invert the KDE directly will be much less efficient. $\endgroup$ – whuber Jan 4 '18 at 17:18
  • $\begingroup$ @Sycorax But I indeed follow the inverse transform sampling procedure as described in Wiki. Please see the code: p = rand; [~, idx] = sort(abs(cdf(:, 2) - p)); rndval(i, 1) = cdf(idx(1), 1); $\endgroup$ – emberbillow Jan 5 '18 at 0:27
  • $\begingroup$ @whuber I am not sure whether my understanding of your idea is correct or not. Please help check: first resample a value from the observations; and then draw a value from kernel, say standard normal distribution; finally, add them together? $\endgroup$ – emberbillow Jan 5 '18 at 2:06
12
$\begingroup$

A kernel density estimator (KDE) produces a distribution that is a location mixture of the kernel distribution, so to draw a value from the kernel density estimate all you need do is (1) draw a value from the kernel density and then (2) independently select one of the data points at random and add its value to the result of (1).

Here is the result of this procedure applied to a dataset like the one in the question.

Figure

The histogram at the left depicts the sample. For reference, the black curve plots the density from which the sample was drawn. The red curve plots the KDE of the sample (using a narrow bandwidth). (It is not a problem, or even unexpected, that the red peaks are shorter than the black peaks: the KDE spreads things out, so peaks will get lower to compensate.)

The histogram at the right depicts a sample (of the same size) from the KDE. The black and red curves are the same as before.

Evidently, the procedure used to sample from the density works. It's also extremely fast: the R implementation below generates millions of values per second from any KDE. I have commented it heavily to assist in the porting to Python or other languages. The sampling algorithm itself is implemented in the function rdens with the lines

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

rkernel draws n iid samples from the kernel function while sample draws n samples with replacement from the data x. The "+" operator adds the two arrays of samples component by component.


For those who would like a formal demonstration of correctness, I offer it here. Let $K$ represent the kernel distribution with CDF $F_K$ and let the data be $\mathbf{x}=(x_1, x_2, \ldots, x_n)$. By the definition of a kernel estimate, the CDF of the KDE is

$$F_{\mathbf{\hat{x}};\, K}(x) = \frac{1}{n}\sum_{i=1}^n F_K(x-x_i).$$

The preceding recipe says to draw $X$ from the empirical distribution of the data (that is, it attains the value $x_i$ with probability $1/n$ for each $i$), independently draw a random variable $Y$ from the kernel distribution, and sum them. I owe you a proof that the distribution function of $X+Y$ is that of the KDE. Let's start with the definition and see where it leads. Let $x$ be any real number. Conditioning on $X$ gives

$$\eqalign{ F_{X+Y}(x) &= \Pr(X+Y \le x) \\ &= \sum_{i=1}^n \Pr(X+Y \le x \mid X=x_i) \Pr(X=x_i) \\ &= \sum_{i=1}^n \Pr(x_i + Y \le x) \frac{1}{n} \\ &= \frac{1}{n}\sum_{i=1}^n \Pr(Y \le x-x_i) \\ &= \frac{1}{n}\sum_{i=1}^n F_K(x-x_i) \\ &= F_{\mathbf{\hat{x}};\, K}(x), }$$

as claimed.


#
# 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))
$\endgroup$
  • $\begingroup$ Hi @whuber, I want to cite this idea in my paper. Do you have some papers that have been published for this? Thank you. $\endgroup$ – emberbillow Mar 3 '18 at 0:40
2
$\begingroup$

You sample from the CDF first by inverting it. The inverse CDF is called the quantile function; it is a mapping from [0,1] to the domain of the RV. You then sample random uniform RVs as percentiles and pass them to the quantile function to obtain a random sample from that distribution.

$\endgroup$
  • 2
    $\begingroup$ This is the hard way: see my comment to the question. $\endgroup$ – whuber Jan 4 '18 at 17:19
  • 2
    $\begingroup$ @whuber good point. Without being too enmeshed in the programmatic aspects, I was assuming we must be working with a CDF in this instance. No doubt the internals to such a function take a kernel smoothed density and then integrate it to obtain a CDF. At that point it is probably better and faster to use inverse transform sampling. However, your suggestion to just use the density and sample straight from the mixture is better. $\endgroup$ – AdamO Jan 4 '18 at 17:36
  • $\begingroup$ @AdamO Thank you for your answer. But my code indeed follows the same idea as you said here. I don't know why the tri-modal patterns can not be reproduced. $\endgroup$ – emberbillow Jan 5 '18 at 0:39
  • $\begingroup$ @AdamO Here whether the word "internals" in your comment should be "intervals"? Thank you. $\endgroup$ – emberbillow Jan 5 '18 at 2:19
  • $\begingroup$ Ember, "internals" makes perfect sense to me. Such a function has to integrate the mixture density and construct an inverse: that's a messy, numerically complicated process as AdamO hints, and so would be buried within the function--its "internals." $\endgroup$ – whuber Jan 5 '18 at 14:32
1
$\begingroup$

Here, I also want to post the Matlab code following the idea described by whuber, to help those who are more familiar with Matlab than 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;

The following is the result: results

Please tell me if anyone find any problem with my understanding and the code. Thank you.

$\endgroup$
  • 1
    $\begingroup$ Addtionally, I found that my code in the question is right. The observation that the pattern can not be reproduced is largely due to the choice of bandwidth. $\endgroup$ – emberbillow Jan 6 '18 at 3:34
0
$\begingroup$

Without looking too close into your implementation, I do not fully get your indexing procedure to draw from the ICDF. I think you draw from the CDF, not it's inverse. Here is my implementation:

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);
$\endgroup$
  • 2
    $\begingroup$ If you have a cdf F it is ttrue that F(X) is uniform. So you do get X by taking the inverse cdf of a random number from a uniform distribution. The problem I think is how to determine the inverse when you are producing a kernel density. $\endgroup$ – Michael Chernick Jan 4 '18 at 15:46
  • $\begingroup$ Thank you for your answer. I did not sample directly from the CDF. The code shows that I indeed did the same thing as inverse transform sampling. p = rand; % this line gets a uniform random number as the cumulative probability. [~, idx] = sort(abs(cdf(:, 2) - p)); rndval(i, 1) = cdf(idx(1), 1);% these two lines are to determine the quantile corresponding to the cumulative probability $\endgroup$ – emberbillow Jan 5 '18 at 3:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.