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(¶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;
