MLEs for Two-Parameter Weibull Distribution

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: IMNLPEX4                                            */
/*   TITLE: MLEs for Two-Parameter Weibull Distribution         */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswmh/sasghg               UPDATE: frwick 5/2015   */
/*     REF: Technical Report Examples                           */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

proc iml;
carcin = { 143  164  188  188  190  192  206
           209  213  216  220  227  230  234
           246  265  304  216  244 };
p = ncol(carcin);
m = p - 2;
minx = carcin[><];
rang = carcin[<>] - minx;

/*--- 1. Two parameter estimation: for thet = 0 ---*/
start f_weib2(x) global(carcin,thet);
   /* use x[1]=sig and x[2]=c */
   p = ncol(carcin); m = p - 2;
   temp = carcin - thet;
   sum1 = sum(log(temp[1:m]));
   sum2 = sum((temp / x[1])##x[2]);
   f = m*log(x[2]) - m*x[2]*log(x[1]) + (x[2]-1)*sum1 - sum2;
   return(f);
finish f_weib2;

start g_weib2(x) global(carcin,thet);
   /* use x[1]=sig and x[2]=c */
   p = ncol(carcin); m = p - 2;
   g = j(1,2,0.);
   temp = carcin - thet;
   sum1 = sum(log(temp[1:m]));
   sum2 = sum((temp / x[1])##x[2]);
   sum3 = sum(((temp / x[1])##x[2]) # (log(temp / x[1])));
   g[1] = -m * x[2] / x[1] + sum2 * x[2] / x[1];
   g[2] = m / x[2] - m * log(x[1]) + sum1 - sum3;
   return(g);
finish g_weib2;

n = 2; thet = 0.;
x0 = j(1,n,.5);
optn = {1 2};
con = { 1.e-6 1.e-6 ,
     .     .    };
CALL NLPTR(rc,xopt,"f_weib2",x0,optn,con,,,,"g_weib2");
fopt = f_weib2(xopt);

/*--- Compute confidence interval for x at optimum ---*/
/* compute Hessian at optimum */
xopt = xopt`;
call nlpfdd(f,g,hes2,"f_weib2",xopt,,"g_weib2");
hin2 = inv(hes2);

/* quantile of normal distribution */
prob = 0.05;
stderr = sqrt(abs(vecdiag(hin2)));
z = quantile("Normal", 1 - prob/2);
xlb = xopt - z * stderr;
xub = xopt + z * stderr;
results = xlb || xopt || xub;
print results[L="Normal Distribution Confidence Interval"
              c={"Lower Bound" "Estimate" "UpperBound"}
              r={"sigma" "c"}];