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

Contents

Introduction

This document runs the numerical example of generalized Gauss-Markov regression (GGMR) described in Clause 10 (Model for uncertainties and covariances associated with the x_i and the y_i) and Annex C EXAMPLE 1 using the orthogonal factorization approach described in C.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 = [50.4 99.0 149.9 200.4 248.5 299.7 349.1]';
m = length(x);

Assign y-values.

y = [52.3 97.8 149.7 200.1 250.4 300.9 349.2]';

Assign factor Bx of covariance matrix Ux associated with x-values.

Bx = [0.5  0    0    0    0    0    0    0.5  0  0
      0    0.5  0    0    0    0    0    0    1  0
      0    0    0.5  0    0    0    0    0.5  1  0
      0    0    0    0.5  0    0    0    0    0  1
      0    0    0    0    0.5  0    0    0.5  0  1
      0    0    0    0    0    0.5  0    0    1  1
      0    0    0    0    0    0    0.5  0.5  1  1];

Assign factor By of covariance matrix Uy associated with y-values.

By = [2  0  0  0  0  0  0  1
      0  2  0  0  0  0  0  1
      0  0  2  0  0  0  0  1
      0  0  0  2  0  0  0  1
      0  0  0  0  2  0  0  1
      0  0  0  0  0  2  0  1
      0  0  0  0  0  0  2  1];

Build complete factor for x- and y-values (assuming the correlation associated with each (x, y)-pair is zero).

B = blkdiag(Bx, By);

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

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

Step 1. Initial approximation using weighted least squares: see algm_wls_steps_1_to_8.m.

[ai, bi] = algm_wls_steps_1_to_8(x, y, sqrt(diag(By*By')));

Round approximations to parameters in step 1 to four decimal places. (This step is included to produce the results given in ISO/TS 28037:2010: the step would not generally be performed.)

ai = round(10000*ai)/10000;
at{1} = ai;
bi = round(10000*bi)/10000;
bt{1} = bi;

Initial approximation.

tt{1} = [x; at{1}; bt{1}];

Loop through steps 2 to 9 until convergence criteria are satisfied. (In this example, convergence is considered to have been achieved when the magnitudes of all increments are no greater than 1e-7. For general user data, it is suggested that convergence is considered to have been achieved when the magnitudes of all increments relative to the initial approximations to the straight-line parameters are no greater than 1e-7. In this case, the tolerance can be assigned using the command tol = 1e-7*norm([ai, bi]);)

Assign tolerances and initialize variables.

tol = 1e-7;
dt{1} = []; f{1} = []; J{1} = []; Q{1} = []; R10{1} = []; R1{1} = [];
T{1} = []; Z{1} = []; ft{1} = []; ft1{1} = []; ft2{1} = []; T11{1} = [];
T12{1} = []; T22{1} = []; et2{1} = [];

Steps 2 to 8: see algm_ggmr_orthogonal_steps_2_to_8.m.

ind = 1;
[tt, dt, f, J, Q, R1, T, Z, ft1, ft2, T11, T12, T22, et2] = ...
  algm_ggmr_orthogonal_steps_2_to_8(x, y, B, tt, dt, f, J, Q, R1, T, Z, ...
    ft1, ft2, T11, T12, T22, et2, ind);

while any(abs(dt{ind}) > tol)

Update iteration number.

  ind = ind + 1;

Step 9. Repeat steps 2 to 8 until convergence has been achieved: see algm_ggmr_orthogonal_steps_2_to_8.m.

  [tt, dt, f, J, Q, R1, T, Z, ft1, ft2, T11, T12, T22, et2] = ...
    algm_ggmr_orthogonal_steps_2_to_8(x, y, B, tt, dt, f, J, Q, R1, T, Z, ...
      ft1, ft2, T11, T12, T22, et2, ind);
end
a = tt{ind+1}(m+1);
b = tt{ind+1}(m+2);

Step 10. Solve upper-triangular system and evaluate uncertainties.

Ra = R1{ind}(m+1:m+2, m+1:m+2);
Ta = T11{ind}(end-1:end, end-1:end);
Ka = Ra\Ta;
Ua = Ka*Ka';
u2a = Ua(1, 1);
u2b = Ua(2, 2);
uab = Ua(1, 2);

Step 11. Form observed chi-squared value and degrees of freedom.

chi_sq_obs = sum(et2{ind}.*et2{ind});
nu = m - 2;

Step 12. Calculate 95 % quantile of chi-squared distribution: see calc_chi_sq_95_percent_quantile.m.

chi_sq = calc_chi_sq_95_percent_quantile(nu);

Display information on screen

Measurement model.

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

ISO/TS 28037:2010(E) ANNEX C 
EXAMPLE 1 

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 7 measurement points 
    50.4     52.3 
    99.0     97.8 
   149.9    149.7 
   200.4    200.1 
   248.5    250.4 
   299.7    300.9 
   349.1    349.2 

Covariance matrix Ux.

fprintf('Covariance matrix Ux \n')
disp(Bx*Bx')
Covariance matrix Ux 
    0.5000         0    0.2500         0    0.2500         0    0.2500
         0    1.2500    1.0000         0         0    1.0000    1.0000
    0.2500    1.0000    1.5000         0    0.2500    1.0000    1.2500
         0         0         0    1.2500    1.0000    1.0000    1.0000
    0.2500         0    0.2500    1.0000    1.5000    1.0000    1.2500
         0    1.0000    1.0000    1.0000    1.0000    2.2500    2.0000
    0.2500    1.0000    1.2500    1.0000    1.2500    2.0000    2.5000

Factor Bx of Ux.

fprintf('Factor Bx of Ux \n')
disp(Bx)
Factor Bx of Ux 
  Columns 1 through 7

    0.5000         0         0         0         0         0         0
         0    0.5000         0         0         0         0         0
         0         0    0.5000         0         0         0         0
         0         0         0    0.5000         0         0         0
         0         0         0         0    0.5000         0         0
         0         0         0         0         0    0.5000         0
         0         0         0         0         0         0    0.5000

  Columns 8 through 10

    0.5000         0         0
         0    1.0000         0
    0.5000    1.0000         0
         0         0    1.0000
    0.5000         0    1.0000
         0    1.0000    1.0000
    0.5000    1.0000    1.0000

Covariance matrix Uy.

fprintf('Covariance matrix Uy \n')
disp(By*By')
Covariance matrix Uy 
     5     1     1     1     1     1     1
     1     5     1     1     1     1     1
     1     1     5     1     1     1     1
     1     1     1     5     1     1     1
     1     1     1     1     5     1     1
     1     1     1     1     1     5     1
     1     1     1     1     1     1     5

Factor By of Uy.

fprintf('Factor By of Uy \n')
disp(By)
Factor By of Uy 
     2     0     0     0     0     0     0     1
     0     2     0     0     0     0     0     1
     0     0     2     0     0     0     0     1
     0     0     0     2     0     0     0     1
     0     0     0     0     2     0     0     1
     0     0     0     0     0     2     0     1
     0     0     0     0     0     0     2     1

Initial approximations to parameters.

fprintf('Initial approximation to intercept \n'), fprintf('%9.4f \n\n', ai)
fprintf('Initial approximation to slope \n'), fprintf('%9.4f \n\n', bi)
Initial approximation to intercept 
   0.2707 

Initial approximation to slope 
   1.0011 

Summary of iterative procedure: see write_ggmr_tableau.m.

write_ggmr_tableau(tt, dt, ind, '%14.4e ');
Summary of iterative procedure (initial approximations, corrections and final estimates) 
    5.0400e+01     1.7253e-01     1.2580e-04     3.0782e-06     2.9044e-09     5.0573e+01 
    9.9000e+01    -4.3150e-01    -3.2145e-04    -6.3201e-06    -7.1009e-09     9.8568e+01 
    1.4990e+02    -2.9164e-01    -3.9604e-04    -3.8889e-06    -7.5640e-09     1.4961e+02 
    2.0040e+02     2.9677e-02    -1.0763e-03    -6.0240e-07    -1.7165e-08     2.0043e+02 
    2.4850e+02     2.4039e-01    -1.1406e-03     3.2378e-06    -1.7064e-08     2.4874e+02 
    2.9970e+02    -2.2251e-01    -1.5777e-03    -3.3581e-06    -2.6110e-08     2.9948e+02 
    3.4910e+02    -2.0619e-01    -1.6622e-03    -3.3805e-06    -2.7429e-08     3.4889e+02 
    2.7070e-01     7.5040e-02    -3.3396e-03     1.0058e-07    -5.3018e-08     3.4240e-01 
    1.0011e+00     1.0962e-04     2.1135e-05     7.6429e-09     3.3717e-10     1.0012e+00 

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

Estimate of slope 
   1.0012 

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 
   2.0569 

Standard uncertainty associated with estimate of slope 
   0.0090 

Covariance associated with estimates of intercept and slope 
  -0.0129 

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 
   5 

Observed chi-squared value 
    1.772 

95 % quantile of chi-squared distribution with 5 degrees of freedom
   11.070 

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

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