Univariate effect

[takafusui]

Dear Joanthan,

Thank you very much at first implementing a nice and reliable tool. I started to learn uncertainty quantification recently and accessed to chaospy through the Uncertainpy package.

I am interested in computing a univariate effect, following Eqs. (31) and (32) in the paper below: Deman et al. (2018), Using sparse polynomial chaos expansions for the global sensitivity analysis of groundwater lifetime expectancy in a multi-layered hydrogeological model, DOI: 10.1016/j.ress.2015.11.005

According to the paper, we can compute a univariate effect directly after post-processing the coefficients of PCE, which is similar to the Sobol' indexes.

Is this function already implemented in chaospy? I think no. Is it easy to implement by myself?

Thank you in advance for your help and time.

Best regards, Takafumi

[jonathf]

My time is a bit limited, so I've just skimmed the paper. I apologies if I misunderstood something. I fond your eq (31), which indeed looks like Sobol indices. I will assume that is what you are after.

Assuming you are already familiar with chaospy, the following script I just made will give you an idea of how to calculate the indices as described in (31).

import numpy
import chaospy

def foo(q):
    """Your solver here."""
    return q[0]*numpy.e**-q[1]

dist = chaospy.Iid(chaospy.Uniform(-1, 1), 3) # or whatever
poly_order = 5
phi = chaospy.orth_ttr(poly_order, dist)

# collocation method:
samples = dist.sample(500)
evals = [foo(sample) for sample in samples.T]
approx_model, y = chaospy.fit_regression(phi, samples, evals, retall=True)

# Gather some indices matching that of phi
alpha = numpy.array(chaospy.bertran.bindex(start=0, stop=poly_order, dim=3, sort="GR"))
i = 1

# calculating (31):
indices = (alpha[:, i] > 0) & (numpy.sum(alpha, -1) == alpha[:, i])
Mi_1 = sum(y[indices]*phi[indices])
# remove (non-existant) other variables from equation:
Mi_1 = Mi_1(q0=0, q2=0)
print(Mi_1)

In addition, I do have tools for calculate Sobol indices directly:

chaospy.Sens_m   # 1st order main
chaospy.Sens_m2  # 2nd order main
chaospy.Sens_t   # total

Neither of these uses the formula in your paper though, but in principle they should give you the same answers.

[takafusui]

Thank you very much Jonathan for your very detailed explanation. I am a novice user of chaospy, thus, I need some time to study chaospy and your paper. Sorry for my late replay.

I am still studying chaospy and please let me clarify whether I understood your code correctly or not.

First, dist = chaospy.Iid(chaospy.Uniform(-1, 1), 3) # or whatever Why do you define three uncertain parameters, instead of two, though foo requires two parameter arguments?

Second, I do not understand the below part: # Gather some indices matching that of phi alpha = numpy.array(chaospy.bertran.bindex(start=0, stop=poly_order, dim=3, sort="GR")) i = 1 # calculating (31): indices = (alpha[:, i] > 0) & (numpy.sum(alpha, -1) == alpha[:, i]) You look like you correct indices that satisfy Eq. 32 of the paper. Is it correct? Please note that the set A_{i} is as same as that of the first-order Sobol' index (please see Eq. 28 in the paper). I also do not understand the argument "GR". What does "GR" stand for?

Third, len(Mi_1) = poly_order is strange for me. I think len(Mi_1) should be 1 since Eq. 32 is PCE. You take Mi_1 = sum(y[indices]*phi[indices]) but still the length is not 1. I am confusing this. I think we should evaluate PCE something similar manner to Mi_1 = approx_model(q1)

I am sorry to bother you but I really appreciate your help on this matter. Thank you.

[jonathf]

1. I have a fast trigger finger and do faster than I think. Sometime that means inconsistencies like 2 vs. 3 dimensions.

2. Yes, that should be the equation, but it says (31) in my version. I guess I might have found a different revision or something.

3. The sort argument references the order of the polynomials. GR implies graded and reverse. See the docs of chaospy.basis and you will see what that means.

4. Yeah, Mi_1 might need a second sum (Mi_1 = sum(Mi_1)) to get the scalar polynomials. Again, the speed of my typing gets ahead of me thought.