Automatic bandwidth calculation valid only for normalized data #26658
Comments
|
For the record, here is the generated plot: 
|
|
Related issue #25623 |
|
I agree it looks like a bug, please feel free to open a PR and rework your reproducer as a non-regression test to check that scikit-learn's KDE approximately match scipy's implementation on this toy data. Please check that your solution is a fix for #25623. |
|
Hi @ogrisel, not sure who exactly to pin but after some investigation I want to compare two approaches here:
Correctness
Multi-dimensional data
Effects on other public interfacesOne public interface that I notice is the
Other issuesFor both approaches, especially Approach 2, I think we should change the default to SummaryPersonally I would prefer Approach 2, and I tried a "hack" implementation locally and seems that it works well and almost completely agrees with the scipy result when using Gaussian kernel. However I cannot guarantee that my analysis is completely correct. What do maintainers think about this? P.S. If any further information is needed, please let me know. AttachmentsTesting code (1D)import os
from itertools import product
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde, norm
from sklearn.neighbors import KernelDensity
if __name__ == "__main__":
savepath = os.path.join(os.path.dirname(__file__), "../figures")
if not os.path.exists(savepath):
os.mkdir(savepath)
for scale, bw in product([0.01, 1, 100], ["scott", 1.0]):
data1 = np.random.normal(5 * scale, scale, 10000)
data2 = np.concatenate(
[np.random.normal(0, scale, 3000), np.random.normal(5 * scale, scale, 7000)]
)
xs1 = np.linspace(start=0, stop=10 * scale, num=1000)
xs2 = np.linspace(start=-5 * scale, stop=10 * scale, num=1000)
ref1 = norm(5 * scale, scale).pdf(xs1)
ref2 = 0.3 * norm(0, scale).pdf(xs2) + 0.7 * norm(5 * scale, scale).pdf(xs2)
for i, (data, xs, ref) in enumerate(
zip([data1, data2, data3], [xs1, xs2, xs3], [ref1, ref2, ref3])
):
kd_sklearn = KernelDensity(kernel="gaussian", bandwidth=bw).fit(data.reshape(-1, 1))
sklearn_result = np.exp(kd_sklearn.score_samples(xs.reshape(-1, 1)))
sklearn_bw = kd_sklearn.bandwidth_
kd_scipy = gaussian_kde(data, bw_method=bw)
scipy_result = kd_scipy.pdf(xs)
scipy_bw = kd_scipy.factor
plt.plot(xs, scipy_result, label=f"scipy ({scipy_bw:.3f})", lw=5, alpha=0.3)
plt.plot(xs, sklearn_result, label=f"sklearn ({sklearn_bw:.3f})")
plt.fill(xs, ref, alpha=0.1, color="black")
plt.legend()
plt.title(f"scale={scale}, bw={bw}")
plt.savefig(os.path.join(savepath, f"{scale}-{bw}-data{i + 1}.png"))
# plt.show()
print("\033[92m.\033[0m", end="", flush=True)
plt.clf()Figures produced by my hack implementation (using the testing code above, 1D)
Figures produced by the current main branch (using the testing code above, 1D, zipped)Download: main-1d-test.zip 2D exampleThe data for this example is simply stacking the two univariate data, the first with
To further validate the result, scipy has the
|
|
@Charlie-XIAO thanks for the summary and experiments, and sorry for the late feedback. I agree with your conclusions. Would you mind opening a PR with approach 2? |
|
Note that the 2 approaches are not necessarily mutually exclusive. In the future we could decide to implement the variance based approach as well if someone can highlight a real-life case where approach 1 would perform significantly better than approach two. |
|
@ogrisel Thanks for your reply! Unfortunately I’m kinda busy recently, but will revisit and make a PR in a month or so. |
Describe the bug
sklearn.neighbors.KernelDensitysupports automatic (optimal) bandwidth calculation viabandwidth = 'silverman'andbandwidth = 'scott'. The algorithm computes the appropriate observation-weighted bandwidth factors (proportional to nobs^0.2) but does not adjust for the standard deviation or interquartile range of the dataset. Roughly, the algorithm should scale the dataset's standard error by the algorithmic bandwidth factors.See, e.g., Wikipedia. The implementation in
scipy.stats._kdeis correct.Steps/Code to Reproduce
Expected Results
Automatic SKLearn bandwidth curve should approximately match SciPy bandwidth curve, roughly the shape of the underlying data histogram.
Actual Results
Automatic SKLearn bandwidth curve generates a flat PDF.
Versions
System: python: 3.9.5 (default, Nov 23 2021, 15:27:38) [GCC 9.3.0] executable: /local_disk0/.ephemeral_nfs/envs/pythonEnv-67c47d19-3f15-49a2-ab8f-bcf25a2bc29f/bin/python machine: Linux-5.15.0-1038-azure-x86_64-with-glibc2.31 Python dependencies: sklearn: 1.2.2 pip: 21.2.4 setuptools: 58.0.4 numpy: 1.20.3 scipy: 1.7.1 Cython: 0.29.24 pandas: 1.3.4 matplotlib: 3.4.3 joblib: 1.2.0 threadpoolctl: 2.2.0 Built with OpenMP: True threadpoolctl info: filepath: /databricks/python3/lib/python3.9/site-packages/numpy.libs/libopenblasp-r0-5bebc122.3.13.dev.so prefix: libopenblas user_api: blas internal_api: openblas version: 0.3.13.dev num_threads: 6 threading_layer: pthreads architecture: Haswell filepath: /local_disk0/.ephemeral_nfs/envs/pythonEnv-67c47d19-3f15-49a2-ab8f-bcf25a2bc29f/lib/python3.9/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0 prefix: libgomp user_api: openmp internal_api: openmp version: None num_threads: 6 filepath: /databricks/python3/lib/python3.9/site-packages/scipy.libs/libopenblasp-r0-085ca80a.3.9.so prefix: libopenblas user_api: blas internal_api: openblas version: 0.3.9 num_threads: 6 threading_layer: pthreads architecture: HaswellThe text was updated successfully, but these errors were encountered: