KernelDensity incorrect handling of bandwidth

### Describe the bug

I was using kernel density estimator
https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html
using 'silverman' or 'scott' as the bandwidth argument. Then I found that the bandwidth automatically adjusted by the algorithm is independent of the actual scale of the dataset. In fact, I was shocked to find that the calculation of a bandwidth in https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/neighbors/_kde.py for 'silverman' and 'scott' does not check the scales of data at all. 

Suppose I fit the model `kde` to some 2D data `X` and get the bandwidth as `kde.bandwidth_`.
Next, I fit the model `kde` to the same 2D data `X` but with all elements multiplied by, say, 20 and get the bandwidth as `kde.bandwidth_`. 
I found that these two values of `kde.bandwidth_` are equal (it is calculated from the shape of `X`, see the source code). But obviously they should differ by a factor of 20 if the bandwidth is really computed in a truly adaptive manner.

For your reference, I want to mention that scipy's KDE https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html calculates the covariance of data to extract the scale of data. I think this is the right thing to do.

Note that if the bandwidth is incorrect, everything else is incorrect too, including probablities of samples, etc.

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.neighbors import KernelDensity
X = np.random.randn(1000, 2)
kde = KernelDensity(bandwidth='scott')
kde.fit(X)
print(kde.bandwidth_)

kde.fit(X * 20)
print(kde.bandwidth_)
```

### Expected Results

Different bandwidths for data sets with different scales.

### Actual Results

0.31622776601683794
0.31622776601683794

### Versions

```shell
1.2.1
```

COMMENT:
SciPy's factor and `KernelDensity.bandwidth_` are the same, but SciPy also takes into account the covariance when computing density. Here is a simpler snippet comparing just silverman between SciPy and scikit-learn:

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy import stats
from sklearn.neighbors import KernelDensity

N = 200
rng = np.random.default_rng(42)

X = np.concatenate(
    (rng.normal(0, 5, int(0.3 * N)), rng.normal(15, 5, int(0.7 * N)))
)
X_plot = np.linspace(-15, 30, 1000)

fig, ax = plt.subplots()
ax.hist(X, fc="black", alpha=0.2, label="input dist", bins=20, density=True)

kde_sklearn = KernelDensity(bandwidth="silverman").fit(X[:, None])
density_sklearn = np.exp(kde_sklearn.score_samples(X_plot[:, None]))

kde_scipy = stats.gaussian_kde(X, bw_method="silverman")
density_scipy = kde_scipy(X_plot)

for (name, density) in [("sklearn", density_sklearn), ("scipy", density_scipy)]:
    ax.plot(
        X_plot,
        density,
        lw=2,
        linestyle="-",
        label=f"{name} density",
    )

ax.legend(loc="upper left")
```

![Screenshot 2023-02-24 at 12 33 28 PM](https://user-images.githubusercontent.com/5402633/221248469-50ad6052-99d2-40f2-ac9f-141bcbdd0db0.jpg)

Concretely SciPy computes the factor and cholesky decomposition of the covariance [here](https://github.com/scipy/scipy/blob/96adfbae697c7972c6fbf986c3848f8d13fd228d/scipy/stats/_kde.py#L573) and use the decomposition to compute the density [here](https://github.com/scipy/scipy/blob/96adfbae697c7972c6fbf986c3848f8d13fd228d/scipy/stats/_kde.py#L268-L270).

I suspect transforming the data before passing it into `tree_` would be easier than adjusting `bandwidth`.

Labels: Bug, module:neighbors