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