Constants for Two-Sided Normal Tolerance Intervals

From sasCommunity
Jump to: navigation, search

Abstract

Two sided normal tolerance intervals are fundamental for the analysis of product reliability. These intervals are of the form x̄ ± ks. For sample size n the constant k is chosen to give the desired confidence that the interval contains a specified proportion of the true population. The constant k cannot be expressed in terms of standard statistical functions. This paper presents a SAS® macro that will compute values for the constants k for any values of the parameters. The computation uses an exact algorithm; in contrast, most published tables use approximations. SAS PROC CAPABILITY also uses an approximation.

Macro code

%macro twosidek(inset, outset, n=n, gamma= gamma, P = P, k = k);
    %let maxiters= 20; 
    %let epsilon = 0.0000001;
    %let error= %upcase(error);
    %let twopi = 2*3.1415926535;
    data &outset (keep= &n &gamma &P &k);
        set &inset;
        attrib
         &k label = "Tolerance Limit Constant"
         format= 8.3;
/* Three places only, because that is all */
/* that Odeh-Owen give. Formats for the */
/* other variables are inherited */
 
/* Error checking, for cases where */
/* constants aren't meaningful */
if &n < 2 then do;
    baddata = 1;
    put "&error invalid sample size " &n;
end;
 
if (&gamma <= 0) I (&gamma >= 1) then do;
    baddata = 1;
    put "&error invalid confidence level "
      &gamma;
end;
 
if (&P <= 0) | (&P >= 1) then do;
    baddata = 1;
    put "&error invalid reliability " &P;
end;
if baddata = 1 then return;
 /* To top of data step */
 
/* Initial guess for the main loop -- */
/* uses the formula from PROC CAPABILITY */
/* This is easier to compute than the */
/* Odeh-Owen initial guess */
&k = probit((1 + &P)I2)(1 + 11(2*&n))*
     sqrt((&n-1)/cinv(1-&gamma, &n-1));
link actcover;
 /* Compute actual confidence */
 
/* Main Newton's method loop */
kiters = 0;
do while ((kiters < &maxiters) &
  (abs(actcover- &gamma) >&epsilon));
 &k = &k - (actcover-&gamma)/actderiv;
 link actcover;
 kiters = kiters+1;
end;
 
/* If we exited the loop on the iteration */
/* count, then the computed value may not */
/* be so hot. So set it to missing. */
if kiters >= &maxiters then do;
 put "&error Maximum iteration count of "
 "&maxiters exceeded in main loop.";
 &k = .;
 end;
return;  /* To top of data step */
/* This ends the main loop -- all else */
/* is subprograms to be linked. */
 
actcover: /* Determines actual confidence */
/* and its derivative */
/* Integration by Simpson's rule */
%let ystep = 0.05;
 /* The step value from Odeh and Owen */
%let ylimit = 5;
 /* The tail integral is < 10**-6 */
xstep = &ystep/sqrt(&n);
intsteps = &ylimit/&ystep;
 
/* Initialize sums for Simpson's rule */
x = 0; link integrnd;
fof0 = integrnd; derof0 = dergrnd;
sum1 = 0; sum2 = 0;
dersum1 = 0; dersum2 = 0;
 
/* Simpson's rule loop */
do intstep = 0 to intsteps;
 x = 2*xstep*intstep;
 x = x + xstep;
 link integrnd;  /* Middle of interval */
  sum1 = sum1 + integrnd;
  dersum1 = dersum1 + dergrnd;
 x = x + xstep;  /* End of interval */
 link integrnd;
  sum2 = sum2 + integrnd;
  dersum2 = dersum2 + dergrnd;
 end;
actcover = (2*sqrt(&n)/sqrt(&twopi))*
 (xstepl3)*(fof0+4*sum1+2*sum2-integrnd);
/* The - dign is correct, because */
/* this value was included twice in */
/* the definition of sum2 */
actderiv = (2*sqrt(&n)/sqrt(&twopi))*
 (xstepl3)*(derof0 + 4*dersum1 +
  2*dersum2- dergrnd);
return;
 
integrnd:  /* Evaluate the integrands */
 link concover;
 integrnd = concover*exp(-&n*x**2/2);
 dergrnd = dercover*exp(-&n*x**2/2);
 return;
 
concover: /* Conditional coverage */
 link r;
 f = &n-1;
 u = f*(rl&k)**2;
 concover = 1- probchi(u, f);
 dercover = exp(-(fl2)*log(2) - lgamma(fl2)
  - ul2 + (fl2-1)*log(u) + log(2)
  + log(u)- log(&k));
 return;
 
r: /* Use Newton's method to solve for r */
   /* Initial value is the final value   */
   /* from the last time around. Special */
   /* case is x = 0 at start of integral */
 riters = 0;
 if x = 0 then r = probit((l+&P)/2);
 else do;
  link probr;
  do while ((riters < &maxiters) &
    (abs(probr) >&epsilon));
   r = r - probrlderiv;
   link probr;
  riters = riters + 1;
  end; end;
 if riters >= &maxiters then
  put "&error Maximum iteration count of "
     "&maxiters exceeded in main loop.";
return;
 
probr:
 probr = probnorm(x+r) - probnorm(x-r) - &P;
 deriv = (1lsqrt(&twopi)) *
  (exp(-(x+r)**212) + exp(-(x-r)••212));
 return;
 
run; /* Computational data step */
%mend;

Test Driver

/* Driver Program */
Options linesize = 72; 
 
Title "Two Sided Normal Tolerance Limit
  Factors (Control Center)";
 
data parmsin;
 attrib
  n label = "Sample Size"
  P label = "Reliability"
  gamma label = "Confidence"
;
/* Compute values to check against tables */
/* Also desired special values            */
do gamma= 0.500, 0.600, 0.750, 0.800,
  0.900, 0.950, 0.990;
 do P = 0.75, 0.80, 0.90, 0.95, 0.99;
  do n = 2 to 25, 30 to 50 by 5,
   60 to 100 by 10, 200, 300, 500, 1000;
  output;
  end; end; end;
run;
 
%twosidek(parmsin, tolcomp);
 
data tolcomp (keep= n P gamma Plbl k);
 set tolcomp;
 attrib Plbl format = 5.3; Plbl = P;
 attrib gamma format = 5.3;
 run;
proc sort data = tolcomp; by gamma n; run;
proc transpose data = tolcomp out = tables;
 by gamma n; id P; idlabel Plbl; var k; run;
data tables; set tables;
 drop _name __ label_;
 run;
proc print noobs label data = tables;
 by gamma; pageby gamma; run;

Online resources

Download the .pdf of the paper from here.