As the first step in the decommissioning of sasCommunity.org the site has been converted to read-only mode.

Here are some tips for How to share your SAS knowledge with your professional network.

Constants for Two-Sided Normal Tolerance Intervals

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;
put "&error invalid sample size " &n;
end;

if (&gamma <= 0) I (&gamma >= 1) then do;
put "&error invalid confidence level "
&gamma;
end;

if (&P <= 0) | (&P >= 1) then do;
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));
/* Compute actual confidence */

/* Main Newton's method loop */
kiters = 0;
do while ((kiters < &maxiters) &
(abs(actcover- &gamma) >&epsilon));
&k = &k - (actcover-&gamma)/actderiv;
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 */
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 */
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 */
integrnd = concover*exp(-&n*x**2/2);
dergrnd = dercover*exp(-&n*x**2/2);
return;

concover: /* Conditional coverage */
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;
do while ((riters < &maxiters) &
(abs(probr) >&epsilon));
r = r - probrlderiv;
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;```