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.
Don't miss out on SAS Innovate - Register now for the FREE Livestream!
Can't make it to Vegas? No problem! Watch our general sessions LIVE or on-demand starting April 17th. Hear from SAS execs, best-selling author Adam Grant, Hot Ones host Sean Evans, top tech journalist Kara Swisher, AI expert Cassie Kozyrkov, and the mind-blowing dance crew iLuminate! Plus, get access to over 20 breakout sessions.
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.