Software to support ISO/TS 28037:2010(E)

Contents

Introduction

This document runs the numerical example of weighted least squares (WLS) with known unequal weights described in Clause 6 (Model for uncertainties associated with the y_i), and performs the prediction described in 11.1 EXAMPLE 2.

Users of MATLAB may run the code in the corresponding M-file directly to obtain the results given in ISO/TS 28037:2010(E) and may also modify the data and run the code on the new data.

For users who do not have access to MATLAB, the software may be used as the basis for preparing implementations in other programming languages.

The software is provided with a software licence agreement (REF: DS/L/24/001) and the use of the software is subject to the terms laid out in that agreement. By running the M-code, the user accepts the terms of the agreement.

close all
clear
format short

Assign measurement data

Assign x-values.

x = [1.0 2.0 3.0 4.0 5.0 6.0]';
m = length(x);

Assign y-values.

y = [3.2 4.3 7.6 8.6 11.7 12.8]';

Assign uncertainties associated with y-values.

uy = [0.5 0.5 0.5 1.0 1.0 1.0]';

Obtain estimates of the straight line calibration function parameters and associated standard uncertainties and covariance

Solve the weighted least squares problem to obtain best fit straight-line parameters.

Step 1.

w = ones(m, 1)./uy;
F2 = sum(w.*w);

Step 2.

g0 = (sum(w.*w.*x))/F2;
h0 = (sum(w.*w.*y))/F2;

Step 3.

g = w.*(x - g0);
h = w.*(y - h0);

Step 4.

G2 = sum(g.*g);

Step 5.

b = (sum(g.*h))/G2;
a = h0 - b*g0;

Step 6.

u2a = 1/F2 + g0^2/G2;
u2b = 1/G2;
uab = -g0/G2;

Step 7.

r = w.*(y - a - b*x);

Step 8.

chi_sq_obs = sum(r.*r);
nu = m - 2;

Step 9.

chi_sq = calc_chi_sq_95_percent_quantile(nu);

Display information on screen and generate figures

Measurement model.

fprintf('\nMODEL FOR UNCERTAINTIES ASSOCIATED WITH THE YI \n\n')
fprintf('ISO/TS 28037:2010(E) CLAUSE 6 \n')
fprintf('EXAMPLE (UNEQUAL WEIGHTS) \n\n')
MODEL FOR UNCERTAINTIES ASSOCIATED WITH THE YI 

ISO/TS 28037:2010(E) CLAUSE 6 
EXAMPLE (UNEQUAL WEIGHTS) 

Measurement data.

fprintf('FITTING \n')
fprintf(['Data representing ', num2str(m),' measurement points, unequal weights \n'])
for i = 1:m
  fprintf('%8.1f %8.1f %8.1f \n', [x(i), y(i), uy(i)])
end
fprintf('\n')
FITTING 
Data representing 6 measurement points, unequal weights 
     1.0      3.2      0.5 
     2.0      4.3      0.5 
     3.0      7.6      0.5 
     4.0      8.6      1.0 
     5.0     11.7      1.0 
     6.0     12.8      1.0 

Calculation tableau: see write_wls_tableau.m.

write_wls_tableau(x, y, w, g0, h0, g, h, a, b, r, '%8.3f ');
Calculation tableau 
   0.000    0.000    0.000    0.000    2.600    6.233    0.000    0.000    0.885    0.000 
   2.000    4.000    4.000   12.800   -3.200   -6.067   10.240   19.413    0.516    0.266 
   2.000    4.000    8.000   17.200   -1.200   -3.867    1.440    4.640   -1.398    1.955 
   2.000    4.000   12.000   30.400    0.800    2.733    0.640    2.187    1.088    1.183 
   1.000    1.000    4.000    8.600    1.400    2.367    1.960    3.313   -0.513    0.263 
   1.000    1.000    5.000   11.700    2.400    5.467    5.760   13.120    0.530    0.281 
   1.000    1.000    6.000   12.800    3.400    6.567   11.560   22.327   -0.427    0.182 
   0.000   15.000   39.000   93.500    0.000    0.000   31.600   65.000    2.057    4.131 

Solution estimates.

fprintf('Estimate of intercept \n'), fprintf('%8.3f \n\n', a)
fprintf('Estimate of slope \n'), fprintf('%8.3f \n\n', b)
Estimate of intercept 
   0.885 

Estimate of slope 
   2.057 

Standard uncertainties associated with solution estimates.

fprintf('Standard uncertainty associated with estimate of intercept \n'), fprintf('%8.3f \n\n', sqrt(u2a))
fprintf('Standard uncertainty associated with estimate of slope \n'), fprintf('%8.3f \n\n', sqrt(u2b))
fprintf('Covariance associated with estimates of intercept and slope \n'), fprintf('%8.3f \n\n', uab)
Standard uncertainty associated with estimate of intercept 
   0.530 

Standard uncertainty associated with estimate of slope 
   0.178 

Covariance associated with estimates of intercept and slope 
  -0.082 

Validation of the model.

fprintf('VALIDATION \n')
fprintf('Degrees of freedom \n'), fprintf('%4u \n\n', nu)
fprintf('Observed chi-squared value \n'), fprintf('%9.3f \n\n', chi_sq_obs)
fprintf('95 %% quantile of chi-squared distribution with %u degrees of freedom', nu), fprintf('\n%9.3f \n\n', chi_sq)
if chi_sq_obs > chi_sq
  fprintf('CHI-SQUARED TEST FAILED - STRAIGHT-LINE MODEL IS REJECTED \n\n')
else
  fprintf('CHI-SQUARED TEST PASSED - STRAIGHT-LINE MODEL IS ACCEPTED \n\n')
end
VALIDATION 
Degrees of freedom 
   4 

Observed chi-squared value 
    4.131 

95 % quantile of chi-squared distribution with 4 degrees of freedom
    9.488 

CHI-SQUARED TEST PASSED - STRAIGHT-LINE MODEL IS ACCEPTED 

Figures.

set(0, 'DefaultLineLineWidth', 2)
set(0, 'DefaultAxesFontSize', 12)
set(0, 'DefaultAxesFontWeight', 'bold')
figure, hold on
errorbar(x, y, uy, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 6)
plot(x, a + b*x, 'b-')
xlabel('\it x', 'FontName', 'Times', 'FontSize', 14)
ylabel('\it y', 'FontName', 'Times', 'FontSize', 14)
axis1 = axis;
figure, hold on
for i = 1:m
  plot([x(i), x(i)], [0, r(i)], 'k-')
  plot(x(i), r(i), 'ko', 'MarkerFaceColor', 'w', 'MarkerSize', 6)
end
plot([axis1(1), axis1(2)], [0, 0], 'b--')
xlabel('\it x', 'FontName', 'Times', 'FontSize', 14)
ylabel('\it r', 'FontName', 'Times', 'FontSize', 14)

Prediction

yp = 10.5;
uyp = 1.0;
xp = (yp - a)/b;
ca = -1/b;
cb = -(yp - a)/(b^2);
cy = 1/b;
u2xp = ca^2*u2a + cb^2*u2b + 2*ca*cb*uab + cy^2*uyp^2;
fprintf('ISO/TS 28037:2010(E) 11.1 \n')
fprintf('EXAMPLE 2 \n\n')
fprintf('PREDICTION \n')
fprintf('Measured y-value \n'), fprintf('%8.3f \n\n', yp)
fprintf('Standard uncertainty associated with measured y-value \n'), fprintf('%8.3f \n\n', uyp)
fprintf('Estimate of x \n'), fprintf('%8.3f \n\n', xp)
fprintf('Sensitivity coefficient wrt a \n'), fprintf('%8.3f \n\n', ca)
fprintf('Sensitivity coefficient wrt b \n'), fprintf('%8.3f \n\n', cb)
fprintf('Sensitivity coefficient wrt y \n'), fprintf('%8.3f \n\n', cy)
fprintf('Standard uncertainty associated with estimate of x \n'), fprintf('%8.3f \n\n', sqrt(u2xp))
ISO/TS 28037:2010(E) 11.1 
EXAMPLE 2 

PREDICTION 
Measured y-value 
  10.500 

Standard uncertainty associated with measured y-value 
   1.000 

Estimate of x 
   4.674 

Sensitivity coefficient wrt a 
  -0.486 

Sensitivity coefficient wrt b 
  -2.272 

Sensitivity coefficient wrt y 
   0.486 

Standard uncertainty associated with estimate of x 
   0.533 

Acknowledgements

The work described here was supported by the National Measurement System (NMS) delivered through the UK Government’s Department for Science, Innovation and Technology (DSIT).