# Mauntz-Kucherenko formula **URL:** \ **Category:** Development\ **Tags:** sensitivity-analysis\ **Created:** [December 1, 2020, 4:13pm UTC](https://openturns.discourse.group/t/mauntz-kucherenko-formula/68 "2020-12-01T16:13:19Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\ **Post date:** [December 1, 2020, 4:13pm UTC](https://openturns.discourse.group/t/mauntz-kucherenko-formula/68/1 "2020-12-01T16:13:19Z")
This post follows the discussion started on Github: > > > If I understand the code correctly, the implemented formula for the MauntzKucherenko V\_{-i} estimator is: > where $\\tilde{G}$ is the "centered" model, i.e.... In the following, I use the notations defined on this page of the doc: > **[SobolIndicesAlgorithm — OpenTURNS 1.16 documentation](https://openturns.github.io/openturns/1.16/user_manual/_generated/openturns.SobolIndicesAlgorithm.html)** Applying Equation (12) of the original paper > **[sobol2007.pdf](https://openturns.github.io/openturns/papers/sobol2007.pdf)** > > 213.93 KB the Mauntz-Kucherenko estimator of the first order Sobol index is: \hat{V}\_i = \frac{1}{N} \sum\_{k=1}^N \tilde{G}(\boldsymbol{B}\_k) \left[\tilde{G}(\boldsymbol{E}\_k^{i}) - \tilde{G}(\boldsymbol{A}\_k) \right]. Logically, the estimator \hat{V}\_{-i} should then be \hat{V}\_{-i} = \frac{1}{N} \sum\_{k=1}^N \tilde{G}(\boldsymbol{B}\_k) \left[\tilde{G}(\boldsymbol{C}\_k^{i}) - \tilde{G}(\boldsymbol{A}\_k) \right]. If we cannot do that (for example because we have not implemented \boldsymbol{C}\_k^{i} 😉 ), then we could define \hat{V}\_{-i}^{alt} = \frac{1}{N} \sum\_{k=1}^N \tilde{G}(\boldsymbol{A}\_k) \left[\tilde{G}(\boldsymbol{E}\_k^{i}) - \tilde{G}(\boldsymbol{B}\_k) \right]. and that would remain consistent with the paper. In that case, we would define \widehat{VT}\_{i}^{alt} as \frac{1}{N} \sum\_{k=1}^N \tilde{G}(\mathcal{A}\_k)^2 - \hat{V}\_{-i}^{alt} (keep in mind that \sum\_{k=1}^N \tilde{G}(\boldsymbol{A}\_k) is supposed to be null). Therefore \widehat{VT}\_{i}^{alt} = \frac{1}{N} \sum\_{k=1}^N \tilde{G}(\boldsymbol{A}\_k) \left[\tilde{G}(\boldsymbol{A}\_k) + \tilde{G}(\boldsymbol{B}\_k) - \tilde{G}(\boldsymbol{E}\_k^{i}) \right]. Of course, in OpenTURNS, we typically use unbiased variance estimators so \frac{1}{N} should be replaced in all formulas by \frac{1}{N-1}. Not that it matters with the large N involved in Sobol estimation. ---
**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\ **Post date:** [December 2, 2020, 10:10am UTC](https://openturns.discourse.group/t/mauntz-kucherenko-formula/68/2 "2020-12-02T10:10:28Z")
Here is the formula for \widehat{VT}\_i^{alt} actually implemented in the code: > `VTi(q, p) = referenceVariance_[q] + (size * muA[q] * muA[q] - yEDotyA[q]) / (size - 1.0);` In the code, `q` is the output dimension. For the sake of simplicity, let us assume that the output is one-dimensional, so `q` is `0`. The input dimension (denoted by i in the post above) is here denoted by `p`: `p` =i. Currently, \tilde{G} is centered with respect to the full (G(\boldsymbol{A}), G(\boldsymbol{B}), G(\boldsymbol{E}^1),...,G(\boldsymbol{E}^{n\_X}))^T: > # Case 1 : \tilde{G} centerered with respect to the sample \boldsymbol{A} In this case, we assume that the code snippet above is changed in order to center the output sample with respect to \boldsymbol{A}. That is to say that \tilde{G}(\cdot) = G(\cdot) - \frac{1}{N} \sum\_{k=1}^N G(\boldsymbol{A}\_k). In this case we would have: `referenceVariance_[q]` = \frac{1}{N-1} \sum\_{k=1}^{N} \tilde{G}(\boldsymbol{A}\_k)^2, `size` = N, `muA[q]` = \frac{1}{N} \sum\_{k=1}^N \tilde{G}(\boldsymbol{A}\_k) = 0, `yEDotyA[q]` = \sum\_{k=1}^N \tilde{G}(\boldsymbol{E}\_k^i) \tilde{G}(\boldsymbol{A}\_k). Therefore, `VTi(q, p)` = \frac{1}{N-1} \sum\_{k=1}^{N} \tilde{G}(\boldsymbol{A}\_k) \left( \tilde{G}(\boldsymbol{A}\_k) - \tilde{G}(\boldsymbol{E}\_k^i) \right) # Case 2 : \tilde{G} centerered with respect to the full sample This is what is currently implemented, so the formulas below are what is really computed in OT version 1.16. With \tilde{G}(\cdot) = G(\cdot) - \frac{1}{(2+n\_X)N} \left( \sum\_{k=1}^N G(\boldsymbol{A}\_k) + \sum\_{k=1}^N G(\boldsymbol{B}\_k) + \sum\_{i=1}^{n\_X} \sum\_{k=1}^N G(\boldsymbol{E}\_k^i) \right), we have `referenceVariance_[q]` = \frac{1}{N-1} \sum\_{k=1}^{N} \left( \tilde{G}(\boldsymbol{A}\_k) - \frac{1}{N} \sum\_{k'=1}^N \tilde{G}(\boldsymbol{A}\_{k'})\right)^2, `size` = N, `muA[q]` = \frac{1}{N} \sum\_{k=1}^N \tilde{G}(\boldsymbol{A}\_k) \neq 0, `yEDotyA[q]` = \sum\_{k=1}^N \tilde{G}(\boldsymbol{E}\_k^i) \tilde{G}(\boldsymbol{A}\_k). Therefore, `VTi(q, p)` = \frac{1}{N-1} \sum\_{k=1}^{N} \left( \tilde{G}(\boldsymbol{A}\_k) - \frac{1}{N} \sum\_{k'=1}^N \tilde{G}(\boldsymbol{A}\_{k'})\right)^2 + \frac{N}{N-1}\left( \frac{1}{N} \sum\_{k=1}^{N} \tilde{G}(\boldsymbol{A}\_k) \right)^2 - \frac{1}{N-1} \sum\_{k=1}^N \tilde{G}(\boldsymbol{E}\_k^i) \tilde{G}(\boldsymbol{A}\_k). # Conclusion The formulas of both Case 1 and Case 2 do not fit the theory. The doc is wrong too, but in a different way. ---
**Author:** ![sofianehaddad](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/sofianehaddad/32/46_2.png) [@sofianehaddad](https://openturns.discourse.group/u/sofianehaddad)\ **Post date:** [December 2, 2020, 1:45pm UTC](https://openturns.discourse.group/t/mauntz-kucherenko-formula/68/3 "2020-12-02T13:45:33Z")
@josephmure Using what your wrote and analysing a bit theory/other sources, we can refine the analysis. Starting from the estimate of first order: \hat{V}\_i=\frac{1}{N}\sum\_{k=1}^N \tilde{G(\boldsymbol{B}\_k)}\left[\tilde{G(\boldsymbol{E}\_k^{i})} - \tilde{G(\boldsymbol{A}\_k)}\right] We get: \hat{V}\_{-i}=\frac{1}{N}\sum\_{k=1}^N \tilde{G(\boldsymbol{B}\_k)}\left[\tilde{G(\boldsymbol{C}\_k^{i}}) - \tilde{G(\boldsymbol{A}\_k)}\right] As C\_k are missing, we use (as you mentioned): \widehat{V}\_{-i}^{alt}=\frac{1}{N}\sum\_{k=1}^N\tilde{G(\boldsymbol{A}\_k)}\left[\tilde{G(\boldsymbol{E}\_k^{i})} - \tilde{G(\boldsymbol{B}\_k)}\right]. Then to get \widehat{VT}\_{i}^{alt}, we use: \widehat{VT}\_{i}^{alt}= Var( \tilde{G(\boldsymbol{A}\_k)} ) - \widehat{V}\_{-i}^{alt} The tip here consists, as for first order, to replace \left(\frac{1}{N}\sum\_{k=1}^N \tilde{G(\boldsymbol{A}\_k)}\right)^2 by \frac{1}{N}\sum\_{k=1}^N \tilde{G(\boldsymbol{A}\_k)} \tilde{G(\boldsymbol{B}\_k)}. Thus we get: \widehat{VT}\_{i}^{alt}= \frac{1}{N}\sum\_{k=1}^N\tilde{G(\boldsymbol{A}\_k)}\left[(\tilde{G(\boldsymbol{A}\_k}) - \tilde{G(\boldsymbol{E}\_k^{i})})\right] If we have a look for example to `R-sensitivity`, the estimator is implemented like that. However, in OpenTURNS, we replace the first term \frac{1}{N}\sum\_{k=1}^N\tilde{G(\boldsymbol{A}\_k)}(\tilde{G(\boldsymbol{A}\_k}) by var(y\_a) + \mu\_a^2. In other words, we should implement the good one with the previous equation