BookmarkSubscribeRSS Feed
trevi83
Calcite | Level 5

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

2 REPLIES 2
trevi83
Calcite | Level 5

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?

trevi83
Calcite | Level 5

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-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

Multiple Linear Regression in SAS

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.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 2 replies
  • 1017 views
  • 0 likes
  • 1 in conversation