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:
@lorentzenchr Sorry, how is having a worse transformation, see the above plots, for sample sizes >> 1e4 a feature?
The current choice is a worse estimator for a performance gain , shouldn't the default not be best estimator and opt-in for worse but faster estimator?

Labels: Enhancement