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:
> However, there seems to be a constant factor problem that would need to be investigated.

I think that's because the paper compares "unnormalized" total Sobol' sensitivity indices to (conditional) permutation importance values. The normalization factor here is the variance. In our example, this is

$$Var(Y) = \mathbb{E}[X_1^2]\mathbb[{E}[X_2^2] = 1/3 \times 1/3 = 1/9.$$

I get the same values if I rescale all Sobol' sensitivity indices by the variance and plot raw permutation importance values instead of normalized ones. So I think the relation we're looking for is

$$PI_i = 2 \frac{S_{T_i}}{V}$$

where $V$ is the variance of the output.

<img width="2023" height="1223" alt="Image" src="https://github.com/user-attachments/assets/8b8e7153-b6d4-457f-954b-be82807ae6e3" />

<details>

# %%
# Import statements
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split
from sklearn.polynomial_chaos import PolynomialChaosRegressor

import itertools
import numpy as np
from scipy.stats import sobol_indices, uniform
import matplotlib.pyplot as plt

# %%
# Set seed
seed = 123
np.random.seed(seed)

# %%
# Generate training data
n_sobol = 4096
n_samples = 4 * n_sobol
distributions = [uniform(loc=-1, scale=2)]*2
X = np.vstack([distribution.rvs(n_samples) for distribution in distributions]).T
y = X[:, 0] * X[:, 1]
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

# %%
# Use a sparse solver
model = PolynomialChaosRegressor(distribution=[uniform(loc=-1, scale=2)]*2)

# %%
# Fit the model
model.fit(X_train, y_train)
print(f"R^2 on test set: {model.score(X_test, y_test):.2f}")

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

# %%
# Compute Sobol indices (from model)
si = sobol_indices(func=lambda x: x[0] * x[1], dists=distributions, n=n_sobol, rng=seed)

# %%
# Assume model has d input features
d = X.shape[1]
feature_names = [f"x{i+1}" for i in range(d)]

# Collect all main effects
indices_first_pce = []
labels = []
for i in range(d):
    indices_first_pce.append(model.joint_sens(i) * model.var())
    labels.append(feature_names[i])

# Collect all joint effects (pairwise and higher order if desired)
for r in range(2, d+1):
    for combo in itertools.combinations(range(d), r):
        indices_first_pce.append(model.joint_sens(*combo) * model.var())
        labels.append(", ".join(feature_names[j] for j in combo))

indices_first_pce = np.array(indices_first_pce).ravel()

indices_total_pce = np.zeros(len(labels))
indices_total_pce[:d] = model.total_sens() * model.var()

# Pad external importances (only for main effects)
normalized_pi_full = np.zeros(len(labels))
normalized_pi_full[:d] = normalized_pi  # assumes normalized_pi is length d

sobol_first_full = np.zeros(len(labels))
sobol_first_full[:d] = si.first_order * np.var(y) # assumes si.first_order is length d

sobol_total_full = np.zeros(len(labels))
sobol_total_full[:d] = si.total_order * np.var(y) # assumes si.total_order is length d

# Collect all importances
importances = [
    ("1/2 * PI", 0.5*normalized_pi_full),
    ("Unnormalized Sobol first order", sobol_first_full),
    ("Unnormalized Sobol total order", sobol_total_full),
    ("PCE Unnormalized Sobol first order", indices_first_pce),
    ("PCE Unnormalized Sobol total order", indices_total_pce),
]

# Plot
y_pos = np.arange(len(labels))
width = 0.5
fig, ax = plt.subplots(figsize=(10, 6), constrained_layout=True)

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

ax.set_yticks(y_pos, labels=labels)
ax.invert_yaxis()
ax.set_xlabel("Importance")
ax.set_title("Feature Importance: PI vs Sobol (sampling vs PCE)")
ax.legend(loc="upper left", bbox_to_anchor=(1, 1), frameon=False)
plt.show()

</details>

Labels: module:inspection, Needs Decision - Include Feature