Chapter Contents 
Previous 
Next 
Nonlinear Optimization Examples 
This example calculates confidence intervals based on the profile likelihood for the parameters estimated in the previous example. The following introduction on profilelikelihood methods is based on the paper of Venzon and Moolgavkar (1988).
Let be the maximum likelihood estimate (MLE) of a parameter vector and let be the loglikelihood function defined for parameter values .
The profilelikelihood method reduces to a function of a single parameter of interest, , where , by treating the others as nuisance parameters and maximizing over them. The profile likelihood for is defined as
The following code defines the module for the system of n=2 nonlinear equations to be solved:
start f_plwei2(x) global(carcin,ipar,lstar); /* x[1]=sigma, x[2]=c */ like = f_weib2(x); grad = g_weib2(x); grad[ipar] = like  lstar; return(grad`); finish f_plwei2;
The following code implements the LevenbergMarquardt algorithm with the NLPLM subroutine to solve the system of two equations for the left and right endpoints of the interval. The starting point is the optimizer , as computed in the previous example, moved toward the left or right endpoint of the interval by an initial step (refer to Venzon and Moolgavkar 1988). This forces the algorithm to approach the specified endpoint.
/* quantile of chi**2 distribution */ chqua = cinv(1prob,1); lstar = fopt  .5 * chqua; optn = {2 0}; do ipar = 1 to 2; /* Compute initial step: */ /* Choose (alfa,delt) to go in right direction */ /* Venzon & Moolgavkar (1988), p.89 */ if ipar=1 then ind = 2; else ind = 1; delt =  inv(hes2[ind,ind]) * hes2[ind,ipar]; alfa =  (hes2[ipar,ipar]  delt` * hes2[ind,ipar]); if alfa > 0 then alfa = .5 * sqrt(chqua / alfa); else do; print "Bad alpha"; alfa = .1 * xopt[ipar]; end; if ipar=1 then delt = 1  delt; else delt = delt  1; /* Get upper end of interval */ x0 = xopt + (alfa * delt)`; /* set lower bound to optimal value */ con2 = con; con2[1,ipar] = xopt[ipar]; call nlplm(rc,xres,"f_plwei2",x0,optn,con2); f = f_plwei2(xres); s = ssq(f); if (s < 1.e6) then xub[ipar] = xres[ipar]; else xub[ipar] = .; /* Get lower end of interval */ x0 = xopt  (alfa * delt)`; /* reset lower bound and set upper bound to optimal value */ con2[1,ipar] = con[1,ipar]; con2[2,ipar] = xopt[ipar]; call nlplm(rc,xres,"f_plwei2",x0,optn,con2); f = f_plwei2(xres); s = ssq(f); if (s < 1.e6) then xlb[ipar] = xres[ipar]; else xlb[ipar] = .; end; print "ProfileLikelihood Confidence Interval"; print xlb xopt xub;
The results, shown in Output 11.5.1, are
close to the results shown in Output 11.4.2.

Chapter Contents 
Previous 
Next 
Top 
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.