stacksize(5000000)
x=3;//niveau x pour lequel on calcule P(G>x)
cdfnor("PQ",-x,0,1)//valeur de P(G>x)

N=1000;//nombre de runs independants
M=100;//nb de particules
xbasc();

//Monte Carlo standard
X=sum((rand(M,N,'g')>x),'r')/M;
histplot(10,X,style=2);
mean(X)
variance(X)
1.96*sqrt(variance(X))/mean(X)


//Avec fonction d'importance
Y=rand(M,N,'g');
Xi=(Y>0).*exp(-x*Y-x^2/2);
Xi=sum(Xi,'r')/M;
histplot(10,Xi,style=3)
mean(Xi)
variance(Xi)
1.96*sqrt(variance(Xi))/mean(Xi)

//par système de particules avec trajectoires et sans resampling à l'instant initial
n=5;//nombre de pas de temps
sp=[];
for i=1:N,
Y=rand(1:M,'g')/sqrt(n);
X=Y;
norm=1;
for p=1:(n-1);
 //selection suivant la valeur du dernier accroissement
  norm=norm*mean(exp(x*Y));
  w=cumsum(exp(x*Y))/sum(exp(x*Y));
  Z=zeros(1:M);
  k=1; 
  l=1;
  r=rand(1,'u');
  while (k<=M)
    while ((k-r)/M<w(l))
      Z(k)=X(l);
      k=k+1;
      r=rand(1,'u');//commenter cette ligne pour passer de stratifié à systématique
    end;
    l=l+1;
  end;
  X=Z;
  //mutation
  Y=rand(1:M,'g')/sqrt(n);
  X=X+Y;
  end;
sp=[sp,norm*sum((X>x).*exp(-x*(X-Y)))/M];
end;
mean(sp)
variance(sp)
1.96*sqrt(variance(sp))/mean(sp)


histplot(10,sp);

cdfnor("PQ",-x,0,1)