'Returns the kernel density estimate for point x based on x-values
xs'

def density(self, xs, x):

    xs = np.asarray(xs)
    n = len(xs)
    if (self.weights is not None):
        (xs, weights) = self.in_domain(xs, self.weights, x)
    else:
        xs = self.in_domain(xs, xs, x)[0]
    xs = np.asarray(xs)
    if (xs.ndim == 1):
        xs = xs[:, None]
    if (len(xs) > 0):
        h = self.h
        if (self.weights is not None):
            w = ((1 / h) * np.sum((self(((xs - x) / h)).T * weights), axis=1))
        else:
            w = ((1.0 / (h * n)) * np.sum(self(((xs - x) / h)), axis=0))
        return w
    else:
        return np.nan
