Sampling uncertainty on precision-recall and ROC curves

### Describe the workflow you want to enable

We would like to add the possibility to plot sampling uncertainty on precision-recall and ROC curves.


### Describe your proposed solution

We (@mbaak, @RUrlus, @ilanfri and I) published a paper in [AISTAT 2023](http://aistats.org/aistats2023/accepted.html) called [Pointwise sampling uncertainties on the Precision-Recall curve](https://github.com/RUrlus/ModelMetricUncertaintyResearch/tree/sklearn_pull_request/notebooks/pr_ellipse_validation/1411.pdf), where we compared multiple methods to compute and plot them.

We found out that a great way to compute them is to use profile likelihoods based on Wilks’ theorem.
It consists of the following steps:

1. Get the curve
2. Get the confusion matrix of each point of the curve
3. For each observed point of the curve, estimate a surrounding 6 (i.e. more than the desired number) sigmas uncertainty grid rectangle (based on first-order approximation of the covariance matrix, with the bivariate normal distribution assumption) 
4. For each of these hypothesis point in the grid, compute the test static with the observed point, called the profile log likelihood ratio (using the fact that the confusion matrix follows a multinomial distribution).
5. Plot the 3 sigmas contour (i.e. isoline) for the observed points (using Wilks’ theorem stating that the profile log likelihood ratio is described asymptotically by a chi2 distribution)


We have a minimal pure Python implementation:
https://github.com/RUrlus/ModelMetricUncertaintyResearch/blob/sklearn_pull_request/notebooks/pr_ellipse_validation/demo_ROC_PR_curves_sklearn_pull_request.ipynb 

And a C++ implementation: the paper is supported by our package [ModelMetricUncertainty](https://github.com/RUrlus/ModelMetricUncertainty) which has a C++ core with, optional, OpenMP support and Pybind11 bindings. Note that this package contains much more functionality than the above notebook. The core is binding agnostic allowing a switch to Cython if needed. Upside is that it is much faster (multiple orders) than the above Python implementation at the cost of complexity.

The pure Python implementation would look like this:
![image4](https://user-images.githubusercontent.com/1184396/225160909-f03bc592-ffb0-4b04-8f6a-eadd889c65a1.png)

I’m also suggesting other visual improvements:
1. Add x and y axis limit: [0, 1], in sklearn axes currently start at ~-0.1
2. Modify the plotting frame: either remove the top and right lines to see the curve better when values are close to 1, or plot the frame with a dotted line
3. Fix aspect ratio to squared, since the two axes are the same scale.

With those it can look like this:
![download-3](https://user-images.githubusercontent.com/1184396/225161339-eb794d06-073e-4c4b-932e-618227e4925a.png)
![download-2](https://user-images.githubusercontent.com/1184396/225161357-e56c8077-271b-4ae2-b92d-1770c17f538c.png)

Remark: I set the contour color to lightblue, let me know if it is fine.

We need to align on the API integration. I suggest adding some parameters in `PrecisionRecallDisplay` and in `RocCurveDisplay` called:

- `uncertainty=True` to enabel plot uncertainty band (or `plot_uncertainty_style=` ?)
- `uncertainty_n_std=3` to decide how +/- standard deviation the band should be
- `uncertainty_n_bins=100` to decide how fine-grained the band should be (see remark about running time)

### Describe alternatives you've considered, if relevant

Other ways to compute uncertainties are evaluated in our paper.

We have noticed that there is open pull request on related topic: https://github.com/scikit-learn/scikit-learn/pull/21211
That is great, however cross-validation covers different sources of uncertainties, and has some limitations (a bias is introduced by overlapping training folds, introducing a correlation in the trained models. In addition, this uncertainty depends on the size of a fold, and is likely larger than on the test set, see [ref](https://www.jmlr.org/papers/volume5/grandvalet04a/grandvalet04a.pdf).)


### Additional context

##### Running time discussion
Here is an analysis of the running time of this pure Python method:

<img width="676" alt="image5" src="https://user-images.githubusercontent.com/1184396/225161990-44fae526-e25d-4129-b03f-0ce24e2df530.png">


The execution time depends on the number of points (i.e. thresholds) plotte and on `uncertainty_n_bins`.
With a surrounding grid of `uncertainty_n_bins=100` per point it is fast enough and fine enough.
There is barely any noticeable visual difference between 50 and 100 (or more) points (at least in this example), [see curves](https://github.com/RUrlus/ModelMetricUncertaintyResearch/tree/sklearn_pull_request/notebooks/pr_ellipse_validation/curves).
For let’s say for a 100k set, it is too slow for ROC, because there is much more thresholds, but this is going to be fixed soon here https://github.com/scikit-learn/scikit-learn/pull/24668 . But anyway, in this case, the uncertainties are really small, so plotting them doesn’t really make.

COMMENT:
> I think this would be a useful tool to have. As I understand it this tool allows you to estimate how well you know the performance of your fitted estimator, given the test dataset. This means it should be easy to demonstrate in an example that the uncertainty reduces as you increase the size of the test dataset (keeping everything else fixed).
> 
> If the above understanding is correct then I'd vote for adding this.

@betatim you're correct. A quick example:
```Python
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

import mmu
from mmu.viz.utils import _set_plot_style

_ = _set_plot_style()

seeds = mmu.commons.utils.SeedGenerator(2343451)
X, y = make_classification(
    n_samples=2000, n_classes=2, random_state=seeds()
)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.5, random_state=seeds()
)
model = LogisticRegression(solver='lbfgs')
model.fit(X_train, y_train)
y_score = model.predict_proba(X_test)[:, 1]

err_500 = mmu.PRU.from_scores(y=y_test[:500], scores=y_score[:500])
err_1000 = mmu.PRU.from_scores(y=y_test[:1000], scores=y_score[:1000])
ax = err_500.plot(other=err_1000)
```

produces
![size_effect](https://user-images.githubusercontent.com/18510572/231157157-45e02f31-ce35-4a1a-8400-225d89b33125.png)

Labels: New Feature, Needs Decision - Include Feature