Hi all,
in the following link:
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(¶m1/&scale_mu);
%let sigma_init=%sysevalf(¶m2/&scale_sigma);
%let gamma_init=%sysevalf(¶m3/&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;
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?
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.
SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!
Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.
Find more tutorials on the SAS Users YouTube channel.