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

function [res]=F(x,g)
   res= x - g;
endfunction

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

n=10000;

C=1;
beta=1;
X=estimation(n,C,beta);
plot2d([1:n],X);


C=1/10;
beta=1;
X=estimation(n,C,beta);
plot2d([1:n],X);

C=2;
beta=1;
X=estimation(n,C,beta);
plot2d([1:n],X);