Estimating a Three-Parameter Lognormal Curve
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: CAPL3 */
/* TITLE: Estimating a Three-Parameter Lognormal Curve */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Capability Analysis, Lognormal Distribution, */
/* PROCS: CAPABILITY NLP MACRO */
/* DATA: */
/* */
/* REF: */
/* MISC: NLP is in SAS/OR */
/* */
/****************************************************************/
options ps=60 ls=80;
data boxes;
label weight='Weight in grams';
input weight @@;
cards;
30.26 31.23 71.96 47.39 33.93 76.15 42.21
81.37 78.48 72.65 61.63 34.90 24.83 68.93
43.27 41.76 57.24 23.80 34.03 33.38 21.87
31.29 32.48 51.54 44.06 42.66 47.98 33.73
25.80 29.95 60.89 55.33 39.44 34.50 73.51
43.41 54.67 99.43 50.76 48.81 31.86 33.88
35.57 60.41 54.92 35.66 59.30 41.96 45.32
;
/*------------------------------------*/
/* Lognormal ( Three Parameter ) Case */
/*------------------------------------*/
/* Calculate Intial Guess for Parameters */
proc univariate data=boxes noprint;
var weight;
output out=stats n=ss min=minx range=range;
run;
proc transpose data=boxes out=row;
var weight;
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 Lognormal */
proc capability data=boxes noprint;
var weight;
hist / lnorm(theta=&theta1 noprint)
outfit=parms noplot;
run;
data stats;
merge intvals parms;
run;
/* INEST= Dataset For NLP */
data parms(type=est);
set stats;
keep _type_ sig zeta theta;
_type_='parms'; sig = _shape1_;
zeta = _scale_;
theta = _locatn_; output;
_type_='lb'; sig = 0.;
zeta = .;
theta = .; output;
_type_='ub'; sig = .;
zeta = .;
theta = minx; output;
run;
proc nlp data=stats inest=parms tech=nrridg fd phes outest=values;
array col{*} col1-col&ss;
min logf;
parms sig, zeta, theta;
sum1 = 0;
sum2 = 0;
pi = 3.141592654;
do i = 1 to ss;
coli = col{i} - theta;
sum1 = sum1 + log(coli);
sum2 = sum2 + (log(coli)-zeta)**2;
end;
logf = ss*log(sqrt(2*pi)*sig) + sum1 + (1/(2*sig**2))*sum2;
run;
data values;
set values;
if _type_="PARMS" then output; else delete;
keep sig zeta theta;
call symput('sig', left(sig));
call symput('zeta', left(zeta));
call symput('theta', left(theta));
run;
title '3-Parameter Lognormal Estimates';
options nodate;
proc print data=values noobs; run;
title 'Three Parameter Lognormal Fit';
proc capability data=boxes noprint;
var weight;
hist / lognormal(scale=&zeta sigma=&sig theta=&theta
color=yellow)
cframe = gray
cfill = blue
cbarline = white
nolegend;
inset lognormal / format = 6.3
cframe = white
ctext = white
cfill = blue
pos = ne;
run;
goptions reset=all;