QuantileTransformer's default subsampling introduces artefacts for unbounded distributions

### Describe the bug

The default behaviour of subsampling in the QuantileTransformer introduces artefacts when the input data originates from an unbounded distribution and the transformed dataset is (significantly) larger than the fitted one.

For any unbounded distribution the largest (smallest) percentiles will in expectation be underestimated (overestimated).
This effect is non-trivial for the selected subsample size of 10_000.

Here we draw a Multivariate Normal distribution of size 100_000
![mvn_input](https://github.com/scikit-learn/scikit-learn/assets/18510572/65676640-0674-46c8-a712-6e73e3a45b93)

Running `QuantileTransformer(output_distribution="normal").fit_transform(X)` i.e. with subsampling results in the below plot. Note the artefacts in the corners and the noisy marginal distributions.
![mvn_qt_defaults](https://github.com/scikit-learn/scikit-learn/assets/18510572/ddb42fb5-252f-4d62-80f4-2bbc2d12abf8)

Turning off the subsampling by specifying the subsampling size to be >= than the size of `X`. (this is undocumented btw)
`QuantileTransformer(output_distribution="normal", subsample=X.shape[0] + 1).fit_transform(X)`
![mvn_qt_no_sampling](https://github.com/scikit-learn/scikit-learn/assets/18510572/461601d5-9562-4230-a6f3-9a1081460cd2)

The effect of subsampling can be easily demonstrated by computing the difference over the subsample percentiles and the percentiles as you would estimate over the complete sample.
```Python
percentiles = np.arange(1, 100)
x = rng.normal(100_000)
xs = rng.choice(x, size=10_000, replace=False)
deltas = np.nanpercentile(x, percentiles) - np.nanpercentile(xs, percentiles)
```
The plot below shows the mean delta per percentile over 100_000 simulations of the above:
![standard_normal_subsample_delta](https://github.com/scikit-learn/scikit-learn/assets/18510572/c333c349-25f3-4386-8d2d-15149528d8a5)


### Steps/Code to Reproduce

```python
import numpy as np
import pandas as pd
import scipy.stats as sps
import seaborn as sns
from sklearn.preprocessing import QuantileTransformer

import matplotlib
import matplotlib.pyplot as plt

matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42
plt.rcParams['text.color'] = 'black'
plt.rcParams['figure.max_open_warning'] = 0
colors = [i['color'] for i in plt.rcParams['axes.prop_cycle']]
markers = ['o', 's', 'p', 'x', '^', '+', '*', '<', 'D', 'h', '>']


def plot_joint(X, ax=None):
    df = pd.DataFrame(X, columns=["x", "y"])
    axs = sns.JointGrid(data=df, x='x', y='y')
    axs.ax_joint.scatter(data=df, x='x', y='y', c=colors[0], marker='x', s=0.3)
    cnt, bins, _ = axs.ax_marg_x.hist(df['x'], bins=100, color=colors[0])
    cnt, bins, _ = axs.ax_marg_y.hist(df['y'], bins=100, color=colors[0], orientation='horizontal')
    axs.ax_joint.set_xlabel('')
    axs.ax_joint.set_ylabel('')
    return axs


def mvn_sample(rng, size=100_000, mu_x = 0., mu_y = 0., sigma_x = 1., sigma_y = 1., rho = 0.7):
    return rng.multivariate_normal(
        [mu_x, mu_y],
        [
            [sigma_x * sigma_x, rho * sigma_x * sigma_y],
            [rho * sigma_x * sigma_y, sigma_y * sigma_y]
        ],
        size=size
    )

rng = np.random.Generator(np.random.PCG64DXSM(753242))
X = mvn_sample(rng)

plot_joint(X).fig.show()
plot_joint(QuantileTransformer(output_distribution="normal").fit_transform(X)).fig.show()
plot_joint(QuantileTransformer(output_distribution="normal", subsample=X.shape[0] + 1).fit_transform(X)).fig.show()
```

### Expected Results

No artefacts introduced by subsampling, subsampling should not be on by default.

### Actual Results

See above.

### Versions

```shell
System:
    python: 3.11.3 (main, May  1 2023, 15:33:25) [Clang 14.0.3 (clang-1403.0.22.14.1)]
executable: ~/.pyenv/versions/3.11.3/envs/SKB/bin/python3.11
   machine: macOS-13.5.2-arm64-arm-64bit

Python dependencies:
      sklearn: 1.3.0
          pip: 22.3.1
   setuptools: 65.5.0
        numpy: 1.25.2
        scipy: 1.11.1
       Cython: None
       pandas: 2.0.3
   matplotlib: 3.7.2
       joblib: 1.3.2
threadpoolctl: 3.2.0

Built with OpenMP: True

threadpoolctl info:
       user_api: openmp
   internal_api: openmp
    num_threads: 12
         prefix: libomp
       filepath: ~/.pyenv/versions/3.11.3/envs/SKB/lib/python3.11/site-packages/sklearn/.dylibs/libomp.dylib
        version: None

       user_api: blas
   internal_api: openblas
    num_threads: 12
         prefix: libopenblas
       filepath: ~/.pyenv/versions/3.11.3/envs/SKB/lib/python3.11/site-packages/numpy/.dylibs/libopenblas64_.0.dylib
        version: 0.3.23.dev
threading_layer: pthreads
   architecture: armv8

       user_api: blas
   internal_api: openblas
    num_threads: 12
         prefix: libopenblas
       filepath: ~/.pyenv/versions/3.11.3/envs/SKB/lib/python3.11/site-packages/scipy/.dylibs/libopenblas.0.dylib
        version: 0.3.21.dev
threading_layer: pthreads
   architecture: armv8
```

COMMENT:
In my humble opinion, it's most certainly an unwanted behavior. We can not guarantee that sub-sampling does not change quantile estimates for sure (honestly I don't think that there is a mathematical sense in calculating percentiles on sampled data, they are not percentiles anymore). Artifacts presented in OP are just the part of the whole problem.

Simplest thing would be to compare whole transformed data and sub-sampled one.

```
import numpy as np
from scipy.stats import ks_2samp
from sklearn.preprocessing import QuantileTransformer

rng = np.random.default_rng(1)

LEN = 100_000
X = np.zeros(shape=(LEN, 3))
X[:,0] = rng.normal(0, 1, size=LEN)
X[:,1] = rng.normal(10, 0.1, size=LEN)
X[:,2] = rng.normal(-0.1, 10, size=LEN)

qt_subsample = QuantileTransformer(n_quantiles=1000, output_distribution='normal', subsample=10_000, random_state=1)
out_subsample = qt_subsample.fit_transform(X)

qt_full = QuantileTransformer(n_quantiles=1000, output_distribution='normal', subsample=100_000, random_state=1)
out_full = qt_full.fit_transform(X)

print('Mean abs percentage deviations')
print(np.abs((out_full - out_subsample) / out_full).mean(0))

print('For N(0, 1)')
res0 = ks_2samp(out_subsample[:,0], out_full[:,0])
print(f'KS test: {res0}')

print('For N(10, 0.1)')
res1 = ks_2samp(out_subsample[:,1], out_full[:,1])
print(f'KS test: {res1}')

print('For N(-0.1, 10)')
res2 = ks_2samp(out_subsample[:,2], out_full[:,2])
print(f'KS test: {res2}')
```

I think the least thing we can do is to change default to `None` and describe a trade-off between performance and estimation of real percentiles of the data.

Labels: Enhancement