Sensitivity Analysis function

## Proposal
Add a Sensitivity Analysis (SA) function.

The function would compute _Sobol'_ indices [1,2]. Consider a function `f` with parameters `x1`, `x2` and `x3`. Hence `y=f(x1,x2,x3)`. We are interested to know which parameter has the most impact, in terms of variance, on the value `y`.

I believe scikit-learn, and the wider scientific community, would greatly benefit to have such tool. I believe scikit-learn has something related with `feature_importances_` in some regressors.

As an expert in the field, I would propose to develop the initial functionalities and provide long term support.

## Background
_Sobol’_ indices are the cornerstone of SA and UQ and its application is not restricted to any field. It’s a non intrusive method which makes the only assumption that the variables are independent (this constraint can be alleviated).

Being able to compute sensitivity indices allows to reduce the dimensionality of a problem, better understand the importance of each factors and also see how parameters are interacting with each other. As such, it’s an important engineering tool. If you have 2 or more variables and only have the budget to improve your knowledge on one of them, this can help to make an informed choice. There are a lot of successful usage of SA in the literature and in real world applications. The EU (through the JRC), is now requiring to conduce uncertainty analysis when evaluating a system. They are recommending the use of _Sobol’_ indices.

## Example
Let’s take an example with the _Ishigami_ function:

`y = sin(x1) + 7*sin(x2)^2 + 0.1*x3^4*sin(x1)`

It’s not obvious to tell which variable would impact more `y`.

The Sobol’ indices are bounded from 0 to 1, with 1 meaning more important. Here they would be:

| Variable | First order Sobol’ | Total Sobol’ |
| -------- | ------------------ | ------------ |
|   x1     | 0.31               | 0.56         |
|   x2     | 0.41               | 0.44         |
|   x3     | 0.00               | 0.24         |


The difference between the first and total indices indicate an interaction between variables. The total indices allow to rank the variable by importance. `x1` is the most important. Looking at the first orders, `x3` by itself does not have an impact on the variance of the output. It’s its combination with another variable which makes it have a total impact of 0.24. `x1` in this case also have a difference in first and total indices while `x2` is the same. We can say that `x1` and `x2` have a second order interaction.

Computing the indices requires a large sample size, to alleviate this constraint, a common approach is to construct a surrogate model with Gaussian Process or Polynomial Chaos (to name the most used strategies). This allows to conduct a SA with a low computational budget (we see lots of engineering applications with expensive HPC codes taking advantage of this strategy).

## Implementation

Here is an example implementation I did a while back:

https://gist.github.com/tupui/09f065d6afc923d4c2f5d6d430e11696

Note that this implementation lack a few things such as higher order indices, other methods, input validation, doc, tests, other wrapping, etc. This is just to show how it works.

## Alternatives
Practitioners might be more familiarized with gradient based technics. Keywords include: gradient, adjoint. These methods are _local_ sensitivity analysis methods. And as such they are only capturing the sensitivity of the model in specific areas of the parameter space.

_Sobol’_ indices are variance based indices. There are other indices using higher moments, namely: _moment independant_ based sensitivity analysis. These methods are very attractive and provide lot of information while being simple to compute/analyse. They have been less studied but there is an increasing interest in the community. I believe this would make a good second step.

## References
Upon request, I can provide more information. Following are some fundational references–first 2 cited a few thousand times.

* [1] Sobol',I.M. (2001), _Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates_, MATH COMPUT SIMULAT,55(1–3),271-280, [doi:10.1016/S0378-4754(00)00270-6](http://dx.doi.org/10.1016/S0378-4754(00)00270-6)

* [2] Saltelli, A. et al., (2008), _Global Sensitivity Analysis. The Primer_, John Wiley & Sons, [doi:10.1002/9780470725184](http://dx.doi.org/10.1002/9780470725184)

* [3] Saltelli, A. et al., (2020), _The Future of Sensitivity Analysis: An essential discipline for systems modeling and policy support_, Environmental Modelling & Software, [doi:10.1016/j.envsoft.2020.104954](http://dx.doi.org/10.1016/j.envsoft.2020.104954)

For some examples of industrial applications (Part. 3), a visual explanation of some methods (Chap. 2.2) or more bibliography, you can also see my thesis:

https://github.com/tupui/PHD-Thesis/blob/master/phd_thesis_pamphile_roy.pdf

COMMENT:
For the record, I was too curious and gave it a shot for a simple low-dimensional, non-linear regression problem on numeric-only inputs:


<details>

```python
# %%
from sklearn.inspection import permutation_importance
from sklearn.datasets import make_friedman1, make_friedman2, make_friedman3
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.kernel_approximation import Nystroem
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.polynomial_chaos import PolynomialChaosRegressor
import numpy as np
from scipy.stats import rv_continuous, sobol_indices
import matplotlib.pyplot as plt


model = make_pipeline(Nystroem(kernel="poly", degree=2, n_components=500), Ridge())
model = RandomForestRegressor(n_estimators=100, random_state=0)
model = PolynomialChaosRegressor()

X, y = make_friedman1(n_samples=10_000, random_state=0)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
%time model.fit(X_train, y_train)
print(f"R^2 on test set: {model.score(X_test, y_test):.2f}")


# %%
pi = permutation_importance(
    model, X_test, y_test, n_repeats=10, random_state=0, scoring="r2"
)
normalized_pi = pi.importances_mean / np.sum(pi.importances_mean)


# %%
class EmpiricalCDF(rv_continuous):
    """Marginal distribution defined by the empirical CDF of the data.

    This class was generated by Mistral.
    """

    def __init__(self, data, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.data = np.sort(data)
        self.n = len(self.data)

    def _ppf(self, q):
        # Calculate the percent point function (inverse of CDF)
        indices = np.floor(q * self.n).astype(int)
        indices = np.clip(indices, 0, self.n - 1)
        return self.data[indices]


dists = [EmpiricalCDF(X_train[:, f_idx]) for f_idx in range(X_train.shape[1])]
si = sobol_indices(func=lambda XT: model.predict(XT.T).T, dists=dists, n=2**13)


# %%
feature_names = [f"Feature {i}" for i in range(X.shape[1])]
y_pos = np.arange(len(feature_names))
width = 0.5
fig, ax = plt.subplots(figsize=(10, 6), constrained_layout=True)

importances = [
    ("Normalized PI", normalized_pi),
    ("Sobol first order", si.first_order),
    ("Sobol total order", si.total_order),
]

if hasattr(model, "feature_importances_"):
    importances.append(("RF feature importance", model.feature_importances_))

if hasattr(model, "main_sens") and hasattr(model, "total_sens"):
    importances.append(("PCE Sobol main", model.main_sens().ravel()))
    importances.append(("PCE Sobol total", model.total_sens().ravel()))

for i, (label, values) in enumerate(importances):
    ax.barh(
        y_pos + i * width / len(importances),
        values,
        width / len(importances),
        label=label,
    )

ax.set_yticks(y_pos, labels=feature_names)
ax.invert_yaxis()  # labels read top-to-bottom
ax.set_xlabel("Importance")
ax.set_title(f"Feature Importance: Normalized PI vs Sobol Indices")
_ = ax.legend()

```

</details>

Here are the results for the RF model on a few datasets:


- [`make_friedman1`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman1.html):

<img width="1011" height="611" alt="Image" src="https://github.com/user-attachments/assets/f69c8259-9183-4d5b-9d71-3c37a18b76f9" />

- [`make_friedman2`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman2.html):

<img width="1011" height="611" alt="Image" src="https://github.com/user-attachments/assets/52bd154f-d7ca-4939-a799-40adbfa1ef0a" />

- [`make_friedman3`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman3.html):

<img width="1011" height="611" alt="Image" src="https://github.com/user-attachments/assets/378aa627-7a8d-4645-8e02-ecbb7f143969" />

and the same for the PCE model:

- [`make_friedman1`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman1.html):

<img width="1011" height="611" alt="Image" src="https://github.com/user-attachments/assets/1a340e09-94c7-489d-a471-64ba6eec3177" />

- [`make_friedman2`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman2.html):

<img width="1011" height="611" alt="Image" src="https://github.com/user-attachments/assets/ca797fb2-3d8b-46e2-93e1-8b564b7f0c63" />

- [`make_friedman3`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman3.html):

<img width="1011" height="611" alt="Image" src="https://github.com/user-attachments/assets/b2e45e88-d84f-4986-a9e5-d71e9f0d687a" />

So on those low dim datasets (with independent input features) the results are qualitatively similar as expected.

EDIT: I changed the script to fit an RF model and include the MDI feature importances of the forest.

EDIT 2: I added support to plot the analytical Sobol indices from the Polynomial Chaos Expansion model from #27842. Note that both methods to estimate the Sobol indices on the PCE model match as expected.

Note that RF feature importances (MDI) are computed on the training set only can therefore give strictly positive importance values to uninformative features while this is not the case for PI (computed on a test set) and neither for Sobol indices (which are computed on a synthetic samples).

I suspect that Sobol indices would become untractable on higher dim datasets (or slower to predict models).

Labels: module:inspection, Needs Decision - Include Feature