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"}];