Calcite | Level 5

## Profile-Likelihood-Based Confidence Intervals

Hi all,

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);
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;
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;

2 REPLIES 2
Calcite | Level 5

## Profile-Likelihood-Based Confidence Intervals

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?

Calcite | Level 5

## Profile-Likelihood-Based Confidence Intervals

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.

From The DO Loop