
%linear additive
clear all
close all
  disp('RandR -- Impedance2')
  set(0, 'DefaultAxesFontSize', 16);
  lw =2.5;  fs = 16;  msize = 10;

 model -- no interaction

impedance = [...
     37 38 37 41 41 40 41 42 41 42 41 43 42 42 42 43 42 43  ...
     30 31 31 31 31 31 29 30 28 42 43 42 43 43 43 42 42 42  ...
     28 30 29 29 30 29 31 29 29 42 42 43 45 45 45 44 46 45  ...
     25 26 27 28 28 30 29 27 27 40 40 40 43 42 42 43 43 41  ...
     25 25 25 27 29 28 26 26 26 35 34 34 35 35 34 35 34 35]' ;

% forming part and operator vectors.
a = repmat([1:10],9,1);  part = a(:);
b = repmat([1:3], 3,1);  operator = repmat(b(:),10,1);

[p table stats terms] =  anovan( impedance,{part, operator},...
'model','linear', 'varnames',{'part','operator'} )

MSE = table{4,5}
MSO = table{3,5}
MSP = table{2,5}
p = stats.nlevels(1); o = stats.nlevels(2); 
n=length(impedance)/(p * o);

s2Part = (MSP - MSE)/(o * n)
s2Oper = (MSO - MSE)/(p * n)
s2 = MSE

s2Repeat = s2

s2Repro = s2Oper

s2Gage = s2Repeat + s2Repro

s2Tot = s2Gage  + s2Part

sqrt(s2Gage/s2Tot)

%percent variation due to part
ps2Part = s2Part/s2Tot

%ps2Part = 0.9674

% signal-to-noise ratio. > 5 measuring system is capable
snr = sqrt( 2 * ps2Part/(1-ps2Part))
%snr = 7.6991
% rounded snr is NDC (number of distinct categories).
ndc = round(snr) %8

% percent of R&R variability
prr = sqrt(s2Gage/s2Tot)  


%MATLAB's built in function
gagerr(impedance,{part, operator},'model','linear')


