# Constants for Two-Sided Normal Tolerance Intervals

From sasCommunity

Revision as of 08:42, 8 June 2016 by Paulkaefer (Talk | contribs)

## 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 " γ 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.