UX `CalibrationDisplay`'s naive use can lead to very confusing results

The naive use of `CalibrationDisplay` parameter silently leads to degenerate, noisy results when some bins have with a few data points.

For instance, look at the variability obtained by displaying for calibration curve of a fitted model evaluated on various resampling with 50 data points in total using the uniform strategy when using `n_bins=10` and the default `strategy="uniform"`:

![Image](https://github.com/user-attachments/assets/0a925e43-0466-47aa-b57b-519da3b61b5a)


<details>

```python
# %%
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.calibration import CalibrationDisplay
from sklearn.model_selection import train_test_split


X, y = make_classification(n_samples=10_000, n_features=200, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# %%
clf = LogisticRegression(C=10).fit(X_train, y_train)
# %%

fig, ax = plt.subplots()
for seed in range(5):
    # resample the test set
    rng = np.random.RandomState(seed)
    indices = rng.choice(X_test.shape[0], size=50, replace=True)
    X_test_sample = X_test[indices]
    y_test_sample = y_test[indices]
    CalibrationDisplay.from_estimator(clf, X_test_sample, y_test_sample, n_bins=10, ax=ax, label=None)
```

</details>


This problem can easily happen with the default `strategy="uniform"` if the test data is not large enough. I think this class should warn the user whenever it generates bins with lower than 10 data points per bin.

A typical user will only get one of the curves above and not suspect that it's just noise without manually plotting the others by random resampling. Note that I chose a minimal test set to make the problem catastrophic above, but it can happen with larger sample sizes, in particular with the uniform strategy, in particular on imbalanced datasets.

## Updated recommendations

EDIT: based on the discussion below, here are my recommendations to address this issue:

- display a confidence interval by default to reflect our uncertainty on the calibration curve estimate (using the Clopper-Pearson analytical CI from scipy looks both reliable and cheap to compute)
- we set `n_bins=15` by default but also introduce `min_samples_per_bins=50` to automatically merge the smallest pair of neighboring bins until that constraint it met.
- we would not warn, unless the resulting merging leads to less than 3 bins for instance,
- we could keep `strategy="uniform"` if we like, but I wouldn't mind switching to `"quantile"` either.

## Original description and suggestion

In such cases, this warning could recommend to:

- use `strategy="quantile"` (which is a bit less likely to generate near empty, unstable bins);
- decrease the number of bins;
- increase the number of data points in the test set.
  
Beyond warnings, we could also offer the user extra-options to assess the stability of the calibration curve.

There are two cases:

- for a fixed trained estimator, we could estimate the uncertainty in the estimation of its calibration curve induced by the finite size of the evaluation set (irrespective of whether the training procedure is stable or not) by resampling the test set (with replacement):

```python
fig, ax = plt.subplots()
for seed in range(10):
    rng = np.random.RandomState(seed)
    indices = rng.choice(X_test.shape[0], size=X_test.shape[0], replace=True)
    X_test_resampled = X_test[indices]
    y_test_resampled = y_test[indices]

    CalibrationDisplay.from_estimator(
        clf,
        X_test_resampled,
        y_test_resampled,
        strategy="quantile",
        n_bins=10,
        ax=ax,
        label=None,
    )
```
 We could optionally offer an option to aggregate those curves by displaying a shaded region derived from the vertical quantiles (e.g. 0.05-0.95) of the interpolated y-curve values at a fine grid of x values. See also the related issue on sampling uncertainties in performance curves: #25856.

- to assess the uncertainty of the calibration of a full fitting procedure (uncertainty induced by both the finite training set size, the training algorithm itself (impacted by the choice of hyperparams) and the finite evaluation set size), we could also offer to plot calibration curves from CV results (see #21211 to implement `CalibratedDisplayCV.from_cv_results`). 


On the opposite, if all the bins have more than 1000 data points and `n_bins <= 10`, we could suggest the user to increase n_bins to 10 or more to get a finer estimation of the calibration curve. Maybe this recommendation is better to be set in a docstring rather than a warning.


I think we could start simple by opening a first PR to issue warnings for the nearly empty bin problem (with an option to silence them).

COMMENT:
I found that [`scipy.stats.binomtest`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binomtest.html) implements various methods to estimate [CI intervals](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats._result_classes.BinomTestResult.proportion_ci.html#scipy.stats._result_classes.BinomTestResult.proportion_ci) for binomial trials:

```python
>>> for y_obs in [[0], [1], [0, 1], [0, 0, 1], [0, 1, 1]]:
...     print(f"{y_obs=}")
...     for method in ["exact", "wilson", "wilsoncc"]:
...         print(method.rjust(10), binomtest(sum(y_obs), len(y_obs)).proportion_ci(confidence_level=0.5, method=method))
...     print()
y_obs=[0]
     exact ConfidenceInterval(low=0.0, high=0.75)
    wilson ConfidenceInterval(low=0.0, high=np.float64(0.3126847440825833))
  wilsoncc ConfidenceInterval(low=0.0, high=np.float64(0.7795911050456467))

y_obs=[1]
     exact ConfidenceInterval(low=0.25, high=1.0)
    wilson ConfidenceInterval(low=np.float64(0.6873152559174167), high=1.0)
  wilsoncc ConfidenceInterval(low=np.float64(0.22040889495435334), high=1.0)

y_obs=[0, 1]
     exact ConfidenceInterval(low=0.13397459621556132, high=0.8660254037844387)
    wilson ConfidenceInterval(low=np.float64(0.2847588845819623), high=np.float64(0.7152411154180377))
  wilsoncc ConfidenceInterval(low=np.float64(0.10425347838136545), high=np.float64(0.8957465216186346))

y_obs=[0, 0, 1]
     exact ConfidenceInterval(low=0.0914397035839695, high=0.6736481776669304)
    wilson ConfidenceInterval(low=np.float64(0.18281725015861316), high=np.float64(0.5277418231406089))
  wilsoncc ConfidenceInterval(low=np.float64(0.06837944703876814), high=np.float64(0.681436779552594))

y_obs=[0, 1, 1]
     exact ConfidenceInterval(low=0.32635182233306964, high=0.9085602964160305)
    wilson ConfidenceInterval(low=np.float64(0.47225817685939114), high=np.float64(0.8171827498413868))
  wilsoncc ConfidenceInterval(low=np.float64(0.3185632204474061), high=np.float64(0.9316205529612319))
```

The `"wilsoncc"` method seems to be the most conservative for small sample sizes followed by the `"exact"` method.

Labels: Enhancement, module:inspection