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:
> @RUrlus Do you have a concrete example where the current default is hurting the performance of predictive models?

@glemaitre It's a bit of work but I can make an example. The reason behind this issue is that we were revamping an approach we used for [data synthesis](https://proceedings.mlr.press/v151/baak22a/baak22a.pdf) and noticed that performance had degraded significantly after updating dependencies. The cause turned out to be the change in the default number of samples to use for subsampling, which decreased from 1e5 to 1e4 about a year ago.

Labels: Enhancement