/////////////////////////////////////////////////////////////////////
//  
//  Simulation de la loi d'une v.a. X restreinte a l'évènement  X>a
//  Estimation de P(X>a) par des arbres genealogiques
//  Application/comparaisons pour X de loi N(0,1)
//
//////////////////////////////////////////////////////////////////////
//QQ Valeurs N=2000, T=300 Dalpha=0.01 a=0.9 estime N(0,1) restreinte a >=3=alpha(T+1)=T*Dalpha 
//(Rmq: alpha(1)=0, alpha(2)=Dalpha, alpha(3)=2*Dalpha,...)
//La configuration X(:,T+1) approche la loi restreinte a >=3=alpha(T+1)=T*Dalpha
// Autre jeu de valeurs : N=2000, T=500, Dalpha=.01, a=0.9 estime les niveaux >=5.
//////
//"Bonne" estimation de l'arbre pour N=10, T=300,a=.99, Dalpha=0.01 (erreur de 0.001)
//N=10, T=500,a=.99, Dalpha=0.01 (erreur de 2.10^(-7))
//N=50 ou 100 T=100 ou 300 a=.999, Dalpha=0.01 
//N=100, a=0.9, T=1000, Dalpha=0.01
//N=500 ou 1000, T=500 ou 1000, a=0.9, Dalpha=0.025, 0.02 ou 0.01 (precision 10^(-24))
//////////////////////////////////

clear
stacksize(90000000)
N=2000;
T=500;
X=rand(N,T+1,'n');
WW=rand(N,T+1,'n');
a=.9;
A=a*ones(N,1);
B=sqrt(1-a^2)*ones(N,1);
alpha=zeros(1,T);
Dalpha=.01;


for t=2:T+1
alpha(t)=alpha(t-1)+Dalpha;
end

w0=find(X(:,1)< alpha(1));
l0=length(w0);

v0=find(X(:,1)>= alpha(1));
q0=length(v0);

proba=(q0/N)*ones(1,T);


Y0=X(v0,1);
U0=floor(1+q0*grand(1,l0,'def'));

     if w0<>[] then
     X(w0,1)=Y0(U0,1);
     end
     
    


for t=1:T
X(:,t+1)=A.*X(:,t)+B.*WW(:,t+1);

W=find(X(:,t+1)<alpha(t));
L=length(W);
V=find(X(:,t+1)>= alpha(t));
Q=length(V);

X(W,t+1)=X(W,t);
    




w=find(X(:,t+1)< alpha(t+1));
l=length(w);

v=find(X(:,t+1)>= alpha(t+1));
q=length(v);



proba(t+1)=proba(t)*(q/N);

Y=X(v,1:t+1);
U=floor(1+q*grand(1,l,'def'));

     if w<>[] then
     X(w,1:t+1)=Y(U,1:t+1);
     end
     
     
     
//Evolution de l'arbre
     //xset("window",1);
     //xbasc();
     //xselect();
     //plot2d(1:t+1,X(:,1:t+1)');
     //xtitle("Evolution de l arbre genealogique : "+string(N),"","");
     

end





xset("window",2);
xbasc();
xselect();
histplot([Dalpha*T-1:.1:Dalpha*T+2],X(:,T+1)',style=2);
deff('[y]=g(x)','y=(sqrt(2/%pi)*exp(-x^2/2)/erfc(Dalpha*T/sqrt(2)));'); 
x=[Dalpha*T:.1:Dalpha*T+2];
plot2d(x,g(x),5)
xtitle("Histogramme gaussien et empirique restreints, Nombre de particules : "+string(N),"","")


 
   


xset("window",3);
xbasc();
xselect();
xx=[Dalpha*T-.5:Dalpha:Dalpha*T+.5];
yy=1/2*erfc(xx/sqrt(2));
yypart=proba(T+1)*ones(xx);
plot2d(xx,yy,5)
plot2d(xx,yypart,3)
xtitle("Estimation part. (en vert), et repartition gaussiennes (reelle en rouge) Nombre de particules : "+string(N),"","")



xset("window",4);
xbasc();
xselect();
yy=1/2*erfc(alpha(:)/sqrt(2));
plot2d(alpha(:),proba(:),5)
plot2d(alpha(:),yy(:),2)
xtitle("Proba. particulaire de depassement des niveaux (theorique en bleu), nombre de particules : "+string(N),"","")