Incorrect results in kernel density estimation (```KernelDensity```) when the kernel width is small

### Describe the bug

I have met an issue when using the class ```sklearn.neighbors.KernelDensity```. The **kernel density estimation (KDE)** works unusually when the kernel width is very small (less than 0.01),e.g., 0.00001. The output kde scores from the function ```score_samples``` are all big constants and almost the same, which is obviously very strange, not likely to be the correct densities. 
For example in the sample code, I generated 1000 Gaussian distributed data, used it to fit the KDE model with a small kernel width (0.00001), and again evaluated the model on the same data. Then, the outputs from ```score_samples``` are almost 3.6862 and the prob are 39.8942 after applying a ```np.exp()``` on top of it. The smaller the kernel width is, the worse KDE scores would be.  
So is the KDE itself still correct when the kernel width is very small, or perhaps is there any issue with the normalization in the implementation? Thank you very much for checking this issue and the help!

### Steps/Code to Reproduce

```python
from sklearn.neighbors import KernelDensity
import numpy as np

# generate training data
rand = np.random.RandomState(seed=100)
data = rand.normal(0, 1, 1000).reshape(-1, 1)

# fit the training data and evaluate on it
# kde scores tend to be probmatic when bandwidth is very small, e.g., 0.00001
model = KernelDensity(bandwidth=0.00001, kernel='gaussian')
model.fit(data)
log_dens = model.score_samples(data)
print(log_dens)
print(np.exp(log_dens))

```

### Expected Results

The results shouldn't be like such big constants and almost the same.

### Actual Results

When the kernel width is 0.00001, the output from ```score_samples``` for the first 20 data:
 ``` 
[3.68623165 3.68623165 3.68623165 3.68623165 3.68623165 3.68623165
 3.68623165 3.68623165 3.68623165 3.68623165 3.68623165 3.68623165
 3.68623165 3.68623165 3.68623165 3.68623165 3.68623165 3.68623165
 3.97824858 3.68623165]
```
After applying ```np.exp()``` on top of the above results:
```
[39.89422804 39.89422804 39.89422804 39.89422804 39.89422804 39.89422804
 39.89422804 39.89422804 39.89422804 39.89422804 39.89422804 39.89422804
 39.89422804 39.89422804 39.89422804 39.89422804 39.89422804 39.89422804
 53.42338568 39.89422804]
```
The results is even worse when the kernel width is smaller, e.g., when the kernel width is 0.000001, the results would be:
Output from ```score_samples```:
```
[5.98881675 5.98881675 5.98881675 5.98881675 5.98881675 5.98881675
 5.98881675 5.98881675 5.98881675 5.98881675 5.98881675 5.98881675
 5.98881675 5.98881675 5.98881675 5.98881675 5.98881675 5.98881675
 5.98881675 5.98881675]
```
After ```np.exp()```:
```
[398.9422804 398.9422804 398.9422804 398.9422804 398.9422804 398.9422804
 398.9422804 398.9422804 398.9422804 398.9422804 398.9422804 398.9422804
 398.9422804 398.9422804 398.9422804 398.9422804 398.9422804 398.9422804
 398.9422804 398.9422804]
```


### Versions

```shell
System:
    python: 3.9.16 (main, Jan 11 2023, 10:02:19)  [Clang 14.0.6 ]
executable: /Users/tongtong/opt/anaconda3/envs/torch/bin/python
   machine: macOS-10.16-x86_64-i386-64bit

Python dependencies:
      sklearn: 1.2.1
          pip: 22.3.1
   setuptools: 65.6.3
        numpy: 1.23.5
        scipy: 1.10.1
       Cython: None
       pandas: 1.5.2
   matplotlib: 3.6.2
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: openmp
   internal_api: openmp
         prefix: libomp
       filepath: /Users/tongtong/opt/anaconda3/envs/torch/lib/python3.9/site-packages/sklearn/.dylibs/libomp.dylib
        version: None
    num_threads: 8

       user_api: blas
   internal_api: mkl
         prefix: libmkl_rt
       filepath: /Users/tongtong/opt/anaconda3/envs/torch/lib/libmkl_rt.1.dylib
        version: 2021.4-Product
threading_layer: intel
    num_threads: 8

       user_api: openmp
   internal_api: openmp
         prefix: libomp
       filepath: /Users/tongtong/opt/anaconda3/envs/torch/lib/libomp.dylib
        version: None
    num_threads: 8

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /Users/tongtong/opt/anaconda3/envs/torch/lib/python3.9/site-packages/scipy/.dylibs/libopenblas.0.dylib
        version: 0.3.18
threading_layer: pthreads
   architecture: Nehalem
    num_threads: 8
```

COMMENT:
> What makes you think this is a bug? If you fit a KDE model with near-zero bandwidth you fit a near mixture of Diracs model where the log-likelihood is not defined anymore (on the training data), right?

Thank you very much for the reply! So you are not really estimating the density (non-negative and integrates to 1) but the likelihood of the observed data as a function of parameter (the integral does not necessarily equal 1); if so, why not normalize the latter to form a probability density function? Sorry for any potential misunderstanding.

Labels: Needs Info