?  Gradient Check code
?  Requires Procs  dologl  and  derivs  .

Proc GRADCHEC;

supres @smpl;

mmake b bnames;
mat novar = nrow(b);

? compute numeric derivatives
NUMDRV g;
?print g;

?  analytic derivatives of LogL
dologl loglold;
const ifhess,0;
derivs;

smpl smpall;
msd(noprint) dnames;
?print @sum;   ? analytic gradient

? print comparison of analytic and numeric
smpl 1,novar;
unmake g numeric;
unmake @sum analytic;
error = (analytic-numeric)/analytic;
print analytic numeric error;

endproc;

?--------------------------------------
? Proc for numeric derivatives
Proc NUMDRV g;
local b novar bold bi eps1 eps2 eps loglp loglm ivar;

mmake b bnames;
mat novar = nrow(b);

? compute numeric derivatives

copy b g;
copy b bold;
do ivar=1,novar;
  copy bold b;
  set bi = b(ivar);
  ?set eps1 = 1.e-4;
  set eps1 = 1.e-6;
  ?set eps2 = abs(bi*1.e-3);
  set eps2 = abs(bi*1.e-5);
  set eps = eps1 + pos(eps2-eps1);  ? eps = max(eps1,eps2)
  set b(ivar) = bi+eps;
  unmake b bnames;
  dologl loglp;
  set b(ivar) = bi-eps;
  unmake b bnames;
  dologl loglm;
  set g(ivar) = (loglp-loglm)/(2*eps);
  set b(ivar) = bi;
  unmake b bnames;
enddo;

endproc;