/*
**   maxlik6.e
**
**   This command file illustrates maximum likelihood estimation
**   of a nonlinear least squares problem using MAXLIK.
**   The data for this example have been created using nlls.sim.
**
**   This example also illustrates the use of a user provided
**   numerical gradient.  This example uses a procedure that
**   comes with MAXLIK that computes a forward difference
**   numerical gradient using the Richardson extrapolation method.
**
**   Both the built-in forward difference method is compared with
**   an analytical gradient.
*/

library maxlik;
#include maxlik.ext;

@ -- procedure to compute log-likelihood -- @

proc lnlk(b,z);
    local dev,s2;
    dev = z[.,1] - (b[1] * exp(-b[2]*z[.,2]));
    s2 = dev'dev/rows(dev);
    retp(lnpdfn2(dev,s2));
endp;

proc grdlk(b,z);
    local d,s2,dev,r;
    d = exp(-b[2]*z[.,2]);
    dev = z[.,1] - b[1]*d;
    s2 = dev'dev/rows(dev);
    r = dev.*d/s2;
    retp(r~(-b[1]*z[.,2].*r));
endp;

proc hslk(b,z);
    local d,s2,dev,r, hss;
    d = exp(-b[2]*z[.,2]);
    dev = z[.,1] - b[1]*d;
    s2 = dev'dev/rows(dev);
    r = -z[.,2].*d.*(b[1].*d + dev);
    hss = -d.*d/s2~r~b[1].*z[.,2].*r;
    retp(xpnd(sumc(hss)));
endp;


startv = { 2, 1 };
output file = maxlik6.out reset;

maxset;
__title = "built-in forward difference numerical gradient";
__output = 5;
_max_Options = { forward };

{ x,f0,g,cov,retcode } = maxprt(maxlik("nlls",0,&lnlk,startv));

maxset;
__title = "built-in central difference numerical gradient";

_max_Options = { central };

{ x,f0,g,cov,retcode } = maxprt(maxlik("nlls",0,&lnlk,startv));

maxset;
__title = "Richardson extrapolation numerical gradient";

_max_UserGrad = &mlgradre;

{ x,f0,g,cov,retcode } = maxprt(maxlik("nlls",0,&lnlk,startv));


maxset;
__title = "analytical gradient";

_max_GradProc = &grdlk;

{ x,f0,g,cov,retcode } = maxprt(maxlik("nlls",0,&lnlk,startv));

maxset;
__title = "analytical Hessian and analytical gradient";

_max_GradProc = &grdlk;
_max_HessProc = &hslk;

{ x,f0,g,cov,retcode } = maxprt(maxlik("nlls",0,&lnlk,startv));

output off;