new; cls; /* Marquart method  */
x={4.2,4.8,6.0,3.7,4.5,3.2,5.3,4.9,4.8,6.8};
y={5.4,5.9,5.9,6.7,3.7,2.3,4.2,3.2,7.3,6.6};
x1=ones(rows(x),1)~x;
b0={1,1};
fn f(b)=x1*b;
print "b*=" marquart(&f,b0)';



proc marquart(&f,beta);
   local maxiter,tol,lambda,i,Zs,nb_1,nb,ys_1,ys,betas,f:proc;
   maxiter=1e+3; tol=1e-8; nb=1;
   lambda=sumc(sumc(gradp(&f,beta)));
   nb=sqrt(beta'beta);
   ys=y-f(beta);
   i=1;
   do while i<=maxiter;
      nb_1=nb;
      ys_1=ys;
      Zs=gradp(&f,beta);
      betas=inv(Zs'Zs+lambda*eye(rows(beta)))*Zs'ys;
      beta=beta+betas;
      nb=sqrt(beta'beta);
      ys=y-f(beta);
      if abs((nb-nb_1)/nb_1)<=tol;
         break;
      endif;
      if ys'ys<ys_1'ys_1;
         lambda=sqrt(0.5)*lambda;
      else;
         lambda=sqrt(2)*lambda;
      endif;
      i=i+1;
   endo;
   retp(beta);
endp;