Test

The integrals defining Sobol indices cannot be evaluated analytically for most computational models. Monte Carlo estimation provides a practical approach. The first-order index can be expressed as
$$S_i=\frac{{\text{Var}}_{{\theta }_i}\left[{\mathbb{E}}_{{\mathbf{\theta }}_{\sim i}}(Y|{\theta }_i)\right]}{\text{Var}(Y)},$$
where ${\mathbb{E}}_{{\mathbf{\theta }}_{\sim i}}(Y|{\theta }_i)$ denotes the expected output when ${\theta }_i$ is fixed and all other parameters vary. The outer variance operator measures how this conditional expectation varies as ${\theta }_i$ changes. Similarly, the total-effect index is
$$S_{T_i}=\frac{{\mathbb{E}}_{{\mathbf{\theta }}_{\sim i}}\left[{\text{Var}}_{{\theta }_i}(Y|{\mathbf{\theta }}_{\sim i})\right]}{\text{Var}(Y)},$$
representing the expected variance remaining when all parameters except ${\theta }_i$ are fixed.

Direct Monte Carlo evaluation of these conditional expectations requires nested sampling loops, which rapidly becomes computationally prohibitive. Saltelli proposed an efficient sampling scheme that circumvents nested loops. The method generates two independent sample matrices $\mathbf{A}$ and $\mathbf{B}$, each of size $N\times d$, where each row represents one parameter set sampled from the joint distribution. For each parameter $i$, a hybrid matrix ${\mathbf{AB}}_i$ is constructed by copying $\mathbf{A}$ and replacing its $i$-th column with the corresponding column from $\mathbf{B}$. The model is evaluated at all sample points to obtain response vectors ${\mathbf{Y}}_A$, ${\mathbf{Y}}_B$, and ${\mathbf{Y}}_{A{B}_i}$ for $i=1,\dots ,d$.

Using these evaluations, the first-order index is estimated as
$$S_i=\frac{\frac{1}{N}{\sum }_{j=1}^N{Y}_B^{(j)}\left(Y_{A{B}_i}^{(j)}-Y_A^{(j)}\right)}{\text{Var}(Y)},$$
and the total-effect index as
$$S_{T_i}=\frac{\frac{1}{2N}{\sum }_{j=1}^N{\left(Y_A^{(j)}-Y_{A{B}_i}^{(j)}\right)}^2}{\text{Var}(Y)}.$$
Computing indices for all $d$ parameters requires $N(d+2)$ model evaluations, see \citet{Saltelli2002MakingBest}. When $N$ is chosen to be several thousand, this becomes expensive for models requiring substantial computational time per evaluation. Surrogate models that approximate the original model with negligible evaluation cost make Sobol analysis tractable for such cases, \citep{Saltelli2008GlobalSA}.

Many computational models produce time-dependent or spatially distributed outputs rather than scalar responses. A dynamic simulation generates a time series $Y(t)$ for $t\in [0,T]$. A spatial model produces a field $Y(\mathbf{x})$ over a domain. The standard Sobol indices defined for scalar outputs do not directly apply.

The most straightforward extension computes Sobol indices pointwise. For a time series, indices are calculated separately at each time instant:
$$S_i(t)=\frac{{\text{Var}}_{{\theta }_i}\left[{\mathbb{E}}_{{\mathbf{\theta }}_{\sim i}}(Y(t)|{\theta }_i)\right]}{\text{Var}(Y(t))}.$$
This produces time-varying sensitivity profiles revealing which parameters dominate at different stages of the process. Early-time sensitivity may differ from late-time behavior. Pointwise indices are computed independently at each time step using the same sample matrices $\mathbf{A}$, $\mathbf{B}$, and ${\mathbf{AB}}_i$, with model evaluations providing response time series rather than scalars.

For applications requiring a single global importance measure across time, pointwise indices can be aggregated. Common approaches include averaging indices over time, integrating them against a weight function that emphasizes critical time windows, or extracting maximum values to identify peak sensitivity. The choice depends on whether early-time, late-time, or transient behavior matters most for the application. In cases where the entire time evolution matters equally, simple averaging provides a reasonable summary.