Q:

Sci-kit Learn KDE variance increases by one

I'm trying to fit a kernel density estimate for some normal data with sklearn's kde.  Here's an example:
import numpy as np
from sklearn.neighbors.kde import KernelDensity as kde

x = np.random.normal(scale = 2, size = [10000, 1])
np.var(x) # 4.0
test_kde = kde()
test_kde.fit(x)

np.var(test_kde.sample(10000)) # 5.0

The variance increases by one.  Am I doing something incredibly stupid here? 

A:

The problem is that you don't specify the correct bandwidth to scale the individual density functions which is why you are over-smoothing the estimated density function.  Since your example data follows a normal distribution, a bandwidth of
>>> h = ((4 * np.std(x)**5) / (3 * len(x)))**(1/5)
>>> h
0.33549590926904804

would be optimal.  An explanation can be found on Wikipedia. 
>>> test_kde = kde(bandwidth=h)
>>> test_kde.fit(x)
>>> samples = test_kde.sample(10000)
>>> np.var(samples)
4.068727474888099   # close enough to 4

But why do I need such scaling? 
Kernel Density Estimation works by using a kernel function (often the density function of the normal distribution) to estimate the density of your data distribution.  The general idea is that by summing many density functions parametrized by your sample will eventually, given enough samples, approximate the original density function:

We can visualize this for your data:
from matplotlib.colors import TABLEAU_COLORS

def gauss_kernel(x, m=0, s=1): 
    return (1/np.sqrt(2 * np.pi * s**2) * np.e**(-((x - m)**2 / (2*s**2))))

from matplotlib.colors import TABLEAU_COLORS

x_plot = np.linspace(-2, 2, 10)
h = 1

for xi, color in zip(x_plot, TABLEAU_COLORS.values()):
    plt.plot(xi, gauss_kernel(xi, m=0, s=2) * 0.001, 'x', color=color)
    plt.plot(x, 1 / (len(x) * h) * gauss_kernel((xi - x) / h), 'o', color=color)
    plt.plot(xi, (1 / (len(x) * h) * gauss_kernel((xi - x) / h)).sum() * 0.001, 'o', color=color)

This plot shows the estimated and true density at some points in [-2; 2] along with the kernel function for each point (same colored curves).  The estimated density is simply the sum of the corresponding kernel function. 
One can see that the farther right/left the individual kernel functions are, the lower their sum is (and consequently the density).  To explain this you have to remember that our original data points are centered around 0 since they are sampled from a normal distribution with mean 0 and variance 2.  Therefore, the farther from the center, the less data points > 0 are there.  Consequently, this means that a Gaussian kernel function that takes these points as input will eventually have all data points in one of its flat tail sections and weight them very close to zero, which is why the sum of this kernel function will be very small there.  One could also say that we are windowing our data points using a Gaussian density function. 
You can see the impact of the bandwidth parameter clearly by setting h=2:
h = 2

for xi, color in zip(x_plot, TABLEAU_COLORS.values()):
    plt.plot(xi, gauss_kernel(xi, m=0, s=2) * 0.001, 'x', color=color)
    plt.plot(x, 1 / (len(x) * h) * gauss_kernel((xi - x) / h), 'o', color=color)
    plt.plot(xi, (1 / (len(x) * h) * gauss_kernel((xi - x) / h)).sum() * 0.001, 'o', color=color)

The individual kernel functions are a lot smoother and, as a consequence, the estimated 
density is a lot smoother, too.  The reason for this is in the formulation of the 
smoothing operator.  The kernel is called as 
1/h K((x - xi)/h)

which in case of a Gaussian kernel means to compute density of a normal distribution with a mean of xi and a variance of h. Therefore: the higher h, the smoother each density estimate! 
In case of sklearn, the best bandwidth can be estimated by using, for example, grid search by measuring the quality of the density estimate.  This example shows you how.  If you have chosen a good bandwidth you can estimate the density function pretty well:

