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:
Personally, it would be great to see the scikit package take some direction on what a future multi quantile api might look like. Scikit has such a great influence it is basically sets the standard for python ml best practice.

Here are a few packages with differing APIs based on scikit learn:
https://github.com/StatMixedML/XGBoostLSS/blob/master/xgboostlss/model.py#L468C47-L468C47
https://github.com/StatMixedML/LightGBMLSS/blob/master/lightgbmlss/model.py#L447
https://xgboost.readthedocs.io/en/latest/parameter.html#parameter-for-using-quantile-loss-reg-quantileerror
https://catboost.ai/en/docs/concepts/loss-functions-regression#MultiQuantile

Other well known packages using quantiles such as MAPIE and LightGBM often resort looping over a list of quantiles. [NGBoost](https://github.com/stanfordmlgroup/ngboost), [PGBM](https://github.com/elephaint/pgbm) and [XGBoost-Distribution](https://github.com/CDonnerer/xgboost-distribution) could all easily add multi-quantile predictions to their scikit APIs - my guess is the main reason they haven't is that there is not standard api coming from scikit on how best to do this!

I believe a `.predict_quantile()` method would be very nice but equally something like XGBoostLSS and LightGBMLSS's `predict(X, pred_type="quantiles")` works nicely.

I would love to hear your in put on this as I am looking to looking to work on some of these projects and think it would be good to conform on an somewhat standard Multi Quantile API.

Labels: API, Needs Decision