//(x0,a,sigma,C,na)=(5,0,1,2,1/2),
//(x0,a,sigma,C,na)=(-3.795,0,1,10,0.000001),
//(x0,a,sigma,C,na)=(2.5,0,1,1,0.99),
//(x0,a,sigma,C,na)=(3,0,1,100,0.999),
//(x0,xi0,a,b,sigma,C,CI,na)=(4,2.8,0,2.8,1,170,0.008,0.999),

clear
x0=3;
xi0=2.8;
sigma=1;
C=100;
CI=0.008;
a=0;
b=2.8;
na=0.999;
T=100000;
X=zeros(1,T)+x0;

N=grand(1,T,'nor',a,sigma); 

for n=2:T

             if N(n-1)<=X(n-1) then U=1;
                               else U=0; 
             end

X(n)=X(n-1)+C*(na-U)/n;

end

xset("window",1);
xbasc();
xselect();
xx=[1:0.01:4];

yymax=1-((xx+xx^(-1))*(sqrt(2*%pi)))^(-1).*exp(-xx^2/2);
yymin=1-(xx*(sqrt(2*%pi)))^(-1).*exp(-xx^2/2);
plot2d(xx,yymin,1)
plot2d(xx,yymax,2)



xset("window",2);
xbasc();
xselect();
plot2d(1:T,X)







XI=zeros(1,T)+xi0;

NI=grand(1,T,'nor',b,sigma); 

for n=2:T

             if NI(n-1)<=X(n-1) then UI=exp(-NI(n-1)^2/2+(NI(n-1)-b)^2/2);
                               else UI=0; 
             end

XI(n)=XI(n-1)+CI*(na-UI)/n;

end
xset("window",3);
xbasc();
xselect();
plot2d(1:T,XI)