How can I generate a CDF using Kernel Density Estimation in Python?

Question

There are a few methods I have come across that can do kernel density estimation which will provide a PDF for a sample of data:

KDEpy
sklearn.neighbors.KernelDensity
scipy.stats.gaussian_kde

Using any of the above I can generate a PDF however I want to know how I can get the CDF for the PDF I am generating. In math I know you can integrate on the PDF to get the CDF, however the issue is that these methods are only supplying x and y points and not a function to integrate on.

I'm wondering how I could transform the data being given into a CDF plot or alternatively find the PDF function for the data to then integrate on to get the CDF. Or use an alternative method where the output is a CDF instead of PDF.

user_3067485 · Answer

MCVE

Let's create some dummy data to shoulder the discussion:

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

np.random.seed(123)
data = stats.norm(loc=0, scale=1).rvs(10**4)

Here is the baseline idea with the scipy.stats package.

Gaussian KDE

We can estimate KDE using dedicated tools such as gaussian_kde:

kde = stats.gaussian_kde(data)

Which exposes a PDF function to evaluate at every x but is missing the CDF.

Checking samples with Kolmogorov-Smirnov Test we cannot reject the null hypothesis (two distributions are identical) with the threshold of 10%:

stats.ks_2samp(data, kde.resample(100).squeeze())
# KstestResult(statistic=0.0969, pvalue=0.29163373800871994)

Continuous Variable

The scipy.stats package also exposes a generic class rv_continous to inherit from. As stated in documentation:

New random variables can be defined by subclassing the rv_continuous class and re-defining at least the _pdf or the _cdf method (normalized to location 0 and scale 1).

So we can use this on purpose logic to fill the gap. Without any performance consideration it boils down to:

class KDEDist(stats.rv_continuous):
    
    def __init__(self, kde, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self._kde = kde
    
    def _pdf(self, x):
        return self._kde.pdf(x)

Then we create the underlying object with our experimental KDE.

X = KDEDist(kde)

stats.ks_2samp(data, X.rvs(size=100))  # This call is kind of intensive
# KstestResult(statistic=0.0625, pvalue=0.8113077271721811)

Now we can naturally - at least in term of API call - evaluate the PDF and CDF as well:

fig, axe = plt.subplots()
axe.hist(data, density=1)
axe.plot(x, X.pdf(x))
axe.plot(x, X.cdf(x))

It returns:

Performance considerations

Notice than this methodology answers your question but is not performant. KDE computation are expensive mainly because the kernel spans the whole data space (Gaussian reaches zero at infinity). Therefore, without cut-off feature computations are based on all observations of the dataset at each evaluation.

Changing the window function can drastically improve the performance. Eg.: triangular window will have fixed span over the whole dataset and reduce computation w.r.t. dataset extent and size.

Implementation considerations

Reading the doc, it seems rv_continuous is initially designed to implement new continuous variable with analytical definition.

Anyway, the class provides automatic resolution/integration for other statistics if underlying methods are not implemented (overridden).

When choosing this methodology, it is up to you to implement missing logic if you wish to make it more performant and robust (numerical stability).

Histogram instead of KDE

If you can relax the KDE needs and is satisfied by an histogram distribution, then you can rely on rv_histogram which essentially does the same based on the binned distribution:

hist = np.histogram(data, bins=100)
hist_dist = stats.rv_histogram(hist)

stats.ks_2samp(data, hist_dist.rvs(size=100))
# KstestResult(statistic=0.0577, pvalue=0.8778871545532821)

Which returns:

KDE Histogram

Provided it is acceptable theoretically, we can mix both strategy by creating the expected histogram from the KDE:

hist = np.histogram(data, bins=1000)
hist_kde = kde.pdf(hist[1][:-1] + np.diff(hist[1]))
hist_dist_kde = stats.rv_histogram([hist_kde, hist[1]])

stats.ks_2samp(data, hist_dist_kde.rvs(size=100))
# KstestResult(statistic=0.1067, pvalue=0.19541766226890545)

Then the CDF has a relative smoothness w.r.t. the KDE (it is still an histogram) and the Continuous Variable object is as performant as rv_histogram can be.