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:
Hi @ogrisel, this is very nice!

Just to clarify, the PCE model does not provide analytic Sobol' indices in the sense of exact values. What it does make possible is to extract the indices directly from the PCE coefficients, without any additional post-processing. The quality of those indices still depends on the quality of the PCE fit (just as with the Nyström and RF fits from your example).

I'd say that for most problems, the sensitivity indices extracted from a PCE will typically converge faster than the sampling-based alternative (see [here](https://output.circle-artifacts.com/output/job/f561f72e-b283-4f3a-8f7d-9a18135ee414/artifacts/0/doc/auto_examples/polynomial_chaos/plot_pcr_ishigami.html#convergence-of-the-sensitivity-indices) for a comparison between both).

I think it's feasible to derive analytic / exact values for the Friedman 1 problem, but it might be too hard for the Friedman 2 and Friedman 3 problem. Using `scipy.stats.sobol_indices`, though, the numerical estimates converge nicely:

- **Friedman 1**
```
Main indices: [0.19731097 0.19731098 0.09326697 0.34975116 0.08743778]
Total indices: [0.27223402 0.27223402 0.09326698 0.34975116 0.08743779]
```
<img width="1578" height="598" alt="Image" src="https://github.com/user-attachments/assets/35188282-0e4b-4cd2-a52d-8f25a7619b2e" />

- **Friedman 2**
```
Main indices: [ 3.15059084e-04  3.74631419e-01  4.92472466e-01 -6.77268968e-12]
Total indices: [6.26005939e-04 5.06992786e-01 6.24975676e-01 7.74579717e-13]
```
<img width="1577" height="598" alt="Image" src="https://github.com/user-attachments/assets/253fb3c1-7f5a-4b05-ae11-effdf5add8cc" />

- **Friedman 3**
```
Main indices: [1.27094912e-01 1.67153099e-01 5.76958325e-01 3.88345491e-09]
Total indices: [2.22827823e-01 2.37028824e-01 6.79756354e-01 2.29890554e-10]
```
<img width="1577" height="598" alt="Image" src="https://github.com/user-attachments/assets/236cdac4-b15c-4b3b-948d-9d61aaca4059" />

<details>
<summary>Source</summary>

```python
# %%
# Import statements
import matplotlib.pyplot as plt

import numpy as np
from scipy.stats import sobol_indices, uniform

# %%
# Default rng
rng = np.random.default_rng()

# %%
# Range of log2 number of points to test
m_range = range(4, 21)

# %%
# Function that evaluates main and total Sobol' indices for a given range of number of points
def sobol_indices_sweep(func, dists, m_range):
    main = []
    total = []
    for m in m_range:
        indices = sobol_indices(func=func, n=2**m, dists=dists, rng=rng)
        main.append(indices.first_order)
        total.append(indices.total_order)
    main = np.vstack(main)
    total = np.vstack(total)

    return main, total

# %%
# Function that plots the convergence of the Sobol' indices
def plot_indices(problem, main, total):
    _, axes = plt.subplots(1, 2, figsize=(8, 3))
    for ax, indices, label in zip(axes.flatten(), (main, total), ("main", "total")):
        for k, index in enumerate(indices.T):
            ax.plot(m_range, index, label=fr"$x_{k}$")
        ax.set_xlabel("log2 number of points")
        ax.set_ylabel(f"{label}-effect Sobol' indices")
        if label == "total":
            ax.legend(frameon=False, bbox_to_anchor=(1, 1), loc="upper left")
    plt.suptitle(problem)
    plt.tight_layout()
    plt.show()

# %%
# Friedman 1
friedman_1 = lambda x: 10*np.sin(np.pi*x[0]*x[1]) + 20*(x[2] - 0.5)**2 + 10*x[3] + 5*x[4]
dists = [uniform()]*5
main, total = sobol_indices_sweep(friedman_1, dists, m_range)

print("Main indices:", main[-1])
print("Total indices:", total[-1])

plot_indices("Friedman 1", main, total)

# %%
# Friedman 2
friedman_2 = lambda x: np.sqrt(x[0]**2 + (x[1]*x[2] - 1/(x[1]*x[3]))**2)
dists = [uniform(0, 100), uniform(40*np.pi, 520*np.pi), uniform(), uniform(1, 10)]
main, total = sobol_indices_sweep(friedman_2, dists, m_range)

print("Main indices:", main[-1])
print("Total indices:", total[-1])

plot_indices("Friedman 2", main, total)

# %%
# Friedman 3
friedman_3 = lambda x: np.arctan((x[1] * x[2] - 1 / (x[1] * x[3])) / x[0])
dists = [uniform(0, 100), uniform(40*np.pi, 520*np.pi), uniform(), uniform(1, 10)]
main, total = sobol_indices_sweep(friedman_3, dists, m_range)

print("Main indices:", main[-1])
print("Total indices:", total[-1])

plot_indices("Friedman 3", main, total)

```

</details>

I think something might be off for the total sensitivity indices in your plot for the Friedman 1 problem using the PCE model. Are you reporting the main sensitivity indices twice by any chance? They don't match with the values reported for the Friedman 1 problem with the RF model.

It's also very interesting to see how the RF feature importance reports non-zero values for features 5 to 9 for the Friedman 1 problem. They should be exactly 0 as the problem does not depend on these features.

For the blog post, it would be very nice if we could find an example where the feature importance or permutation importance and the Sobol' indices do not agree. I'll try to think of something.

EDIT: For completeness, here's how you could compute the Sobol' sensitivity indices analytically using `sympy`. This seems to work for Friedman 1, but it took too long on my computer for Friedman 2 and 3 (I gave up after 10 minutes).

<details>
<summary>Details</summary>

```python
import sympy as sp

def sobol_indices(f_expr, vars, bounds, verbose=False, return_total=True, mean=None):
    """
    Compute analytical expressions for mean, variance, and Sobol indices
    using multi-dimensional integration in one go (no for-loops).
    """
    d = len(vars)
    pdfs = [sp.S.One / (b - a) for (a, b) in bounds]
    density = sp.Mul(*pdfs)
    
    if verbose:
        print("PDFs:", pdfs)
        print("Joint density:", density)
    
    # Define integration ranges once
    ranges = [(vars[i], bounds[i][0], bounds[i][1]) for i in range(d)]
    
    # Mean: E[f]
    if mean is None:
        mean = sp.integrate(f_expr * density, *ranges)
        if verbose:
            print("Mean:", mean)
    
    # Second moment: E[f^2]
    second_moment = sp.integrate(f_expr**2 * density, *ranges)
    if verbose:
        print("Second moment:", second_moment)
    
    variance = second_moment - mean**2
    if verbose:
        print("Variance:", variance)
    
    Si = {}
    STi = {}
    
    for i, xi in enumerate(vars):
        if verbose:
            print(f"\n--- Variable {xi} ---")
        
        # First-order: Var(E[f|Xi])
        cond_ranges = [(vars[j], bounds[j][0], bounds[j][1]) for j in range(d) if j != i]
        cond_exp = sp.integrate(f_expr, *cond_ranges) / sp.Mul(*[bounds[j][1]-bounds[j][0] for j in range(d) if j != i])
        var_cond = sp.integrate(cond_exp**2 * pdfs[i], (xi, bounds[i][0], bounds[i][1])) - mean**2
        Si[xi] = var_cond / variance
        
        if verbose:
            print(f"First-order numerator Var(E[f|{xi}]): {var_cond}")
            print(f"First-order index Si[{xi}]: {Si[xi]}")
        
        # Total-effect: E[Var(f|X_-i)] / Var
        if return_total:
            cond_var = f_expr**2 - sp.integrate(f_expr, (xi, bounds[i][0], bounds[i][1]))**2
            cond_var = sp.integrate(cond_var, *cond_ranges) / sp.Mul(*[bounds[j][1]-bounds[j][0] for j in range(d) if j != i])
            STi[xi] = sp.integrate(cond_var * pdfs[i], (xi, bounds[i][0], bounds[i][1])) / variance
            if verbose:
                print(f"Total-effect numerator E[Var(f|X_-{xi})]: {STi[xi]}")
                print(f"Total-effect index STi[{xi}]: {STi[xi]}")
    
    result = {"mean": mean, "variance": variance, "Si": Si}
    if return_total:
        result["STi"] = STi
    
    return result

# Friedman #1
x1, x2, x3, x4, x5 = sp.symbols("x1:6", real=True, positive=True)
f = 10*sp.sin(sp.pi*x1*x2) + 20*(x3 - sp.Rational(1,2))**2 + 10*x4 + 5*x5
bounds = [(0, 1)]*5

results = sobol_indices(f, [x1, x2, x3, x4, x5], bounds)

# Friedman #2
x1, x2, x3, x4 = sp.symbols("x1:5", real=True, positive=True)
f = sp.sqrt(x1**2 + (x2*x3 - 1/(x2*x4))**2)
bounds = [(0, 100), (40*sp.pi, 560*sp.pi), (0, 1), (1, 11)]

results = sobol_indices(f, [x1, x2, x3, x4], bounds)

# Friedman #3
x1, x2, x3, x4 = sp.symbols("x1:5", real=True, positive=True)
f = sp.atan((x2*x3 - 1/(x2*x4))/x1)
bounds = [(0, 100), (40*sp.pi, 560*sp.pi), (0, 1), (1, 11)]

results = sobol_indices(f, [x1, x2, x3, x4], bounds)

```
</details>

Also, I saw your review comments in [27842](https://github.com/scikit-learn/scikit-learn/pull/27842), thanks! I’ll take a look later this week.

Labels: module:inspection, Needs Decision - Include Feature