API to predict multiple quantiles at once

Classifiers have a `predict_proba` method that makes it possible to quantify probabilistic ally the certainty in the predictions for a given input `X_i`.

Currently most regressors in scikit-learn only predict a conditional expectile E[Y|X], and some have a `return_std` option that makes it also possible to estimate sqrt(VAR[Y|X]), which can be used to quantify the certainty when assuming a Gaussian predictive distribution (typically for Gaussian processes which estimate a Gaussian predictive posterior distribution).

We do have pointwise quantile estimators (linear models, gradient boosting, hist gradient boosting) where the `predict` method returns a single point estimate for target quantile passed as an hyper-parameter instead of estimating an expectile.

Several people have expressed the need to have more generic API that can return an array of quantile estimates for a given input `X_i`.

The goal of this issue is to centralize the discussion of an API extension to be able to do this more uniformly in scikit-learn, either via a meta-estimator that wraps an array of point-wise quantile estimator to turn it into a quantile-array estimator or to directly have the base estimators able to do this directly (and sometimes more efficiently).

Some non-exhausitive list of related PRs and issues (feel free to add or suggest new ones):

- https://github.com/scikit-learn/scikit-learn/issues/19851
- https://github.com/scikit-learn/scikit-learn/pull/19754
- https://github.com/scikit-learn/scikit-learn/issues/20964

Also related:

- conformal predictions: https://github.com/scikit-learn-contrib/MAPIE
- XGBoostLSS (location, scale and shape) https://github.com/StatMixedML/XGBoostLSS

Furthermore, models like Poisson regression that make a specific assumption about the conditional Y|X distribution, it would be possible to estimates of the inverse-CDF values of the estimated Y|X for instance. Those could probably also benefit from an expanded API.

If we do this, then we have the side question of how to evaluate such multi-quantile models. We could probably extend the pinball_loss scorer to average the pinball scores for an array of quantiles for instance.

/cc @GaelVaroquaux @amueller @lorentzenchr

COMMENT:
About names, MAPIE uses the term "prediction intervals" which is nice because it could be used either in a frequentist or a Bayesian context (instead of confidence vs credible intervals) and makes it explicit that this is for a specific prediction.

However "interval" is annoying because we might want to grid the inverse-CDF of the predictive distribution instead of using a pair of quantiles. So maybe we should never speak about intervals in our API (and maybe we should focus on names such as "prediction quantiles" or "grided inverse-CDF of the predictive distribution").

Then we have the problem about when to pass the quantiles: some estimators require them at fit time, and therefore the quantiles should be passed as constructor parameters. Others (e.g. `PoissonRegressor`, MAPIE confirm prediction wrapper) could compute any quantile predictions at predict time without having to specify the grid of quantiles before fitting.

I think we should strive to find an API that can handle both cases.

```python
regressor = FitTimeQuantileAwareRegressor(quantiles=[0.05, 0.5, 0.95])
regressor.fit(X_train, y_train_observations)

y_test_quantiles = regressor.predict(X)

# or with a dedicated method:

y_test_quantiles = regressor.predict_quantiles(X)
```

if the latter, what would `predict` return for this estimator?


```python
regressor = PredictTimeQuantileAwareRegressor()
regressor.fit(X_train, y_train_observations)

y_pred, y_quantiles = regressor.predict(X_test, return_quantiles=[0.05, 0.5, 0.95])

# or with a dedicated method:

y_quantiles = regressor.predict_quantiles(X_test, quantiles=[0.05, 0.5, 0.95])

```

here `y_pred` would typically be the usual expectile.

Labels: API, Needs Decision