new; cls; /* Gaussian 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*=" gausm(&f,b0)';




proc gausm(&f,beta);
   local maxiter,tol,i,Zs,ys,beta_1,beta0,f:proc;
   maxiter=1e+3; tol=1e-8;
   i=1;
   do while i<=maxiter;
      beta_1=beta;
      Zs=gradp(&f,beta);
      ys=y-f(beta);
      beta0=inv(Zs'Zs)*Zs'ys;
      beta=beta+beta0;
      if abs((beta-beta_1)/beta_1)<tol;
         break;
      endif;
      i=i+1;
   endo;
   retp(beta);
endp;