Profile-Likelihood-Based Confidence Intervals

trevi83 · Posted 03-28-2012 09:57 AM

Hi all,

in the following link:

http://support.sas.com/documentation/cdl/en/imlug/59656/HTML/default/viewer.htm#nonlinearoptexpls_se...

I find an example for conficence intervals but it's based on 2 equations.

Instead I have 3 equations to calculate....i have create a code for this case but the results are bad.

I think that the problem is "delt" and its calculation....

I hope someone can help me.

Thanks!

%macro conf_interval_lnorm_3_par;

%let mu_init=%sysevalf(&param1/&scale_mu);
%let sigma_init=%sysevalf(&param2/&scale_sigma);
%let gamma_init=%sysevalf(&param3/&scale_gamma);
%let bound_gamma=%sysevalf(&soglia/&scale_gamma);

proc iml;
use importi_netti_&distribuzione.;
read all var {x} into ing;
ingresso=t(ing);/* ingresso=row vector*/

xopt={&mu_init &sigma_init &gamma_init}; /* inizializes vectors */
xub={0,0,0};
xlb={0,0,0};
scala={&scale_mu 0 0 ,0 &scale_sigma 0,0 0 &scale_gamma};

start funzione(x) global(ingresso);/*function loglikelihood*/
    p=ncol(ingresso);
    if x[2] > 0 & x[3]< &bound_gamma & CDF('LOGNORMAL',&soglia-(x[3]*&scale_gamma),x[1]*&scale_mu,x[2]*&scale_sigma) ^= 1 then do;
     w=0;
     do i=1 to p;
      w=w-log((ingresso-(x[3]*&scale_gamma))*x[2]*&scale_sigma*sqrt(2*CONSTANT('PI'))) -
                  0.5*((log(ingresso-(x[3]*&scale_gamma))-(x[1]*&scale_mu))/(x[2]*&scale_sigma))**2;
     end;
     loglik=w-p*log(1-CDF('LOGNORMAL',&soglia-(x[3]*&scale_gamma),x[1]*&scale_mu,x[2]*&scale_sigma))-0.3*((x[2]*&scale_sigma)**2);
    end;
    else
     loglik=.;
    f=loglik;
    return (f);
finish funzione;

start derivate(x) global(ingresso); /* partial derivates loglikelihood */
   g = j(1,3,0.);
       p=ncol(ingresso);
       w=0;v=0;z=0;
       do i=1 to p;
      w=w+(log(ingresso-x[3]*&scale_gamma)-x[1]*&scale_mu)/(x[2]*&scale_sigma)**2;
      v=v+((log(ingresso-x[3]*&scale_gamma)-x[1]*&scale_mu)**2)/(x[2]*&scale_sigma)**3-1/(x[2]*&scale_sigma);
      z=z+(log(ingresso-x[3]*&scale_gamma)-x[1]*&scale_mu)/(x[2]*&scale_sigma)/((x[2]*&scale_sigma)*(ingresso-x[3]*&scale_gamma))-1/(ingresso-x[3]*&scale_gamma);
       end;
         T=(log(&soglia-x[3]*&scale_gamma)-x[1]*&scale_mu)/(x[2]*&scale_sigma);
      FLN=cdf('lognormal',&soglia-x[3]*&scale_gamma,x[1]*&scale_mu,x[2]*&scale_sigma);
      g[1]=w-p*(1/(x[2]*&scale_sigma))*pdf('normal',T)/(1-FLN);
       g[2]=v+p*pdf('normal',T)/(1-FLN)*(-T/(x[2]*&scale_sigma))-0.6*(x[2]*&scale_sigma);
       g[3]=z-p*pdf('normal',T)/(1-FLN)/((x[2]*&scale_sigma)*(&soglia-x[3]*&scale_gamma));
       return(g);
finish derivate;

start f_plln3(x) global(ingresso,ipar,lstar);
         like = funzione(x);
         grad = derivate(x);
         grad[ipar] = like - lstar;
         return(grad`);
finish f_plln3;

h={&h11 &h12 &h13, &h12 &h22 &h23, &h13 &h23 &h33};

/* quantile of chi**2 distribution */
fopt=-&loglikelihood;
prob=0.32;
chqua = cinv(1-prob,1);
lstar = fopt - .5 * chqua;
optn = {3 0};

con={. 1.e-6 1.e-6 . .,
. . &bound_gamma . .};

tc=j(1,12,.);
tc[1]=100000;
tc[2]=100000;

ddd=j(2,2,.);
v_col=j(2,1,.);

do ipar = 1 to 3;
      if ipar=1 then do;
         ind = 2;
         jnd = 3;
        end;
    else if ipar=2 then do;
         ind = 1;
         jnd = 3;
        end;
    else if ipar=3 then do;
         ind = 1;
         jnd = 2;
        end;
    print ipar ind jnd;
    a=h[ind,ind];
    b=h[jnd,ind];
    c=h[jnd,jnd];
    d=h[ipar,ind];
    e=h[ipar,jnd];
    f=h[ipar,ipar];
    ddd[1,1]=a;
    ddd[1,2]=b;
    ddd[2,1]=b;
    ddd[2,2]=c;
    v_col[1,1]=d;
    v_col[2,1]=e;
         delt = - inv(ddd) * v_col;
         alfa = -( f + delt` * v_col);

        if alfa > 0 then alfa = .5 * sqrt(chqua / alfa);
         else do;
            print "Bad alpha";
            alfa = .1 * xopt[ipar];
      print alfa;
         end;
         if ipar=1 then delt = 1 // delt;
         else if ipar=2 then do;
            delt = 1//delt;
         delt[1]=delt[2];
         delt[2]=1;
        end;
      else if ipar=3 then delt = delt // 1;

/* upper bound */
     x0 = (xopt + alfa * delt);
     con2 = con; con2[1,ipar] = xopt[ipar];
        CALL NLPLM(rc,xresu,"f_plln3",x0,optn,con2,tc);
        fu = f_plln3(xresu);
     su = ssq(fu);
     if xresu[ipar]=xopt[ipar] | su >= 1 then
         xub[ipar]=.;
         else
         xub[ipar]= xresu[ipar];

/* lower bound */
        x0 = (xopt - alfa * delt);
        con2= con; con2[2,ipar] = xopt[ipar];
        CALL NLPLM(rc,xresl,"f_plln3",x0,optn,con2,tc);
      fl = f_plln3(xresl);
     sl = ssq(fl);
     if xresl[ipar]=xopt[ipar] | sl >= 1 then
         xlb[ipar]=.;
         else
       xlb[ipar]= xresl[ipar];

   xopt[ipar]=xopt[ipar]*scala[ipar,ipar];
   xub[ipar]=xub[ipar]*scala[ipar,ipar];
   xlb[ipar]=xlb[ipar]*scala[ipar,ipar];

end;

    /* Output */
if xlb[1]^=. then do;
    CALL SYMPUT("CL_inf_param1",char(xlb[1]));
end;
if xlb[2]^=. then do;
    CALL SYMPUT("CL_inf_param2",char(xlb[2]));
end;
if xlb[3]^=. then do;
    CALL SYMPUT("CL_inf_param3",char(xlb[3]));
end;
if xub[1]^=. then do;
    CALL SYMPUT("CL_sup_param1",char(xub[1]));
end;
if xub[2]^=. then do;
    CALL SYMPUT("CL_sup_param2",char(xub[2]));
end;
if xub[3]^=. then do;
    CALL SYMPUT("CL_sup_param3",char(xub[3]));
end;
s=t(xopt);
print "Profile-Likelihood Confidence Interval";
      print xlb s xub;
quit;

%mend;

trevi83 · Posted 03-29-2012 08:55 AM

I submit this code for many input parameters (mu, sigma and gamma) and my problem is that the confidence

interval (lower or upper, depends of cases) not calculate.

What are the causes of non-computing of confidence intervals?

trevi83 · Posted 03-29-2012 10:10 AM

When I calculate the upper (or lower) bound, I add this condition "if xresu[ipar]=xopt[ipar] then xub[ipar]=.;" because the result (xresu) of CALL NLPLM is equal to optimization value (xopt) but I expect different value.