function [res]=gamma(n,C,Beta)
  res= C / (n^Beta);
endfunction

function [res]=F(x,g,proba)
   if (abs(g) < x) then res=1-proba; 
      else res= - proba; end;
endfunction

function X=estimation(n,C,Beta)
  proba=0.95;
  X=zeros(1,n);
  G=grand(1,n,"nor",0,1);
  X(1) = 2;
  for k=[1:n-1] do
	X(k+1)=X(k)-gamma(k,C,Beta)*F(X(k),G(k),proba);
  end;
endfunction

// La valeur pour avoir le TCL 2c^* alpha > 1
x_star=1.96;
c_star=2*exp(-x_star^2/2)/sqrt(2*%pi);
alpha_star=1/(2*c_star);

n=1000;

C=2.0*alpha_star;
Beta=1.5;
X=estimation(n,C,Beta);
plot2d(1:n,X);

N=1000;
z=0;
for i=[1:N] do
  X=estimation(n,C,Beta);
  z(i)=X(n);
end
//z=z-mean(z);
//z=z/std(z);
histplot(50,z);