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

Contents

Introduction

This document runs the numerical example of Gauss-Markov regression (GMR) described in Clause 9 (Model for uncertainties and covariances associated with the y_i).

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 7.0 8.0 9.0 10.0]';
m = length(x);

Assign y-values.

y = [1.3 4.1 6.9 7.5 10.2 12.0 14.5 17.1 19.5 21.0]';

Assign covariance matrix associated with y-values.

Uy = [2  1  1  1  1  0  0  0  0  0
      1  2  1  1  1  0  0  0  0  0
      1  1  2  1  1  0  0  0  0  0
      1  1  1  2  1  0  0  0  0  0
      1  1  1  1  2  0  0  0  0  0
      0  0  0  0  0  5  4  4  4  4
      0  0  0  0  0  4  5  4  4  4
      0  0  0  0  0  4  4  5  4  4
      0  0  0  0  0  4  4  4  5  4
      0  0  0  0  0  4  4  4  4  5];

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

Solve the Gauss-Markov regression problem to obtain best fit straight-line parameters.

Step 1.

Ly = chol(Uy, 'lower');

Step 2.

e = ones(m, 1);
f = Ly\e;
g = Ly\x;
h = Ly\y;

Step 3.

F2 = sum(f.*f);

Step 4.

g0 = sum(f.*g)/F2;
h0 = sum(f.*h)/F2;

Step 5.

gt = g - g0*f;
ht = h - h0*f;

Step 6.

Gt2 = sum(gt.*gt);

Step 7.

b = sum(gt.*ht)/Gt2;
a = h0 - b*g0;

Step 8.

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

Step 9.

r = ht - b*gt;

Step 10.

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

Step 11.

chi_sq = calc_chi_sq_95_percent_quantile(nu);

Display information on screen and generate figures

Measurement model.

fprintf('\nMODEL FOR UNCERTAINTIES AND COVARIANCES ASSOCIATED WITH THE THE YI \n\n')
fprintf('ISO/TS 28037:2010(E) CLAUSE 9 \n')
fprintf('EXAMPLE \n\n')
MODEL FOR UNCERTAINTIES AND COVARIANCES ASSOCIATED WITH THE THE YI 

ISO/TS 28037:2010(E) CLAUSE 9 
EXAMPLE 

Measurement data.

fprintf('FITTING \n')
fprintf(['Data representing ', num2str(m),' measurement points \n'])
for i = 1:m
  fprintf('%8.1f %8.1f \n', [x(i), y(i)])
end
fprintf('\n')
FITTING 
Data representing 10 measurement points 
     1.0      1.3 
     2.0      4.1 
     3.0      6.9 
     4.0      7.5 
     5.0     10.2 
     6.0     12.0 
     7.0     14.5 
     8.0     17.1 
     9.0     19.5 
    10.0     21.0 

Covariance matrix Uy.

fprintf('Covariance matrix Uy \n')
disp(Uy)
Covariance matrix Uy 
     2     1     1     1     1     0     0     0     0     0
     1     2     1     1     1     0     0     0     0     0
     1     1     2     1     1     0     0     0     0     0
     1     1     1     2     1     0     0     0     0     0
     1     1     1     1     2     0     0     0     0     0
     0     0     0     0     0     5     4     4     4     4
     0     0     0     0     0     4     5     4     4     4
     0     0     0     0     0     4     4     5     4     4
     0     0     0     0     0     4     4     4     5     4
     0     0     0     0     0     4     4     4     4     5

Cholesky factor Ly of Uy.

fprintf('Cholesky factor Ly of Uy \n')
disp(Ly)
Cholesky factor Ly of Uy 
  Columns 1 through 7

    1.4142         0         0         0         0         0         0
    0.7071    1.2247         0         0         0         0         0
    0.7071    0.4082    1.1547         0         0         0         0
    0.7071    0.4082    0.2887    1.1180         0         0         0
    0.7071    0.4082    0.2887    0.2236    1.0954         0         0
         0         0         0         0         0    2.2361         0
         0         0         0         0         0    1.7889    1.3416
         0         0         0         0         0    1.7889    0.5963
         0         0         0         0         0    1.7889    0.5963
         0         0         0         0         0    1.7889    0.5963

  Columns 8 through 10

         0         0         0
         0         0         0
         0         0         0
         0         0         0
         0         0         0
         0         0         0
         0         0         0
    1.2019         0         0
    0.3698    1.1435         0
    0.3698    0.2691    1.1114

Calculation tableau: see write_gmr_tableaux.m.

write_gmr_tableaux(f, g, h, g0, h0, gt, ht, a, b, r, '%9.4f ');
Initial calculation tableau 
   0.7071    0.7071    0.9192 
   0.4082    1.2247    2.8169 
   0.2887    1.7321    4.4167 
   0.2236    2.2361    3.9578 
   0.1826    2.7386    5.6963 
   0.4472    2.6833    5.3666 
   0.1491    1.6398    3.6522 
   0.0925    1.8490    4.4284 
   0.0673    2.2198    5.3208 
   0.0529    2.6463    5.5360 

Main calculation tableau 
   0.0000    0.0000    0.0000    4.1111    8.4044    0.0000    0.0000   -0.6456    0.0000 
   0.5000    0.5000    0.6500   -2.1999   -5.0236    4.8395   11.0514   -0.1809    0.0327 
   0.1667    0.5000    1.1500   -0.4536   -0.6142    0.2058    0.2786    0.3844    0.1477 
   0.0833    0.5000    1.2750    0.5453    1.9906    0.2973    1.0854    0.7902    0.6245 
   0.0500    0.5000    0.8850    1.3168    2.0785    1.7340    2.7370   -0.8202    0.6727 
   0.0333    0.5000    1.0400    1.9880    4.1619    3.9523    8.2739   -0.2145    0.0460 
   0.2000    1.2000    2.4000    0.8447    1.6080    0.7136    1.3583   -0.2516    0.0633 
   0.0222    0.2444    0.5444    1.0269    2.3994    1.0546    2.4640    0.1387    0.0192 
   0.0085    0.1709    0.4094    1.4689    3.6514    2.1578    5.3636    0.4177    0.1745 
   0.0045    0.1493    0.3579    1.9433    4.7555    3.7763    9.2412    0.4777    0.2282 
   0.0028    0.1401    0.2930    2.4287    5.0912    5.8986   12.3650   -0.2552    0.0651 
   1.0714    4.4048    9.0048    0.0000    0.0000   24.6296   54.2185    2.2014    2.0740 

Solution estimates.

fprintf('Estimate of intercept \n'), fprintf('%9.4f \n\n', a)
fprintf('Estimate of slope \n'), fprintf('%9.4f \n\n', b)
Estimate of intercept 
  -0.6456 

Estimate of slope 
   2.2014 

Standard uncertainties associated with solution estimates.

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

Standard uncertainty associated with estimate of slope 
   0.2015 

Covariance associated with estimates of intercept and slope 
  -0.1669 

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 
   8 

Observed chi-squared value 
    2.074 

95 % quantile of chi-squared distribution with 8 degrees of freedom
   15.507 

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, sqrt(diag(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)

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