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:
AFAIU, our API is general enough to support multiple outputs from `predict`, same as `predict_proba` for multiclass classifiers.

If a model predicts quantiles, I guess a keyword like `regressor.predict(X_test, quantile_levels=[0.05, 0.5, 0.95]) -> np.array of shape (n_samples, 3)` makes sense to me.

I would not, however, mix quantiles and point estimates of the expectation/mean because if you specify two quantile levels, say 25% and 75%, the estimated expectation might still lie outside this interval. In this case, I would propose to use [expectiles](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expectile.html#scipy.stats.expectile) of different probability levels.

I guess, the problem is then to pass this meta-info ("hi there, here come 3 quantiles for levels 5%, 50% and 95%") to scorers and other model diagnostic tools.

> Classifiers have a `predict_proba` method that makes it possible to quantify probabilistically the certainty in the predictions for a given input `X_i`.

To be precise, `predict_proba` does quantify the likelihood of `Y=class ...` given `X_i` to occur (not the certainty of the prediction itself, also not the certainty of `predict` in general, only if the threshold 50% is applied (which is a design flaw)).

Labels: API, Needs Decision