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

Contents

Introduction

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

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.3 5.6 7.1 9.3 10.7 12.1]';

Assign uncertainties associated with y-values.

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

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 (EQUAL WEIGHTS) \n\n')
MODEL FOR UNCERTAINTIES ASSOCIATED WITH THE YI 

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

Measurement data.

fprintf('FITTING \n')
fprintf(['Data representing ', num2str(m),' measurement points, equal 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, equal weights 
     1.0      3.3      0.5 
     2.0      5.6      0.5 
     3.0      7.1      0.5 
     4.0      9.3      0.5 
     5.0     10.7      0.5 
     6.0     12.1      0.5 

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    3.500    8.017    0.000    0.000    1.867    0.000 
   2.000    4.000    4.000   13.200   -5.000   -9.433   25.000   47.167   -0.648    0.419 
   2.000    4.000    8.000   22.400   -3.000   -4.833    9.000   14.500    0.438    0.192 
   2.000    4.000   12.000   28.400   -1.000   -1.833    1.000    1.833   -0.076    0.006 
   2.000    4.000   16.000   37.200    1.000    2.567    1.000    2.567    0.810    0.655 
   2.000    4.000   20.000   42.800    3.000    5.367    9.000   16.100    0.095    0.009 
   2.000    4.000   24.000   48.400    5.000    8.167   25.000   40.833   -0.619    0.383 
   0.000   24.000   84.000  192.400    0.000    0.000   70.000  123.000    1.757    1.665 

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 
   1.867 

Estimate of slope 
   1.757 

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.465 

Standard uncertainty associated with estimate of slope 
   0.120 

Covariance associated with estimates of intercept and slope 
  -0.050 

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 
    1.665 

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 = 0.5;
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 1 \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 1 

PREDICTION 
Measured y-value 
  10.500 

Standard uncertainty associated with measured y-value 
   0.500 

Estimate of x 
   4.913 

Sensitivity coefficient wrt a 
  -0.569 

Sensitivity coefficient wrt b 
  -2.796 

Sensitivity coefficient wrt y 
   0.569 

Standard uncertainty associated with estimate of x 
   0.322 

Forward evaluation

xf = 3.5;
uxf = 0.2;
yf = a + b*xf;
ca = 1;
cb = xf;
cx = b;
u2yf = ca^2*u2a + cb^2*u2b + 2*ca*cb*uab + cx^2*uxf^2;
fprintf('ISO/TS 28037:2010(E) 11.2 \n')
fprintf('EXAMPLE \n\n')
fprintf('FORWARD EVALUATION \n')
fprintf('Measured x-value \n'), fprintf('%8.3f \n\n', xf)
fprintf('Standard uncertainty associated with measured x-value \n'), fprintf('%8.3f \n\n', uxf)
fprintf('Estimate of y \n'), fprintf('%8.3f \n\n', yf)
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 x \n'), fprintf('%8.3f \n\n', cx)
fprintf('Standard uncertainty associated with estimate of y \n'), fprintf('%8.3f \n\n', sqrt(u2yf))
ISO/TS 28037:2010(E) 11.2 
EXAMPLE 

FORWARD EVALUATION 
Measured x-value 
   3.500 

Standard uncertainty associated with measured x-value 
   0.200 

Estimate of y 
   8.017 

Sensitivity coefficient wrt a 
   1.000 

Sensitivity coefficient wrt b 
   3.500 

Sensitivity coefficient wrt x 
   1.757 

Standard uncertainty associated with estimate of y 
   0.406 

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).