3
$\begingroup$

Consider this samples set, feel free to look at it by using your own favourite tools (e.g. R, Python, whatever other stats tools).

The problem is that I lack experience in deciding whether any of these PDF estimations are reasonable, or whether there is even a possibly better estimation.

Thus the question in the title:

  • Q1: how to decide whether a KDE is good?

From this question raises another sub question:

  • Q2: is it a problem that my tails aren't touching 0 on the Y axis?

Using the data above, here are various KDEs that I made using 3 different kernels and bandwidths as follows: enter image description here enter image description here enter image description here enter image description here

And the code that generated those plots (replace PATH/TO/SAMPLES/FILE by the correct path to the samples set):

import numpy as np
from sklearn.neighbors.kde import KernelDensity
import matplotlib.pyplot as plt
import math

fle = 'PATH/TO/SAMPLES/FILE'
with open(fle, 'r') as f:
    X = [float(s) for s in f.read().splitlines()]
X.sort()

X_np = np.array(X)
X_np = X_np.reshape(-1,1)

# kernel density estimations
kernel='epanechnikov'
bw=0.001
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='gaussian'
bw=0.001
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='tophat'
bw=0.001
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

plt.show()


kernel='epanechnikov'
bw=0.01
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='gaussian'
bw=0.01
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='tophat'
bw=0.01
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

plt.show()


kernel='epanechnikov'
bw=0.05
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='gaussian'
bw=0.05
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='tophat'
bw=0.05
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

plt.show()


kernel='epanechnikov'
bw=0.09
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='gaussian'
bw=0.09
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

kernel='tophat'
bw=0.09
kde = KernelDensity(kernel=kernel, bandwidth=bw).fit(X_np)
plt.plot(X_np, np.exp(kde.score_samples(X_np)), label='%s, bw=%s' % (kernel, bw))
plt.legend(loc=0)

plt.show()
$\endgroup$
2
$\begingroup$

By eye, I'd say something between plots 2 and 3 would be good, but this is based on prior beliefs about what distributions tend to look like for data that I've encountered. For example, #1 looks far too spiky. But, it's certainly possible to draw data from a distribution that's actually like that, and the spiky estimate would then be the correct choice. I just believe a priori that it's probably not spiky. But in fact, among all possible distributions, the one with highest likelihood is the empirical distribution, which is a set of delta functions over the data points. Of course that doesn't seem reasonable, and the whole point of using a KDE is to impose something like a smoothness prior.

Instead of eyeballing it, it's often better to pick the kernel bandwidth using a principled approach. There are some simple formulas based on distributional assumptions. I prefer cross validation (this should be straightforward since you're using scikit-learn). Fit the KDE on the training set, evaluate the likelihood on the test set. Choose the bandwidth that maximizes the cross validated likelihood.

Regarding question 2, the tails will touch zero if you extend the range over which you evaluate the density.

$\endgroup$
  • $\begingroup$ While not very related (but I find it interesting) is your comment on the set of delta functions: do you mean the Dirac delta functions? If so, I assume the PDF that maximizes the likelihood would look like a comb (many tiny spikes) such that the length of each comb is $1/500$ (since I have $500$ data points) -- is this right? $\endgroup$ – caveman Jun 12 '16 at 10:23
  • 1
    $\begingroup$ Yes, weighted Dirac delta functions, each w/ integral 1/n (which incidentally is equivalent to the limit of kde w/ gaussian kernel as kernel width goes to 0). Just mentioned empirical distribution to say that training set likelihood isn't a good criterion for choosing kernel width $\endgroup$ – user20160 Jun 12 '16 at 15:47
  • 1
    $\begingroup$ Re:cross validation...Fit kde on the training set. Measure summed log probability of the test points using that kde. Sum over CV folds. Repeat for different kernel widths. Choose width that maximizes log probability of test set. $\endgroup$ – user20160 Jun 12 '16 at 16:06
  • 1
    $\begingroup$ Yes exactly right about the log. In your code, numpy.sum(kde.score_samples(xtest)) should give sum of log probabilities of test points (for test set vector xtest) $\endgroup$ – user20160 Jun 12 '16 at 18:34
  • 1
    $\begingroup$ Apologies, I was using incorrect terminology in the previous 2 comments. As you pointed out, they're densities and not probabilities bc we're talking about a continuous distribution. What we want to maximize is log likelihood of the kde given the test set, which is calculated as sum of log of kde evaluated at the test points. $\endgroup$ – user20160 Jun 13 '16 at 6:24

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.