Hi all, in the following link: http://support.sas.com/documentation/cdl/en/imlug/59656/HTML/default/viewer.htm#nonlinearoptexpls_sect19.htm 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;
... View more
