Estimating a Three-Parameter Weibull Curve
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: CAPW3 */
/* TITLE: Estimating a Three-Parameter Weibull Curve */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Capability Analysis, Weibull Distribution, */
/* PROCS: CAPABILITY NLP MACRO */
/* DATA: */
/* */
/* REF: */
/* MISC: NLP is in SAS/OR */
/* */
/****************************************************************/
options ps=60 ls=80;
data wire;
label length='Length in cm';
input length @@;
cards;
3.779 4.649 4.795 4.307 2.983
3.863 4.023 4.079 4.735 4.551
4.009 4.835 4.164 4.721 4.168
5.582 3.491 4.752 3.939 4.824
5.266 2.142 3.903 4.734 3.373
3.645 3.726 3.219 3.726 3.737
3.047 4.622 5.449 4.580 2.922
5.023 4.437 3.432 4.149 2.840
4.062 4.437 3.432 4.149 2.840
5.257 4.062 5.257 4.571 3.565
3.896 4.287 3.285 4.174 4.574
;
/*------------------------------------*/
/* Weibull ( Three Parameter ) Case */
/*------------------------------------*/
/* Calculate Intial Guess for Parameters */
proc univariate data=wire noprint;
var length;
output out=stats n=ss min=minx range=range;
run;
proc transpose data=wire out=row;
var length;
run;
data intvals;
merge row stats;
theta1= minx - .02*range;
call symput('ss',left(ss));
call symput('theta1',left(theta1));
run;
/* Get Intial Guesses From 2-Parameter Weibull */
proc capability data=wire noprint;
var length;
hist / weibull(theta=&theta1 noprint)
outfit=parms1 noplot;
run;
data stats;
merge intvals parms1;
run;
/* INEST= Dataset for PROC NLP */
data parms(type=est);
set stats;
keep _type_ sig c theta;
_type_='parms'; sig = _scale_;
c = _shape1_;
theta = _locatn_; output;
_type_='lb'; sig = 0.;
c = 0.;
theta = .; output;
_type_='ub'; sig = .;
c = .;
theta = minx; output;
run;
/* Use PROC NLP to iteratively estimate parameters */
proc nlp data=stats inest=parms tech=nrridg fd phes outest=values;
array col{*} col1-col&ss;
min logf;
parms sig, c, theta;
sum1 = 0;
sumlogx = 0;
do i=1 to ss;
coli = col{i} - theta;
sumlogx = sumlogx + log(coli);
sum1 = sum1 + ((coli/sig)**c);
end;
logf = sum1 - ss*log(c) + ss*c*log(sig) - (c-1)*sumlogx;
run;
data values;
set values;
if _type_="PARMS" then output; else delete;
keep sig c theta;
call symput('sig',left(sig));
call symput('c',left(c));
call symput('theta',left(theta));
run;
title '3-Parameter Weibull Estimates';
options nodate;
proc print data=values noobs; run;
title1;
title 'Three Parameter Weibull Fit';
proc capability data=wire noprint;
var length;
hist / weibull(c=&c sigma=&sig theta=&theta
color=yellow)
cframe = gray
cfill = blue
cbarline = white
nolegend;
inset weibull / format = 6.3
cframe = white
ctext = white
cfill = blue;
run;
goptions reset=all;