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

Contents

Introduction

This document runs a numerical example to illustrate the algorithm of generalized distance regression (GDR) described in Clause 8 (Model for uncertainties associated with the x_i and the y_i and covariances associated with the pairs (x_i, 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.2 1.9 2.9 4.0 4.7 5.9]';
m = length(x);

Assign y-values.

y = [3.4 4.4 7.2 8.5 10.8 13.5]';

Assign uncertainties associated with x-values.

ux = [0.2 0.2 0.2 0.2 0.2 0.2]';

Assign uncertainties associated with y-values.

uy = [0.2 0.2 0.2 0.4 0.4 0.4]';

Assign covariances associated with (x, y)-pairs of values.

uxy = [0.036 0.036 0.036 0.036 0.036 0.036]';

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

Solve the generalized distance 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, u2ai, u2bi, uabi, wi, g0i, h0i, gi, hi, G2i, ri, Ri] = algm_wls_steps_1_to_8(x, y, uy);
at{1} = ai;
bt{1} = bi;

Loop through steps 2 to 5 until convergence criteria are satisfied. (In this example, convergence is considered to have been achieved when the magnitudes of the increments deltaa and deltab in a and b are no greater than 0.00005. 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 = 0.00005;
da{1} = []; db{1} = [];
t{1} = []; xs{1} = []; z{1} = []; f{1} = []; g{1} = []; h{1} = []; F2{1} = [];
g0{1} = []; h0{1} = []; gt{1} = []; ht{1} = []; Gt2{1} = []; r{1} = [];

Steps 2 to 5: see algm_gdr2_steps_2_to_5.m.

ind = 1;
[at, bt, da, db, t, xs, z, f, g, h, F2, g0, h0, gt, ht, Gt2, r] ...
  = algm_gdr2_steps_2_to_5(x, ux, y, uy, uxy, at, bt, da, db, ...
      t, xs, z, f, g, h, F2, g0, h0, gt, ht, Gt2, r, ind);

while (abs(da{ind}) > tol) || (abs(db{ind}) > tol)

Update iteration number.

  ind = ind + 1;

Step 6. Repeat steps 2 to 5 until convergence has been achieved: see algm_gdr2_steps_2_to_5.m.

  [at, bt, da, db, t, xs, z, f, g, h, F2, g0, h0, gt, ht, Gt2, r] ...
    = algm_gdr2_steps_2_to_5(x, ux, y, uy, uxy, at, bt, da, db, ...
        t, xs, z, f, g, h, F2, g0, h0, gt, ht, Gt2, r, ind);
end
a = at{ind+1};
b = bt{ind+1};

Step 7. Evaluate uncertainties.

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

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

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

Step 9. 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 and generate figures

Measurement model.

fprintf('\nMODEL FOR UNCERTAINTIES ASSOCIATED WITH THE XI AND THE YI AND COVARIANCES ASSOCIATED WITH THE PAIRS (XI, YI) \n\n')
MODEL FOR UNCERTAINTIES ASSOCIATED WITH THE XI AND THE YI AND COVARIANCES ASSOCIATED WITH THE PAIRS (XI, YI) 

Measurement data.

fprintf('FITTING \n')
fprintf(['Data representing ', num2str(m),' measurement points \n'])
for i = 1:m
  fprintf('%8.1f %8.1f %8.1f %8.1f %8.3f\n', [x(i), ux(i), y(i), uy(i), uxy(i)])
end
fprintf('\n')
FITTING 
Data representing 6 measurement points 
     1.2      0.2      3.4      0.2    0.036
     1.9      0.2      4.4      0.2    0.036
     2.9      0.2      7.2      0.2    0.036
     4.0      0.2      8.5      0.4    0.036
     4.7      0.2     10.8      0.4    0.036
     5.9      0.2     13.5      0.4    0.036

Calculation tableau for WLS problem: see write_wls_tableau.m.

fprintf('Initial approximations to parameters \n')
write_wls_tableau(x, y, wi, g0i, h0i, gi, hi, ai, bi, ri, '%9.4f ');
Initial approximations to parameters 
Calculation tableau 
   0.0000    0.0000    0.0000    0.0000    2.5733    6.1867    0.0000    0.0000    0.6583    0.0000 
   5.0000   25.0000   30.0000   85.0000   -6.8667  -13.9333   47.1511   95.6756    0.8186    0.6701 
   5.0000   25.0000   47.5000  110.0000   -3.3667   -8.9333   11.3344   30.0756   -1.7006    2.8920 
   5.0000   25.0000   72.5000  180.0000    1.6333    5.0667    2.6678    8.2756    1.5577    2.4264 
   2.5000    6.2500   25.0000   53.1250    3.5667    5.7833   12.7211   20.6272   -1.8791    3.5310 
   2.5000    6.2500   29.3750   67.5000    5.3167   11.5333   28.2669   61.3189    0.1113    0.0124 
   2.5000    6.2500   36.8750   84.3750    8.3167   18.2833   69.1669  152.0564    0.4163    0.1733 
   0.0000   93.7500  241.2500  580.0000    0.0000    0.0000  171.3083  368.0292    2.1483    9.7052 

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

Initial approximation to slope 
   2.1483 

Calculation tableau for each iteration: see write_gdr2_tableaux.m. (The tableau is similar to that for the problem of generalized distance regression (GDR) described in Clause 7, but includes an additional column (column 5) in the first tableau for each iteration that contains values of the covariances uxy.)

for niter = 1:ind
  write_gdr2_tableaux(x, ux, y, uy, uxy, t, xs, z, f, g, h, g0, h0, gt, ht, ...
                    at, bt, da, db, r, niter, '%9.4f');
end
Iteration 1 to determine fi, gi, and hi 
   0.0000   0.0000   0.0000   0.0000   0.0000   0.6583   2.1483   0.0000   0.0000   0.0000   0.0000
   1.2000   0.2000   3.4000   0.2000   0.0360  14.2991   1.3169   0.1637   3.7814   4.9797   0.6191
   1.9000   0.2000   4.4000   0.2000   0.0360  14.2991   1.6572  -0.3401   3.7814   6.2664  -1.2861
   2.9000   0.2000   7.2000   0.2000   0.0360  14.2991   3.1224   0.3115   3.7814  11.8073   1.1781
   4.0000   0.2000   8.5000   0.4000   0.0360   5.2650   3.8024  -0.7516   2.2946   8.7248  -1.7247
   4.7000   0.2000  10.8000   0.4000   0.0360   5.2650   4.7117   0.0445   2.2946  10.8113   0.1022
   5.9000   0.2000  13.5000   0.4000   0.0360   5.2650   5.9438   0.1665   2.2946  13.6383   0.3821
 
Iteration 1 to determine increments deltaa and deltab 
   0.0000   0.0000   0.0000   2.7822  -0.0156   0.0000   0.0000  -0.0751   0.0000
  14.2991  18.8305   2.3411  -5.5410   0.6780  30.7025  -3.7567   0.7966   0.6346
  14.2991  23.6958  -4.8633  -4.2543  -1.2272  18.0994   5.2211  -1.1362   1.2909
  14.2991  44.6481   4.4548   1.2865   1.2369   1.6551   1.5913   1.2094   1.4626
   5.2650  20.0195  -3.9573   2.3408  -1.6889   5.4795  -3.9536  -1.7391   3.0243
   5.2650  24.8070   0.2344   4.4273   0.1379  19.6010   0.6104   0.0431   0.0019
   5.2650  31.2938   0.8767   7.2544   0.4178  52.6256   3.0308   0.2625   0.0689
  58.6922 163.2948  -0.9137   0.0000   0.0000 128.1633   2.7434   0.0214   6.4832

Iteration 2 to determine fi, gi, and hi 
   0.0000   0.0000   0.0000   0.0000   0.0000   0.5831   2.1697   0.0000   0.0000   0.0000   0.0000
   1.2000   0.2000   3.4000   0.2000   0.0360  13.8714   1.3502   0.2132   3.7244   5.0287   0.7939
   1.9000   0.2000   4.4000   0.2000   0.0360  13.8714   1.6847  -0.3057   3.7244   6.2744  -1.1384
   2.9000   0.2000   7.2000   0.2000   0.0360  13.8714   3.1287   0.3246   3.7244  11.6526   1.2089
   4.0000   0.2000   8.5000   0.4000   0.0360   5.2059   3.7985  -0.7621   2.2816   8.6668  -1.7389
   4.7000   0.2000  10.8000   0.4000   0.0360   5.2059   4.7050   0.0190   2.2816  10.7352   0.0434
   5.9000   0.2000  13.5000   0.4000   0.0360   5.2059   5.9305   0.1153   2.2816  13.5313   0.2632
 
Iteration 2 to determine increments deltaa and deltab 
   0.0000   0.0000   0.0000   2.8068  -0.0008   0.0000   0.0000  -0.0015   0.0000
  13.8714  18.7289   2.9569  -5.4251   0.7971  29.4316  -4.3242   0.7984   0.6374
  13.8714  23.3685  -4.2400  -4.1794  -1.1353  17.4671   4.7446  -1.1343   1.2865
  13.8714  43.3993   4.5025   1.1988   1.2121   1.4372   1.4531   1.2118   1.4684
   5.2059  19.7745  -3.9676   2.2627  -1.7370   5.1197  -3.9302  -1.7375   3.0190
   5.2059  24.4938   0.0991   4.3311   0.0454  18.7583   0.1965   0.0443   0.0020
   5.2059  30.8734   0.6004   7.1272   0.2651  50.7963   1.8894   0.2634   0.0694
  57.2320 160.6386  -0.0486   0.0000   0.0000 123.0102   0.0293   0.0002   6.4827

Iteration 3 to determine fi, gi, and hi 
   0.0000   0.0000   0.0000   0.0000   0.0000   0.5816   2.1700   0.0000   0.0000   0.0000   0.0000
   1.2000   0.2000   3.4000   0.2000   0.0360  13.8668   1.3510   0.2144   3.7238   5.0310   0.7984
   1.9000   0.2000   4.4000   0.2000   0.0360  13.8668   1.6854  -0.3046   3.7238   6.2762  -1.1343
   2.9000   0.2000   7.2000   0.2000   0.0360  13.8668   3.1292   0.3254   3.7238  11.6527   1.2118
   4.0000   0.2000   8.5000   0.4000   0.0360   5.2052   3.7986  -0.7616   2.2815   8.6665  -1.7375
   4.7000   0.2000  10.8000   0.4000   0.0360   5.2052   4.7051   0.0194   2.2815  10.7348   0.0443
   5.9000   0.2000  13.5000   0.4000   0.0360   5.2052   5.9305   0.1155   2.2815  13.5305   0.2634
 
Iteration 3 to determine increments deltaa and deltab 
   0.0000   0.0000   0.0000   2.8075  -0.0000   0.0000   0.0000  -0.0000   0.0000
  13.8668  18.7344   2.9730  -5.4235   0.7984  29.4145  -4.3301   0.7984   0.6375
  13.8668  23.3716  -4.2238  -4.1782  -1.1342  17.4577   4.7391  -1.1342   1.2864
  13.8668  43.3924   4.5125   1.1982   1.2118   1.4357   1.4520   1.2118   1.4685
   5.2052  19.7727  -3.9641   2.2613  -1.7375   5.1136  -3.9291  -1.7375   3.0190
   5.2052  24.4913   0.1012   4.3295   0.0444  18.7449   0.1921   0.0443   0.0020
   5.2052  30.8697   0.6010   7.1253   0.2634  50.7693   1.8769   0.2634   0.0694
  57.2160 160.6321  -0.0003   0.0000   0.0000 122.9356   0.0010   0.0000   6.4827

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

Estimate of slope 
   2.1700 

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 
   0.2856 

Standard uncertainty associated with estimate of slope 
   0.0902 

Covariance associated with estimates of intercept and slope 
  -0.0228 

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 
    6.483 

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
plot_errorbar2(x, ux, y, uy, 'k.', 'k-')
plot(xs{ind}, at{ind} + bt{ind}*xs{ind}, 'ro');
plot(xs{ind}, at{ind} + bt{ind}*xs{ind}, '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{ind}(i)], 'k-')
  plot(x(i), r{ind}(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).