function [k,options] = ToleranceFactor(n,coverage,confidence,options)
%ToleranceFactor computes the exact tolerance factor k for the two-sided
%(optionally also for the one-sided) p-content and gamma-confidence
%tolerance interval
% TI = [Xmean - k * S, Xmean + k * S],
%where Xmean = mean(X), S = std(X), X = [X_1,...,X_n] is a random sample
%from the distribution N(mu,sig2) with unknown mean mu and variance sig2.
%
%The value of the tolerance factor k is determined such that the tolerance
%intervals with the confidence gamma cover at least the fraction p
%('coverage') of the distribution N(mu,sigma^2), i.e.
% Prob[ Prob( Xmean - k * S < X < Xmean + k * S ) >= p ]= gamma,
%for X ~ N(mu,sig2) which is independent with Xmean and S. For more details
%see e.g. Krishnamoorthy and Mathew (2009).
%
%Syntax:
%k = ToleranceFactor(n,coverage,confidence)
%[k,options] = ToleranceFactor(n,coverage,confidence)
%or
%k = ToleranceFactor(n,coverage,confidence,options)
%
%Input:
% - n: samlpe size
% - coverage:
% Prob( Xmean - k * S < X < Xmean + k * S ) >= coverage,
% - confidence:
% Prob[ Prob( Xmean - k * S < X < Xmean + k * S ) >= p ] = confidence.
% - options: structure with further optional parameters:
% - options.DegreesOfFreedom:
% set the degrees of freedom (different from nu = n-1),
% for one sample, with n observations, set nu = n - 1,
% for m samples, each with n observations, set nu = m * (n - 1).
% - options.Onesided:
% logical flag for calculating the factor of (upper) one-sided
% tolerance interval; default value: onesided = false;
% - options.NormalizingConstant
% normalizing constant c; default value: c = 1/n. Alternatively, for
% computing the nonsimultaneous tolerance factors in linear regression
% use use e.g. c = x'*inv(X'*X)*x, here X is the regression design
% matrix and x is the vector of predictors.
%
%Output:
% - k, the calculated tolerance factor for tolerance interval
%
%Example:
%Calculate the tolerance factor k for the two-sided statistical p-content
%and gamma-coverage tolerance interval, with p = 0.80, and gamma = 0.95.
%
%Generate the random sample X of size n = 7 from a normal distribution
%N(mu,sig2) with mu = 5, sig2 = 0.5^2 and estimate the tolerance interval,
%based on the estimated sample mean and the sample standard deviation.
%
%Check the quality of estimatd TI: Generate a new random sample Z of size
%N = 1000 from N(mu,sig2) and calculate proportion of the generated
%observations that are covered by the estimated TI.
%
%n = 7;
%p = 0.80;
%gamma = 0.95;
%k = ToleranceFactor(n,p,gamma)
%
%mu = 5; sig = 0.5; X = mu + sig * randn(1,n);
%Xmean = mean(X);
%S = std(X);
%TI = [Xmean - k * S, Xmean + k * S]
%
%N = 1000;
%Z = mu + sig * randn(1,N);
%prop = sum(TI(1) < Z & Z < TI(2))/N
%
%Dependence:
%stats\norminv.m, stats\chi2inv.m, stats\nctinv.m
%
%See also:
%MATLAB\fzero.m
%
%References:
%
%Krishnamoorthy K., Mathew T.: Statistical Tolerance Regions: Theory,
%Applications, and Computation. Wiley, ISBN: 978-0-470-38026-0, 512 pages,
%May 2009.
%
%ISO 16269-6:2005: Statistical interpretation of data - Part 6:
%Determination of statistical tolerance intervals.
%
%Janiga I., Garaj I.: Two-sided tolerance limits of normal distributions
%with unknown means and unknown common variability. MEASUREMENT SCIENCE
%REVIEW, Volume 3, Section 1, 2003, 75-78.
%
%Cite this algorithm as:
%Viktor Witkovsky: ToleranceFactor - The Matlab algorithm for computing the
%exact tolerance factors of the two-sided non-simultaneous tolerance
%intervals for normal distribution. MATLAB Central File Exchange, May 2009.
%http://www.mathworks.com/matlabcentral/fileexchange/ .
%
%Viktor Witkovsky
%Institute of Mesaurement Science
%Slovak Academy of Sciences
%Dubravska cesta 9
%84104 Bratislava
%Slovak Republic
%E-mail: witkovsky@savba.sk
%http://www.um.sav.sk/en/department-03/viktor-witkovsky.html
%(c) Viktor Witkovsky, 2009, (witkovsky@savba.sk)
%Ver.: 20-Aug-2009 12:02:16
error(nargchk(1, 4, nargin));
if nargin < 2
coverage = 0.95;
end
if nargin < 3
confidence = 0.95;
end
if nargin < 4
options = struct('DegreesOfFreedom',n-1,...
'NormalizingConstant',1/n,'Onesided',false);
% Gauss-Hermite nodes and weights (order 499)
options.GHweights = [2.335254455909132e-020 8.727520180512132e-020 ...
3.192492986472047e-019 1.143079966323183e-018 4.006410689872739e-018 ...
1.374645207412527e-017 4.617487133659739e-017 1.518531475276486e-016 ...
4.889535045861178e-016 1.541555393524301e-015 4.759049063256765e-015 ...
1.438705492841146e-014 4.259268561632164e-014 1.234894400387112e-013 ...
3.506518269122293e-013 9.751989883410978e-013 2.656440833077607e-012 ...
7.087862727356068e-012 1.852493672185278e-011 4.742879166739566e-011 ...
1.189565769188374e-010 2.922891421620524e-010 7.036084851819923e-010 ...
1.659432282079332e-009 3.834542792182359e-009 8.681771070540109e-009 ...
1.926009618049048e-008 4.186761739458309e-008 8.918290399352374e-008 ...
1.861581441292868e-007 3.807962364485096e-007 7.633542883937283e-007 ...
1.499667160723453e-006 2.887425542736802e-006 5.448606677924250e-006 ...
1.007697238150531e-005 1.826650426945954e-005 3.245428369842251e-005 ...
5.651840553955225e-005 9.647588492209113e-005 1.614239536456450e-004 ...
2.647564252726149e-004 4.256601318096663e-004 6.708503311910412e-004 ...
1.036435786038859e-003 1.569714986526230e-003 2.330601619290612e-003 ...
3.392274252602869e-003 4.840564691851943e-003 6.771570420208750e-003 ...
9.287030186384078e-003 1.248716483984685e-002 1.646097132612199e-002 ...
2.127435734221991e-002 2.695697939802246e-002 3.348912916697763e-002 ...
4.079040867403259e-002 4.871214063099133e-002 5.703538422726187e-002 ...
6.547601341067191e-002 7.369756694384425e-002 8.133157079705838e-002 ...
8.800390531735061e-002 9.336473180595961e-002 9.711870844600029e-002 ...
9.905188660428430e-002]';
options.GHnodes = [6.551018453076543e+000 6.449557694700710e+000 ...
6.348166092579258e+000 6.246842375771083e+000 6.145585283625421e+000 ...
6.044393565536637e+000 5.943265980703555e+000 5.842201297893355e+000 ...
5.741198295210179e+000 5.640255759867880e+000 5.539372487966921e+000 ...
5.438547284275569e+000 5.337778962014790e+000 5.237066342647089e+000 ...
5.136408255669010e+000 5.035803538407177e+000 4.935251035817865e+000 ...
4.834749600289886e+000 4.734298091450734e+000 4.633895375975809e+000 ...
4.533540327400811e+000 4.433231825937007e+000 4.332968758289295e+000 ...
4.232750017477094e+000 4.132574502657917e+000 4.032441118953432e+000 ...
3.932348777278103e+000 3.832296394170188e+000 3.732282891625054e+000 ...
3.632307196930786e+000 3.532368242505885e+000 3.432464965739177e+000 ...
3.332596308831606e+000 3.232761218640027e+000 3.132958646522919e+000 ...
3.033187548187789e+000 2.933446883540434e+000 2.833735616535785e+000 ...
2.734052715030375e+000 2.634397150636366e+000 2.534767898577027e+000 ...
2.435163937543635e+000 2.335584249553732e+000 2.236027819810664e+000 ...
2.136493636564379e+000 2.036980690973368e+000 1.937487976967762e+000 ...
1.838014491113478e+000 1.738559232477368e+000 1.639121202493338e+000 ...
1.539699404829372e+000 1.440292845255389e+000 1.340900531511917e+000 ...
1.241521473179509e+000 1.142154681548851e+000 1.042799169491520e+000 ...
9.434539513313326e-001 8.441180427162521e-001 7.447904604907797e-001 ...
6.454702225688125e-001 5.461563478068849e-001 4.468478558777696e-001 ...
3.475437671443922e-001 2.482431025339954e-001 1.489448834125141e-001 ...
4.964813145912062e-002]';
end
% Set the optional parameters
c = options.NormalizingConstant;
onesided = options.Onesided;
nu = options.DegreesOfFreedom;
GHweights = options.GHweights;
GHnodes = options.GHnodes;
% Set values for limit cases
if confidence == 0
k = NaN;
return
elseif confidence == 1
k = Inf;
return
elseif coverage == 1
k = Inf;
return
elseif coverage == 0
k = 0;
return
end
if nu == Inf
k = norminv(1-(1-coverage)/2);
return
elseif nu <= 0
error('VW:ToleranceFactor','Degrees of freedom should be positive ...')
end
% Get the starting point: Howe's approximation
if nu < 2^20,
kH = sqrt( (nu * (1+c) * norminv((1-coverage)/2)^2) ...
/ chi2inv(1-confidence,nu));
else
kH = sqrt(((1+c) * norminv((1-coverage)/2)^2));
end
k0 = kH;
% Compute the tolerance factor
if onesided
k = sqrt(c)*nctinv(confidence,nu,norminv(coverage)/sqrt(c));
else
x = sqrt(2 * c) * GHnodes;
root = FindRoot(x,coverage);
ncx2points = nu * root.^2;
k = fzero(@(k) Integral(k,nu,GHweights,ncx2points)- confidence,k0);
end
%% Auxiliary function Integral
function int = Integral(k,nu,GHweights,ncx2points)
%Auxiliary function Integral for computing two sided tolerance intervals.
%Integral evaluates the integral defined by eq. (1.2.4) in Krishnamoorthy
%and Mathew: Statistical Tolerance Regions, 2009, p.7, by the Gauss-Hermite
%quadrature.
%(c) Viktor Witkovsky (witkovsky@savba.sk)
%Ver.: 12-Aug-2009 11:31:14
x = ncx2points / k^2;
fun = gammainc(x/2,nu/2);
int = 1.0 - 1.128379167095513 * GHweights' * fun;
%% Auxiliry function FindRoot
function r = FindRoot(x,coverage,tol,maxiter)
%ROOTFIDER numerically finds the root(s) r, the solution to the equation
%normcdf(x+r) - normcdf(x-r) = coverage,
%given x and coverage by the Halley's method
%(c) Viktor Witkovsky (witkovsky@savba.sk)
%Ver.: 20-Aug-2009 12:49:42
if nargin < 4
maxiter = 100;
end
if nargin < 3
tol = 3 * eps;
end
r = x;
iter = 0;
while true
iter = iter + 1;
[fun,funD1,funD2] = Coverage(x,r,coverage);
r = r - 2 * fun .* funD1 ./ (2 * funD1.^2 - fun .* funD2);
if iter > maxiter
break
end
if all(abs(fun) < tol)
break
end
end
%% Auxiliary function Coverage
function [fun,funD1,funD2] = Coverage(x,r,p)
%Auxiliary function Coverage: fun = normcdf(x+r) - normcdf(x-r);
%(c) Viktor Witkovsky (witkovsky@savba.sk)
%Ver.: 20-Aug-2009 12:49:42
if nargin < 2
p = 0;
end
sqrt2 = 1.414213562373095;
sqrt2pi = 2.506628274631;
fun = 0.5 * ( erfc(-(x+r)/sqrt2) - erfc(-(x-r)/sqrt2) ) - p;
aux1 = exp(-0.5 * (x + r).^2);
aux2 = exp(-0.5 * (x - r).^2);
funD1 = (aux1 + aux2)/sqrt2pi;
funD2 = ((x - r) .* aux2 - (x + r) .* aux1)/sqrt2pi;
%%