function indexOutliers=BWTL(inputs,confidence)

%
% indexOutliers=BWTL(inputs,confidence)
%
% Function to perform Box and Whisker test and tolerance limit test:
%
% R. W. Walpole, R. H. Myers, A. Myers, and K. Ye, Probability & Statistics for Engineers and
% Scientists, 8th ed. Upper Saddle River, NJ: Prentice-Hall, 2007.
%
% Dongrui WU (dongruiw@usc.edu) 4/25/2009
%
% inputs: input vector
% confidence: confidence level used to determine k
%
% indexOutliers: indices of the outliers

lengInputs=length(inputs);
sortedInputs=sort(inputs);
index25=lengInputs*.25;  %% 1st Quantile
%% if index25 is an integer
if mod(index25,1)==0
    Q25=sortedInputs(index25);
else
    Q25=sortedInputs(floor(index25))*(ceil(index25)-index25)+...
        sortedInputs(ceil(index25))*(index25-floor(index25)); %% linear interpolation
end

index75=lengInputs*.75; %% 3rd quantile
if mod(index75,1)==0
    Q75=sortedInputs(index75);
else
    Q75=sortedInputs(floor(index75))*(ceil(index75)-index75)+...
        sortedInputs(ceil(index75))*(index75-floor(index75)); %% lnear interpolation
end

%% Box and Whisker Test
IQR=Q75-Q25; %% inter-quantile range
indexOutliers=[];
for i=1:lengInputs
    if inputs(i)<Q25-1.5*IQR | inputs(i)>Q75+1.5*IQR
        indexOutliers=[indexOutliers i];
    end
end

%% Tolerance-limit test
tempInputs=inputs;
tempInputs(indexOutliers)=[];
m=mean(tempInputs);
sigma=std(tempInputs);
k=findK(lengInputs,confidence); %% find the parameter k
lowerBound=m-k*sigma;
upperBound=m+k*sigma;

for i=1:length(tempInputs)
    if tempInputs(i)<lowerBound | tempInputs(i)>upperBound
        index=find(inputs==tempInputs(i));
        indexOutliers=[indexOutliers index];
    end
end

%% remove duplicate indices
indexOutliers=unique(indexOutliers);

