PDF of Estimator of Cpk
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: CPKDEN */
/* TITLE: PDF of Estimator of Cpk */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Capability Analysis, Capability Indices, */
/* PROCS: IML GPLOT */
/* DATA: */
/* */
/* REF: Chou, Y. and Owen, D. B. (1989). "On the */
/* Distributions of the Estimated Process Capability */
/* Indices". Communications in Statistics--Theory and */
/* Methods, 18. 4549-4560. */
/* */
/* Guirguis, G. H. and Rodriguez, R. N. (1992). */
/* "Computation of Owen's Q Function With Applications */
/* to Process Capability Analysis". Journal of Quality */
/* Technology 24. To appear. */
/* */
/* NOTES: This example provides a SAS/IML program to compute */
/* the probability density function of the usual */
/* estimator of the process capability index Cpk. The */
/* expression for the density function was derived by */
/* Chou and Owen (1989), and the computational method */
/* is given by Guirguis and Rodriguez (1992). */
/* */
/* The parameters of the density are: ds = d/sigma, */
/* ms = |mu-muo|/sigma and n = sample size. Here */
/* */
/* d = ( USL - LSL ) / 2 */
/* mu0 = ( USL + LSL ) / 2 */
/* */
/* LSL = lower specification limit */
/* USL = upper specification limit */
/* */
/* mu = process mean */
/* sigma = process standard deviation */
/* */
/* The program computes the density for the estimator */
/* in standardized form (scaled so that its expected */
/* value is zero and its standard deviation is one.) */
/* This enables one to compare the density with that */
/* of the standard normal distribution. */
/* */
/* */
/* MISC: Depending on the system that you use to run this *
/* sample, you may find that the run time for this */
/* program is long. */
/* */
/****************************************************************/
options ps=60 ls=80;
proc iml;
/* Assign the parameters and set up title for plot */
n = 30;
cpl = 0.1;
cpu = 0.1;
ran = 4.0;
m = 200;
title1 'Probability Density of Cpk Estimator for n=30';
/***************************************************************/
/* */
/* IML MODULES */
/* */
/***************************************************************/
/* Support of Non-Central-t distribution */
start supnct(df,nc,inf,sup);
ex = sqrt(df/2)*nc*exp(lgamma((df-1)/2)-lgamma(df/2));
va = ((df/(df-2))*(1+(nc*nc)))-(ex*ex);
sd = sqrt(va);
sup = ex+(4*sd);
inf = ex-(4*sd);
finish;
/* Non-central t distribution */
start cdfnct(x,df,nc,inf,sup);
if (x >= sup) then y = 1;
else if (x <= inf) then y = 0;
else y = probt(x,df,nc);
return(y);
finish;
/* Integrand for positive section of density Cpk estimator */
start fun(x) global(ny,nn,nu,nl);
y = exp(((nn-2)*log(x))-(0.5*x*x));
y = Y*(((x*x)-(nn-1))/(nn-1));
y = y*(probnorm((ny*x)-nl)-probnorm(-(ny*x)+nu));
return(y);
finish;
/* Density Cpk estimator */
start dencpk(y,n,cpu,cpl)
global(ny,nn,nu,nl,ng,r,infl,supl,infu,supu,mu,si);
na = 3*sqrt(n);
h = 0.000001;
ny = (na/sqrt(n-1))*y;
nn = n;
if (y <= 0) then do;
zu = (cdfnct((-na*y)+h,n-1,-nu,infu,supu)
-cdfnct(-na*y,n-1,-nu,infu,supu))/h;
zl = (cdfnct((na*y)+h,n-1,nl,infl,supl)
-cdfnct(na*y,n-1,nl,infl,supl))/h;
z = na*(zu+zl);
end;
else do;
int = 0.0||(r/y);
call quad(z,"fun",int) scale = si peak = 5.*mu;
z = (z*ng)/y;
end;
return(z);
finish;
/****************************************************************/
/* */
/* MAIN PROGRAM */
/* */
/****************************************************************/
r = (cpl+cpu)*sqrt(n-1)*(0.5);
ng = (n-1)/exp(lgamma((n-1)/2)+(log(2)*((n-3)/2)));
nu = 3*sqrt(n)*cpu;
nl = 3*sqrt(n)*cpl;
rig = ran;
lef = -ran;
inc = (rig-lef)/(m-1);
run supnct(n-1,-nu,infu,supu);
run supnct(n-1,nl,infl,supl);
do x = lef to rig by inc;
y = exp(-0.5*x*x)/sqrt(2*2*arsin(1));
if(x = lef) then u = 0||x||y;
else u = u//(0||x||y);
end;
ds = (3/2)*(cpl+cpu);
ms = (3/2)*(cpl-cpu);
f1 = (1/3)*sqrt((n-1)/2)*gamma((n-2)/2)*(1/gamma((n-1)/2));
f2 = sqrt(2/n)*(1/gamma(0.5))*exp(-n*0.5*ms*ms);
f3 = ms*(1-(2*probnorm(-sqrt(n)*ms)));
ex = f1*(ds-f2-f3);
va = ((n-1)/(9*(n-3)))*(ds**2-(2*ds*(f2+f3))+ms**2+(1/n));
va = va-(ex*ex);
sd = sqrt(va);
mu = ex;
si = sd;
do x = lef to rig by inc;
z = (sd*x)+ex;
y = dencpk(z,n,cpu,cpl)*sd;
u = u//(1||x||y);
end;
name = {aux x y};
create density from u[colname = name];
append from u;
close density;
data density;
set density;
length curve $ 15;
if (aux = 0) then curve = "N(0,1)";
else curve = "Cpk Estimator";
run;
* Plot the density function;
symbol1 i=join value=none l=1 w=2;
symbol2 i=join value=none l=3 w=2;
axis1 order = (-5.00 to 5.00 by 1.00) label = none;
axis2 order = ( 0.00 to 0.50 by 0.05) label = none;
proc gplot data = density;
plot y*x = curve /
haxis = axis1
vaxis = axis2
vminor = 0
hminor = 0
frame
;
run;
goptions reset=all;