new; cls; /* Davidon-Fletcher-Powell 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};
print "b*=" dfp(&ll,b0)';


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

proc dfp(&f,beta);
   local step,maxiter,tol,i,g,g_1,H,sk,yk,beta_1,d,j,fv,f:proc;
   maxiter=1e+3; tol=1e-4;
   /* main */
   g=gradp(&f,beta)';
   H=eye(rows(beta));
   d=-H*g;	
   /* line search */
   step=seqa(1,1,10000)/10000;
   fv=zeros(10000,1);
   j=1;
   do while j<=10000;
      fv[j]=f(beta+step[j]*d);
      j=j+1;
   endo;
   step=step[minindc(fv)];
   step=step-1/10000+seqa(1,1,2000)/10000/1000;
   fv=zeros(2000,1);
   j=1;
   do while j<=2000;
      fv[j]=f(beta+step[j]*d);
      j=j+1;
   endo;
   step=step[minindc(fv)];
   step=step-1/10000/1000+seqa(1,1,2000)/10000/1000/1000;
   fv=zeros(2000,1);
   j=1;
   do while j<=2000;
      fv[j]=f(beta+step[j]*d);
      j=j+1;
   endo;
   step=step[minindc(fv)];
   /* update */
   g_1=g;
   beta_1=beta;
   beta=beta+step*d;
   i=1;
   do while i<=maxiter;
      if abs(g)<tol; 
         break; 
      endif;
      /* main */
      g=gradp(&f,beta)'; 
      yk=g-g_1;
      sk=beta-beta_1;
      if sk'yk>0 and yk'H*yk>0;
         H=H+(sk*sk')/(yk'sk)-(H*yk*yk'H)/(yk'H*yk);
      else;
         H=H;
      endif;
      d=-H*g;
      /* line search */
      step=seqa(1,1,10000)/10000;
      fv=zeros(10000,1);
      j=1;
      do while j<=10000;
         fv[j]=f(beta+step[j]*d);
         j=j+1;
      endo;
      step=step[minindc(fv)];
      step=step-1/10000+seqa(1,1,2000)/10000/1000;
      fv=zeros(2000,1);
      j=1;
      do while j<=2000;
         fv[j]=f(beta+step[j]*d);
         j=j+1;
      endo;
      step=step[minindc(fv)];
      step=step-1/10000/1000+seqa(1,1,2000)/10000/1000/1000;
      fv=zeros(2000,1);
      j=1;
      do while j<=2000;
         fv[j]=f(beta+step[j]*d);
         j=j+1;
      endo;
      step=step[minindc(fv)];
      /* update */
      g_1=g;
      beta_1=beta;
      beta=beta+step*d;
      i=i+1;
   endo;
   retp(beta);
endp;