new; cls; /* Steepest descent method with different step length */
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};
print "b*=" sd2(&ll,b0)';


proc ll(b);
   local n,e;
   n=rows(x1);
   e=y-x1*b;
   retp( -sumc( -1/2*(1+ln(2*pi)-ln(n))-1/2*e^2 ) );
endp;

proc sd2(&f,beta);
   local step,maxiter,tol,d,i,g,j1,j2,temp,zmin,ztemp,f:proc;
   maxiter=1e+3; tol=1e-5; d=1/10000;
   i=1;
   do while i<=maxiter;
      g=gradp(&f,beta);
      if abs(g)<tol; 
         break; 
      endif;
      step=d|d;
      zmin=f(beta-step.*g');
      j1=1;
      do while j1<=100;
         j2=1;
         do while j2<=100;
            temp=(j1|j2)*d;
            ztemp=f(beta-temp.*g');
            if ztemp<zmin;
               step=temp;
               zmin=ztemp;
            endif;
            j2=j2+1;
         endo;
         j1=j1+1;
      endo;
      beta=beta-step.*g';	
      i=i+1;
   endo;
   retp(beta);
endp;

/* original steepest descent */


proc sd(&f,beta);
   local step,maxiter,tol,i,g,j,fv,f:proc;
   maxiter=1e+3; tol=1e-5;
   i=1;
   do while i<=maxiter;
      g=gradp(&f,beta);
      if abs(g)<tol; 
         break; 
      endif;
      step=seqa(1,1,10000)/10000;
      fv=zeros(10000,1);
      j=1;
      do while j<=10000;
         fv[j]=f(beta-step[j]*g');
         j=j+1;
      endo;
      step=step[minindc(fv)];
      beta=beta-step*g';	
      i=i+1;
   endo;
   retp(beta);
endp;