//
// UQ avec méthode du chaos polynomial (Hermite), toy problem
//  
////////////////////////////////////
///  4 direct computations
//
X1=-sqrt(3+sqrt(6))
X2=-sqrt(3-sqrt(6))
X3=sqrt(3-sqrt(6))
X4=sqrt(3+sqrt(6))
//
  a=X1+2
  N=10000
  S1=1
  for i=1:N
    S1=S1+1/N*(-a*S1) 
  end
///////////////////////////
  a=X2+2
  N=10000
  S2=1
  for i=1:N
    S2=S2+1/N*(-a*S2) 
  end
///////////////////////////
  a=X3+2
  N=10000
  S3=1
  for i=1:N
    S3=S3+1/N*(-a*S3) 
  end
///////////////////////////
  a=X4+2
  N=10000
  S4=1
  for i=1:N
    S4=S4+1/N*(-a*S4) 
  end
// 
////////////////////////////////////
deff('y=H0(x)','y=1')
deff('y=H1(x)','y=x')
deff('y=H2(x)','y=x^2-1')
deff('y=H3(x)','y=x^3-3*x')
deff('y=H4(x)','y=x^4-6*x^2+3')
//
X1=-sqrt(3+sqrt(6))
X2=-sqrt(3-sqrt(6))
X3=sqrt(3-sqrt(6))
X4=sqrt(3+sqrt(6))
// poids 
omega1=(4*3*2)/(4*H3(X1))^2;
omega2=(4*3*2)/(4*H3(X2))^2;
omega3=(4*3*2)/(4*H3(X3))^2;
omega4=(4*3*2)/(4*H3(X4))^2;
//
//
alpha0=omega1*S1*H0(X1)+omega2*S2*H0(X2)+omega3*S3*H0(X3)+omega4*S4*H0(X4);
alpha1=omega1*S1*H1(X1)+omega2*S2*H1(X2)+omega3*S3*H1(X3)+omega4*S4*H1(X4);
alpha2=omega1*S1*H2(X1)+omega2*S2*H2(X2)+omega3*S3*H2(X3)+omega4*S4*H2(X4);
alpha3=omega1*S1*H3(X1)+omega2*S2*H3(X2)+omega3*S3*H3(X3)+omega4*S4*H3(X4);
//
//
disp('mean of S with polynomial chaos');disp(alpha0);
disp('mean of S valeur exacte');disp(%e^(-(3/2)))
//
varappro=((alpha1/1)^2+(alpha2/2)^2+(alpha3/6)^2)
disp('variance of S with polynomial chaos');disp(varappro)
disp('variance of S valeur exacte');disp(%e^(-2)-%e^(-3))
//
disp('CV at entrance');disp(1/2)
disp('CV at exit');disp(sqrt(%e^(-2)-%e^(-3))/(%e^(-(3/2))))
//
//// simulation of the pdf of the polynomial approximation
//
Ntest=100000;
a_mean=2;
a_var=1;
Stot=[]
for j=1:Ntest
 a=a_mean+sqrt(a_var)*rand(1,1,'normal')  
 Srand=(alpha0/1)*H0(a-2)+alpha1/2*H1(a-2)+alpha2/2*H2(a-2)+alpha3/6*H3(a-2)  
 Stot=[Stot,Srand]       
end
figure(2)
clf()
histplot(100,Stot)

