//
// UQ avec méthode Monte Carlo, toy problem
//
Ntest=100000;
amin=1;
amax=3;
Stot=[]
for k=1:Ntest
  a=amin+(amax-amin)*rand()
  N=100
  y=1
  for i=1:N
    y=y+1/N*(-a*y)  // Euler method for the computation of S
  end
  Stot=[Stot,y];
//disp('exact value');disp(exp(-a))
//disp('obtained value');disp(S)
end
Smin=exp(-3)
Smax=exp(-1)
Stest=linspace(Smin,Smax,100)
pdfS=ones(Stest)./(2*Stest)
figure(1)
clf()
histplot(Stest,Stot)
plot2d(Stest,(Smax-Smin)/100*pdfS)
//
disp('exact mean of S');disp(1/2*(exp(-1)-exp(-3)));
disp('mean of S (with MC)');disp(mean(Stot))
//
disp('exact variance of S');disp(1/4*(%e^(-2)-%e^(-6))-(1/2*(exp(-1)-exp(-3)))^2)
disp('variance of S (with MC)');disp((stdev(Stot))^2)
