Errors increasing with polynomial order and calculation of condition mean

[Sumit9013]

Dr. Feinberg,

My procedure is defined below:
I have 12 input variables X = Xi (i = 1,2,..,12) each normally distrbuted.
I have generated random samples and ran vibration simulation using abaqus on generated samples and got 250 evaluation points.
I have then generated 2nd order othogonal expansion and calculated the coefficients using the regression.
approx_function = chaospy.fit_regression(expansion, samples, evaluations) these have points using the data generated by the vibration simulations of 250 samples
using this i have got a approx_function of shape (90,). I have then calculated the sobol indices.
I am checking the error in Mean and variance by:
True mean (mean of 250 evaluation points) - Approx Mean (mean using approx_function and distribution)
True variance (variance of 250 evaluation points) - Approx variance (variance using approx_function and distribution)

My doubts are as follows:
1)Why are errors increasing on increasing the order of polynomial? How to determine which is the best order to be used?
2) How do i calculate Conditional Mean and Sobol indices (without using m_sens function):
E(Y/Xi) = chaospy.E_cond(approx_function, freeze, joint_distribution) because freeze is accepted as a polynomial and not the RV distribution.
Sxi = Var(E(Y/Xi)(This is calculated over distribution of Xi)/Var(Y) (this is calculated over the joint distribution)

[jonathf]

1. Are you increaing the polynomial order without increasing the number of evaluated samples? For regression, try to keep the number of evaluated samples roughtly 2 times the number of parameters solving for.
2. Get your manuel calculations like this:

variable = chaospy.variable(12)
e_y_xi = chaospy.E_cond(approx_function, variable[i], joint_distribution)
sxi = chaospy.Var(e_y_xi, joint_distribution)/chaospy.Var(approx_function, joint_distribution)

Hope this helps.

[Sumit9013]

Dr Feinberg,

Thank you for replying.

1. Yes, i am not increasing the sample points. I am determining the order by comparing the number of co-efficients (that the order gives) vs the sample points such that the former is lesser than the latter. I think you meant that the sample points >= 2 *number of deterministic co-efficients. If not, I would appreciate if you could elaborate on the order number.
2. Thank you I will try it.

[jonathf]

1. yes, that is what I meant.

[Sumit9013]

hello Dr Feinberg,
I have an issue with consistency in the module. My code is given below. I am extracting the main effects of each variable.

variable = chaospy.variable(12)
mean = chaospy.E(approx_function, joint_distribution)
for i in range(len(variable)):
main_effects = chaospy.E_cond(approx_function, variable[i], joint_distribution) - mean

In the above code while debugging i found out that while variable[i] is q2 the main effects is a function of q4 and the condition mean function doesnt freeze the variables correctly. Could you please hint where i am possibly going wrong.

[jonathf]

I am unsure. Could you post the following:

• OS or Distro
• Version of Python
• Version of Chaospy and Numpoly

From your example, do mind posting the output from:

• approx_function.names
• joint_distribution

[Sumit9013]

Here are the information you asked for:
OS : Windows 10 (64 bit)
Version of python : 3.7
Version of chaospy : 4.3.2
Version of numpoly : 1.2.3
approx_function.names
Out[213]: ('q0', 'q1', 'q10', 'q11', 'q2', 'q3', 'q4', 'q5', 'q6', 'q7', 'q8', 'q9')
approx_function
Out[214]: polynomial(3.2571009070176376e-05q9+0.0006414454772852574q8+0.0004933292599055015q7-0.0007797560917218664q6+0.000652856879873117q5+3.6371167984546095e-05q4-2.470233162860938e-05q3+4.199309431183505e-05q2+0.00083850924459076q11-3.0040556898121853e-05q10-0.0007319851968323927q1+6.259819993009027e-06q0-0.024978669449116083)
variable
Out[215]: polynomial([q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11])
JOINT_DIST
Out[216]: J(Normal(mu=4, sigma=1.8), Normal(mu=9.2, sigma=5.3), Normal(mu=3.75, sigma=4.132), Normal(mu=13.332, sigma=2.3112), Normal(mu=5.12, sigma=5.112), Normal(mu=14.2, sigma=5.2), Normal(mu=32, sigma=6.32), Normal(mu=29.53, sigma=7.432), Normal(mu=2, sigma=0.89), Normal(mu=22, sigma=11.32), Normal(mu=2, sigma=0.95), Normal(mu=11, sigma=3.9))
print(cp.E_cond(approx_function,variable[2],JOINT_DIST))
Out[217]: 3.6371167984546095e-05*q4+0.009456457965737242

as shown above variable[2] is q2 and conditional mean is a function of q4

[jonathf]

Ah, the variable names are not sorted for more than 10 variables. That makes sense.

I've made a new release of Chaospy now version 4.3.7. which has a fix for conditional expectations. Let me know if that solves your problem.

[Sumit9013]

Are the sobol indices using Sens_m function also dependent on the sorting of variables in the approx_function? Do they give in the same order as the variables in joint distribution?

[jonathf]

In principle yes as Sens_m is using E_cond underneath the hood. But I don't think so from skimming the code.
Are you getting different results with Sens_m in 4.3.7?

[Sumit9013]

Thanks i have upgraded to 4.3.7.
Whenever i am changing the sequence of input variables in the joint distribution and the samples, the sobol indices are changing for the same RV. For instance, originally, the SI of q5 was 0.4 now when i changed the sequence (i.e. made the RV q5 to q0 in the joint distribution it gave the same RV q0, SI as 0.007. Ideally, it should be 0.4 only right ??given that i have changed the sequence of the sample too. The sample data (inputs) is independent with co-rel coeff being less than 0.01.

[jonathf]

Chaospy is inherently positional. If you rearange the order of something, that is fine, you just need to be consistent about it.

The following two methods should be equivalent for getting s0 and s1:

dist1 = cp.Normal(0, 1)
dist2 = cp.Uniform(-1, 1)
q0, q1 = chaospy.variable(2)

# METHOD 1
joint = cp.J(dist1, dist2)
expansion = cp.create_expansion(4, joint)
samples = joint.sample(1000)
evals_ = [foo(x[0], x[1]) for x in samples.T]
approx = cp.fit_regression(expansion, samples, evals)
s0 = cp.Var(cp.E_cond(approx, q0, joint), joint) / cp.Var(approx, joint)
s1 = cp.Var(cp.E_cond(approx, q1, joint), joint) / cp.Var(approx, joint)

# METHOD 2
joint = cp.J(dist2, dist1)
expansion = cp.create_expansion(4, joint)
samples = joint.sample(1000)
evals_ = [foo(x[1], x[0]) for x in samples.T]
approx = cp.fit_regression(expansion, samples, evals)
s0 = cp.Var(cp.E_cond(approx, q1, joint), joint) / cp.Var(approx, joint)
s1 = cp.Var(cp.E_cond(approx, q0, joint), joint) / cp.Var(approx, joint)

In other words, there is a bit of book keeping needed to ensure that the reordering stays the same. But if you manage to not make any mistake, you should have order invariance as you expect, yes.

[rachelgomez161999]

In statistics, polynomial regression is a form of regression analysis in which the relationship between the independent variable x and the dependent variable y is modelled as an nth degree polynomial in x. Polynomial regression fits a nonlinear relationship between the value of x and the corresponding conditional mean of y, denoted E(y |x). Although polynomial regression fits a nonlinear model to the data, as a statistical estimation problem it is linear, in the sense that the regression function E(y | x) is linear in the unknown parameters that are estimated from the data. For this reason, polynomial regression is considered to be a special case of multiple linear regression.

The explanatory (independent) variables resulting from the polynomial expansion of the "baseline" variables are known as higher-degree terms. Such variables are also used in classification settings.[1]
The goal of regression analysis is to model the expected value of a dependent variable y in terms of the value of an independent variable (or vector of independent variables) x. In simple linear regression, the model

{\displaystyle y=\beta _{0}+\beta _{1}x+\varepsilon ,,}{\displaystyle y=\beta _{0}+\beta _{1}x+\varepsilon ,,}
is used, where ε is an unobserved random error with mean zero conditioned on a scalar variable x. In this model, for each unit increase in the value of x, the conditional expectation of y increases by β1 units.

In many settings, such a linear relationship may not hold. For example, if we are modeling the yield of a chemical synthesis in terms of the temperature at which the synthesis takes place, we may find that the yield improves by increasing amounts for each unit increase in temperature. In this case, we might propose a quadratic model of the form

{\displaystyle y=\beta _{0}+\beta _{1}x+\beta _{2}x^{2}+\varepsilon .,}{\displaystyle y=\beta _{0}+\beta _{1}x+\beta _{2}x^{2}+\varepsilon .,}
In this model, when the temperature is increased from x to x + 1 units, the expected yield changes by {\displaystyle \beta _{1}+\beta _{2}(2x+1).}{\displaystyle \beta _{1}+\beta _{2}(2x+1).} (This can be seen by replacing x in this equation with x+1 and subtracting the equation in x from the equation in x+1.) For infinitesimal changes in x, the effect on y is given by the total derivative with respect to x: {\displaystyle \beta _{1}+2\beta _{2}x.}{\displaystyle \beta _{1}+2\beta _{2}x.} The fact that the change in yield depends on x is what makes the relationship between x and y nonlinear even though the model is linear in the parameters to be estimated.

In general, we can model the expected value of y as an nth degree polynomial, yielding the general polynomial regression model

{\displaystyle y=\beta _{0}+\beta _{1}x+\beta _{2}x^{2}+\beta _{3}x^{3}+\cdots +\beta _{n}x^{n}+\varepsilon .,}{\displaystyle y=\beta _{0}+\beta _{1}x+\beta _{2}x^{2}+\beta _{3}x^{3}+\cdots +\beta _{n}x^{n}+\varepsilon .,}
Conveniently, these models are all linear from the point of view of estimation, since the regression function is linear in terms of the unknown parameters β0, β1, .... Therefore, for least squares analysis, the computational and inferential problems of polynomial regression can be completely addressed using the techniques of multiple regression. This is done by treating x, x2, ... as being distinct independent variables in a multiple regression model.
Regards,
Rachel Gomez

[Sesha02]

Hi Jonath

I used below paper as reference (section 4.4.3 in the paper), calculated the Sobol indices using PCE coefficients. UQLab uses similar methodology in calculating Sobol indices.
https://ethz.ch/content/dam/ethz/special-interest/baug/ibk/risk-safety-and-uncertainty-dam/publications/reports/RSUQ-2015-008.pdf

In chaospy exported PCE coefficients using .todict() function. And calculated first and total order Sobol indices.
As below image shows they kind of reasonable, both first or total order indices sum to one

Below image Sobol first and total order indices calculated using inbuilt chaospy functions chaospy.Sens_m & chaospy.Sens_t respectively. As image shows below, Sobol indices does not sum to one and total order indices looks way off, given there is no two factor interaction of input variables on the response. First and total order indices suppose to be same, given the model behaves linear relationship between input variables and output response.

I am not able to figure it out what causing this when we use chaospy.Sens_m & chaospy.Sens_t functions to calculate Sobol indices. These functions kind of take longer to run and results they are giving seems not accurate.

In my model, I have 11 input variables. Using chaospy 4.3.10 & numpoly 1.2.5

Any help on this greatly appreciated. Thank you so much in advance.

Thanks

Seshagiri

[jonathf]

It is hard for me to answer without looking at code. Do you have a complete minimal snippet that I can have a look at?

[Sesha02]

Hi Jonath

Thank you so much for the reply.

Here is the code I used to calculate the Sobol indices second plot in my previous post for your reference.
chaospy cp.Sens_m and cp.Sens_t functions kind of predicting Sobol indices wrong for additive linear model without any interactions in the model.

import os
import warnings
import chaospy as cp
import pandas as pd
import numpy as np
from sklearn import linear_model as lm

doefilepath=r'C:\Temp'
evalfilename=r'Input_Data.csv'
pathtofiles=os.path.join(doefilepath,evalfilename)
df_eval=pd.read_csv(pathtofiles,sep=",", encoding='cp1252')
df_eval=df_eval.dropna()

inputvariables=df_eval.columns[0:11]
yvarname=df_eval.columns[11]

V1=cp.Uniform(6,40)
V2=cp.Uniform(0,34)
V3=cp.Uniform(20,80)
V4=cp.Uniform(2,30)
V5=cp.Uniform(0,55)
V6=cp.Uniform(2,30)
V7=cp.Uniform(0,55)
V8=cp.Uniform(2,30)
V9=cp.Uniform(0,55)
V10=cp.Uniform(2,30)
V11=cp.Uniform(0,55)

Joint=cp.J(V1,V2,V3,V4,V5,V6,V7,V8,V9,V10,V11)

warnings.filterwarnings("ignore", category=RuntimeWarning)
kws = {"fit_intercept": False}
model=lm.LinearRegression(**kws)

order=2
expansion = cp.generate_expansion(order, Joint, normed=True, cross_truncation=1)
doesamples=df_eval[inputvariables].values.T
evaluations=df_eval[yvarname].values
evaluations=evaluations.reshape(len(evaluations))
approx_solver = cp.fit_regression(expansion, doesamples, evaluations, model=model)

S1_chaos=pd.DataFrame(data=cp.Sens_m(approx_solver, Joint), index=inputvariables,columns=['First Order'])
ST_chaos=pd.DataFrame(data=cp.Sens_t(approx_solver, Joint), index=inputvariables,columns=['Total Order'])

SAll_chaos=pd.concat([S1_chaos,ST_chaos],axis=1)
sobolax=SAll_chaos.plot.bar(title='Chaospy Sobol Indices',width=0.75)
sobolax.set_ylim(0,1)
for pt in sobolax.patches:
sobolax.annotate(str(np.round(pt.get_height(),2)), (pt.get_x()+0.05, pt.get_height()+0.02),rotation=90,fontsize=8)

When I use approximate solver from chaospy PCE and use SALib to calculate Sobol indices, I am getting right answers.
Please refer the plot at the lost that is the correct answer.

from SALib.sample import sobol as sobolsample
from SALib.analyze import sobol as sobolanalyze
from SALib.plotting.bar import plot as barplot

problem = {
'num_vars': len(inputvariables),
'names': list(inputvariables),
'bounds': [[6,40],
[0,34],
[20,80],
[2,30],
[0,55],
[2,30],
[0,55],
[2,30],
[0,55],
[2,30],
[0,55]],
'dists': ['unif']*11,
'outputs':yvarname
}

doesamples_salib = sobolsample.sample(problem, 1024).T
evaluations_pred=approx_solver(*doesamples_salib)

Si = sobolanalyze.analyze(problem, evaluations_pred, print_to_console=False)

total_Si, first_Si, second_Si = Si.to_df()
sobolsalib=pd.concat([first_Si['S1'],total_Si['ST']],axis=1)

sobolax=sobolsalib.plot.bar(title='Sobol Indices - SALib',width=0.75)
sobolax.set_ylim(0,1)
for pt in sobolax.patches:
sobolax.annotate(str(np.round(pt.get_height(),2)), (pt.get_x() * 1.005, pt.get_height()+0.02),rotation=90,fontsize=8)

Thanks

Seshagiri

[Sesha02]

Like in the below paper reference, I tried to calculate the Sobol indices using PCE coefficients (first picture in my first post), that also i am getting wrong answers. That may be something to do with Isoprobabilistic transform of input variable distributions

https://ethz.ch/content/dam/ethz/special-interest/baug/ibk/risk-safety-and-uncertainty-dam/publications/reports/RSUQ-2015-008.pdf

Calculation proposed in the paper similar to below pictures for your reference - UQLab does this way

Is there a way to calculate Sobol Indices from PCE coefficient in chaospy, like UQLab. This method seems much faster

Thanks

Seshagiri

[jonathf]

Could you post the input dataset Input_Data.csv.

[jonathf]

the method you are showing is better, but I never got around to setting it up.

So the method you are proposing could be done in chaospy manually by using approx, a_vec = cp.fit_regression(.., retall=True).
You can also get alpha_vec = chaospy.monomial(...).exponent where the args are the same as expansion creation.

[Sesha02]

Hi Jonathan

Thanks for looking into it. Here is the data you are looking for.

Input_Data.csv

Thanks

Seshagiri

[jonathf]

I just made a new release of chaospy version 4.3.11. It contains an PCE specific implementation of the Sobol indices that you can try out: chaospy.FirstOrderSobol, chaospy.SecondOrderSobol and chaospy.TotalOrderSobel. They can be used with expansion and Fourier coefficients alone:

approx, coefficients = chaospy.fit_regression(order, nodes, evals, retall=True)
first_order_indices = chaospy.FirstOrderSobol(expansion, coefficients)

My time is somewhat limited these days, so I don't have the oportunity to test the implementation really. So consider it experimental in this first iteration. If you try it out, please let me know how it goes. This approach is clearly supperior for PCE and should in time be made the default.

[Sesha02]

Hi Jonathan

Thank you so much for quick implantation of the Sobol Indices suggestion.

I have tested the code implemented in version 4.3.11. This the results I got

The actual answer should be (using SALib, the code I posted earlier, using approx_solver from chaospy)

Thanks

Seshagiri

[Sesha02]

Hi Jonathan

I looked at the code, problem i feel, when variance calculated using coefficients, we should exclude the coefficient value of the constant term (a0)

Thanks

Seshagiri

[jonathf]

Yeah, you are right. The formula is incorrect for the numerator D. Thanks for taking the time to figure out my mistake.

I've made a fix and a dev release. Try it out with:

pip install chaospy==4.3.12.dev0

[Sesha02]

Hi Jonathan

Thank you so much for the quick fix. Tested the chaospy version 4.3.12.dev0 .
Its working fine now. Sobol indices are calculated correctly.
Tested with one case right now, i will do few more tests and update you later.

Thanks

Seshagiri

[Sesha02]

Hi Jonathan

There something wrong with the SecondOrderSobol function, its not giving results correctly.
Code kind of adding first order coefficient's also while calculating the second order Sobol indices using pce coefficient's

Thanks

Seshagiri

[jonathf]

Yeah there was a bug, but I had to stare at the code a bit to spot the error. I'vev made a patch now in 4.3.13. Lets hope all the errors are ironed out now.

Thanks for suggesting the improved way to calculate Sobol indices, btw. It is an obvious improvement to how it was done before.