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

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 covariance matrix associated with x-values.

Ux = [0.5   0     0.25  0     0.25  0     0.25
      0     1.25  1     0     0     1     1
      0.25  1     1.5   0     0.25  1     1.25
      0     0     0     1.25  1     1     1
      0.25  0     0.25  1     1.5   1     1.25
      0     1     1     1     1     2.25  2
      0.25  1     1.25  1     1.25  2     2.5];

Assign covariance matrix associated with y-values.

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];

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

U = blkdiag(Ux, Uy);

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(Uy)));

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} = []; ft{1} = []; Jt{1} = []; g{1} = [];
H{1} = []; M{1} = []; q{1} = [];

Steps 2 to 9: see algm_ggmr_cholesky_steps_2_to_9.m.

ind = 1;
[tt, dt, f, J, ft, Jt, g, H, M, q] ...
  = algm_ggmr_cholesky_steps_2_to_9(x, y, U, tt, dt, f, J, ft, Jt, g, H, M, q, ind);

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

Update iteration number.

  ind = ind + 1;

Step 10. Repeat steps 2 to 9 until convergence has been achieved: see algm_ggmr_cholesky_steps_2_to_9.m.

  [tt, dt, f, J, ft, Jt, g, H, M, q] ...
    = algm_ggmr_cholesky_steps_2_to_9(x, y, U, tt, dt, f, J, ft, Jt, g, H, M, q, ind);
end
a = tt{ind+1}(m+1);
b = tt{ind+1}(m+2);

Step 11. Partition Cholesky factor and evaluate uncertainties.

M22 = M{ind}(m+1:m+2, m+1:m+2);
m11 = M22(1, 1);
m21 = M22(2, 1);
m22 = M22(2, 2);
u2a = (m22^2 + m21^2)/(m11^2*m22^2);
u2b = 1/(m22^2);
uab = -m21/(m11*m22^2);

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

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

Step 13. 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) CLAUSE 10 \n')
fprintf('EXAMPLE \n\n')
MODEL FOR UNCERTAINTIES AND COVARIANCES ASSOCIATED WITH THE XI AND THE YI 

ISO/TS 28037:2010(E) CLAUSE 10 
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 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(Ux)
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

Cholesky factor Lx of Ux.

fprintf('Cholesky factor Lx of Ux \n')
disp(chol(Ux, 'lower'))
Cholesky factor Lx of Ux 
    0.7071         0         0         0         0         0         0
         0    1.1180         0         0         0         0         0
    0.3536    0.8944    0.7583         0         0         0         0
         0         0         0    1.1180         0         0         0
    0.3536         0    0.1648    0.8944    0.7402         0         0
         0    0.8944    0.2638    0.8944    0.2115    0.7319         0
    0.3536    0.8944    0.4286    0.8944    0.3436    0.2928    0.6225

Covariance matrix Uy.

fprintf('Covariance matrix Uy \n')
disp(Uy)
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

Cholesky factor Ly of Uy.

fprintf('Cholesky factor Ly of Uy \n')
disp(chol(Uy, 'lower'))
Cholesky factor Ly of Uy 
    2.2361         0         0         0         0         0         0
    0.4472    2.1909         0         0         0         0         0
    0.4472    0.3651    2.1602         0         0         0         0
    0.4472    0.3651    0.3086    2.1381         0         0         0
    0.4472    0.3651    0.3086    0.2673    2.1213         0         0
    0.4472    0.3651    0.3086    0.2673    0.2357    2.1082         0
    0.4472    0.3651    0.3086    0.2673    0.2357    0.2108    2.0976

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.3019e-08     3.4240e-01 
    1.0011e+00     1.0962e-04     2.1135e-05     7.6429e-09     3.3717e-10     1.0012e+00 

Matrices M and M22 at final iteration.

fprintf('Matrix M at final iteration \n')
disp(M{ind})
fprintf('Matrix M22 at final iteration \n')
disp(M22)
Matrix M at final iteration 
  Columns 1 through 7

    1.7755         0         0         0         0         0         0
    0.1810    1.6959         0         0         0         0         0
   -0.4246   -0.5430    1.4306         0         0         0         0
    0.1810    0.2378    0.6893    1.5312         0         0         0
   -0.4246    0.5053    0.0999   -0.7249    1.2453         0         0
    0.3748   -0.5607   -0.2377   -0.4269   -0.0306    1.3465         0
   -0.2308   -0.2931   -0.8271    0.0932   -0.5828   -0.9711    0.8329
    0.0513    0.0482    0.0971    0.0022    0.0645    0.0927    0.3899
  -10.7652   -3.0381   -0.3815   13.9302   30.1844   38.8412  127.1155

  Columns 8 through 9

         0         0
         0         0
         0         0
         0         0
         0         0
         0         0
         0         0
    0.6762         0
  107.2706  110.9677

Matrix M22 at final iteration 
    0.6762         0
  107.2706  110.9677

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