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;