/*
**   maxlik13.e
**
**  This example illustrates the profile likelihood confidence
**  intervals.  This method does not rely on a quadratic approximation
**  to the log-likelihood surface as do the usual t-statistics which
**  may be refered to as the "Wald-type" confidence intervals.  As 
**  emphasized by many statisticians such as Bates and Watts, the
**  Wald-type confidence intervals may be unreliable for nonlinear
**  models.  This example calls MAXPflClimits which computes confidence
**  intervals using the log-likelihood surface itself rather than a
**  quadratic approximation.  One outcome of this method is that the
**  intervals may not be symmetric about the estimate.
*/

library maxlik;
#include maxlik.ext;
maxset;

@ -- 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 = maxlik13.out reset;

_max_GradProc = &grdlk;
_max_HessProc = &hslk;


_max_Alpha = .05;

__title = "Wald-type confidence limits, alpha = .05";
{ x,f0,g,cov,retcode } = maxlik("nlls",0,&lnlk,startv);
cl1 = maxtlimits(x,cov);
call maxclprt(x,f0,g,cl1,retcode);

__title = "Profile likelihood confidence limits, alpha = .05";
cl2 = maxpflclimits(x,f0,"nlls",0,&lnlk);
call maxclprt(x,f0,g,cl2,retcode);


_max_Alpha = .01;

__title = "Wald-type confidence limits, alpha = .01";
{ x,f0,g,cov,retcode } = maxlik("nlls",0,&lnlk,startv);
cl1 = maxtlimits(x,cov);
call maxclprt(x,f0,g,cl1,retcode);

__title = "Profile likelihood confidence limits, alpha = .01";
cl2 = maxpflclimits(x,f0,"nlls",0,&lnlk);
call maxclprt(x,f0,g,cl2,retcode);


_max_Alpha = .001;

__title = "Wald-type confidence limits, alpha = .001";
{ x,f0,g,cov,retcode } = maxlik("nlls",0,&lnlk,startv);
cl1 = maxtlimits(x,cov);
call maxclprt(x,f0,g,cl1,retcode);

__title = "Profile likelihood confidence limits, alpha = .001";
cl2 = maxpflclimits(x,f0,"nlls",0,&lnlk);
call maxclprt(x,f0,g,cl2,retcode);


output off;