How to use the PCE method?

[Einshane]

Hi, I would like to use PCE to perform the uncertainty quantification. I wonder if Chaospy can help me solve the following problem: For example, I have two input parameters subject to the uniform distribution. The output will be calculated by a black box software. I'd like to use the tensor-product expansion. Can chaospy help me figure out the input parameter values that need to do deterministic calculations, i.e. Gaussian quadrature points ? Then I will take these parameters into the software to calculate the output. After obtaining results, can Chaospy combine the quadrature points and output results to calculate the formula of PCE and get the mean, variance of the output and correlation coefficients between the input and output?
Thank you very much.

[jonathf]

If I understand correctly, the answer is yes. It is one of the main usecases of Chaospy.

See here: http://chaospy.readthedocs.io/en/master/tutorial.html#pseudo-spectral-projection-method

To create Gaussian quadrature points, use rule="golub_welsch" instead of "C". Or if you know that you are sticking to Uniform inputs, use the dedicated and somewhat faster Gauss-Legendre method directly with rule="gauss_legendre". It will automatically create tensor-products for multivariate distributions.

[Einshane]

Hi, Jonathf I did some calculations according to your instructions.

nodes=np.array([].[])
evaluations=np.array([])  
distribution=cp.J(cp.Uniform(),cp.Uniform())  
polynomials=cp.orth_ttr(1,distribution)  
proxy_model=cp.fit_regression(polynomials,nodes,evaluations)  
mean=cp.E(proxy_model,distribution)  
variance=cp.Var(proxy_model,distribution)  
print(mean)  
print(variance)

Here I can get the mean and variance of the response. But how can I print the coefficients of the expansion? And I also want to do some sensitivity analysis: local sensitivity analysis (derivatives with respect to expansion variables)、Global sensitivity analysis (the main effect sensitivity index and the total effect index). Can Chaospy calculate these measures and how?

Thank you very much for your help!

[jonathf]

The coefficients from the expansion:

proxy_model, coeffs = cp.fit_regression(polynomials, nodes, evaluation, retall=True)

Main and total sensitivity index:

main = cp.Sens_m(proxy_model, distribution)
total = cp.Sens_t(proxy_model, distribution)

[Einshane]

Thanks for your prompt reply.

But I can only complete global sensitivity analysis after obtaining main and total sensitivity index. Can Chaospy perform local sensitivity analysis?

And I have one more question, I want to try the Smolyak sparse grids method. How can I obtain the sparse nodes and then return the model evaluations to perform the PCE in the above example?

[jonathf]

Local sensitivity analysis? You mean like derivative based analysis? Sure:

approx_model_prime = cp.gradient(approx_model)
sensitivity = approx_model_prime(*location)

Sparse grid is a methodology for going from 1 dimension to multiple dimensions. It is embedded in the cp.generate_quadrature function:

nodes, _ = cp.generate_quadrature(order, distribution, sparse=True, rule="C")

this create Smolyak sparsegrid samples for the Clenshaw-Curtis sampling scheme. CC with apropriate groth rules is fully nested, and therefore a good fit to the Smolyak algorithm.

[Einshane]

Hi, I have another question about the PCE method. polynomials=cp.orth_ttr(1,distribution) generates the total-order expansion. How does Chaospy generate the tensor-product expansion with the dimension preference?

[jonathf]

The first argument order which you have set to 1, can also receive a vector determining the order along each axis. The length of the vector has to be the same as the number of dimensions in the distribution though.

If you also want to stronger or weaker hyperbolic cross-trunction (e.g. prefer cross-terms over non-cross-terms or vice versa), use the flag cross_truncation.

[Einshane]

I change the order into np.array([2,2]) But there is an error: The truth value of an array with more than one element is ambiguous

[jonathf]

You are right there was a bug there.

I've added a patch to the development branch now. If you install it, you should be able to use the feature right away.

[Einshane]

Thank you for your update.

But can I choose to use the total-order expansion or tensor-product expansion? Even if isotropic expansion is used, the tensor expansion and total-order expansion are different, and the coefficients to be calculated are also different. I want to compare the computational efficiency of these two methods.

And I have another question. In the PCE method, I can obtain the statistical moment of the evaluations, but is there any way to estimate the uncertainty band of the evaluations? that is, the lower bound and the upper bound.

[jonathf]

Full tensor-product can be achieved by having a very low cross-truncation order. Somehting like 1e-100 should garanty that you always output the correct expansion I think.

Moments of polynomial expansions do not have errors as they are calculated analytically. The "source of error" comes from the polynomial approximation. To capture that, I propose you look into non-parametric bootstrapping.

[TsubasaHamda]

Hi, I am a student and want to use chaospy for UQ now.

I firstly tried to use this code to create nodes and weights like bellow: import chaospy as cp import numpy as np distribution=cp.J(cp.Uniform(-1,1)) absissas, weights = cp.generate_quadrature(2, distribution, "gauss_legendre")

and outputs are like this: [-1. 0. 1.] [0.16667, 0.66667,0.16667].

But, I want to obtain nodes based on Gaussian quadrature points, so the results should be like this: node [-0.7745966692414833770359, 0, 0.7745966692414833770359] weight [0.5555555555555555555556, 0.8888888888888888888889, 0.555555555555555555556]

How can I get this result and How should I modify my input code? I am sorry to bother you.

[jonathf]

Hello.

No problem. :)

The argument rule in the function generate_quadrature which you are using to select quadrature rule from is not the third positional argument. To change the quadrature rule, you will have to explicitly use the keyword argument:

abscissas, weights = cp.generate_quadrature(2, distribution, rule="gauss_legendre")

One thing to note though, your node is on target, but because we want to use the quadrature nodes in the context of expected values to solve E(X**k) = sum(weights * abcsissas**k), the weights are normalized, giving the values [0.2777, 0.4444, 0.2777]. This would be equivalent to distribution.pdf(abcsissas)*weights from your numbers.

For future reference, please create "new issues" (button near the top right) for a new line of questions.

[TsubasaHamda]

Hi!

Thank you for your reply. I could solve above question.

But, I have another question. If I use total-order expansion( not tensor expansion), how should I change following?

dist1 = cp.Uniform(-1,1) dist2 = cp.Uniform(-1,1) dist = cp.J(dist1,dist2)

nodes,weights = cp.generate_quadrature(2,dist,rule="G")

[jonathf]

Not quite sure what you are asking.

If you are asking about tensor product polynomal expansions vs. total order expansions, this can be achieved by setting the cross_truncation=0.001 to the function orth_ttr. (Moving forward I am adding a patch to allow for the value 0 as well.)

If you are talking about Smolyak sparse grid, then you pass sparse=True to generate_quadrature. But note that sparse grid isn't really that efficient for sampling scheme that isn't properly nested. Gaussian quadrature isn't nested.

[jonathf]

As there is no activity on this thread, I consider the issue as resolved. Feel free to reopen if further discussion is needed.

[ganeshghimire1986]

Hi,

I also tried using Smolyak sparse grid to obtain quadrature nodes for PCE. However, I found it not that accurate in UQ compared to that obtained from tensor product based PCE. I used the regression based approach to compute the PCE coefficients for sparse grid as we do for tensor based scheme. nodes, _ = cp.generate_quadrature(order, distribution, sparse=True, rule="C") proxy_model=cp.fit_regression(polynomials,nodes,evaluations)

Is it the correct approach to compute coefficients for sparse grid considering it's basically the least square method?

[jonathf]

It is not inherently incorrect, but by using regression, you are not really utilizing the sparse grid properties.

To utilize the sparse grid points, you should switch to the fit_quadrature function.