Software to support ISO/TS 28037:2010(E)
Contents
Introduction
This document runs the numerical example of generalized distance regression (GDR) described in Clause 7 (Model for uncertainties 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 = [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]';
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);
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;
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_gdr1_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_gdr1_steps_2_to_5(x, ux, y, uy, 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_gdr1_steps_2_to_5.m.
[at, bt, da, db, t, xs, z, f, g, h, F2, g0, h0, gt, ht, Gt2, r] ... = algm_gdr1_steps_2_to_5(x, ux, y, uy, 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 \n\n') fprintf('ISO/TS 28037:2010(E) CLAUSE 7 \n') fprintf('EXAMPLE \n\n')
MODEL FOR UNCERTAINTIES ASSOCIATED WITH THE XI AND THE YI ISO/TS 28037:2010(E) CLAUSE 7 EXAMPLE
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 \n', [x(i), ux(i), y(i), uy(i)]) end fprintf('\n')
FITTING Data representing 6 measurement points 1.2 0.2 3.4 0.2 1.9 0.2 4.4 0.2 2.9 0.2 7.2 0.2 4.0 0.2 8.5 0.4 4.7 0.2 10.8 0.4 5.9 0.2 13.5 0.4
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_gdr1_tableaux.m.
for niter = 1:ind write_gdr1_tableaux(x, ux, y, uy, 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.6583 2.1483 0.0000 0.0000 0.0000 0.0000 1.2000 0.2000 3.4000 0.2000 4.4522 1.2626 0.1637 2.1100 2.6642 0.3455 1.9000 0.2000 4.4000 0.2000 4.4522 1.7699 -0.3401 2.1100 3.7345 -0.7176 2.9000 0.2000 7.2000 0.2000 4.4522 3.0192 0.3116 2.1100 6.3706 0.6575 4.0000 0.2000 8.5000 0.4000 2.9019 3.8126 -0.7515 1.7035 6.4947 -1.2802 4.7000 0.2000 10.8000 0.4000 2.9019 4.7111 0.0447 1.7035 8.0253 0.0761 5.9000 0.2000 13.5000 0.4000 2.9019 5.9416 0.1667 1.7035 10.1214 0.2840 Iteration 1 to determine increments deltaa and deltab 0.0000 0.0000 0.0000 3.1239 -0.0437 0.0000 0.0000 -0.0784 0.0000 4.4522 5.6216 0.7290 -3.9273 0.4378 15.4236 -1.7193 0.4814 0.2318 4.4522 7.8799 -1.5141 -2.8570 -0.6253 8.1623 1.7864 -0.5935 0.3523 4.4522 13.4422 1.3874 -0.2209 0.7498 0.0488 -0.1656 0.7523 0.5659 2.9019 11.0636 -2.1807 1.1732 -1.2057 1.3764 -1.4145 -1.2187 1.4852 2.9019 13.6710 0.1297 2.7038 0.1506 7.3108 0.4073 0.1206 0.0145 2.9019 17.2416 0.4838 4.7999 0.3585 23.0387 1.7208 0.3052 0.0931 22.0622 68.9199 -0.9648 0.0000 0.0000 55.3606 0.6152 0.0111 2.7429 Iteration 2 to determine fi, gi, and hi 0.0000 0.0000 0.0000 0.0000 0.5799 2.1594 0.0000 0.0000 0.0000 0.0000 1.2000 0.2000 3.4000 0.2000 4.4146 1.2873 0.2289 2.1011 2.7047 0.4808 1.9000 0.2000 4.4000 0.2000 4.4146 1.7922 -0.2827 2.1011 3.7655 -0.5941 2.9000 0.2000 7.2000 0.2000 4.4146 3.0365 0.3579 2.1011 6.3799 0.7519 4.0000 0.2000 8.5000 0.4000 2.8858 3.8212 -0.7175 1.6988 6.4913 -1.2189 4.7000 0.2000 10.8000 0.4000 2.8858 4.7177 0.0709 1.6988 8.0142 0.1205 5.9000 0.2000 13.5000 0.4000 2.8858 5.9448 0.1796 1.6988 10.0988 0.3051 Iteration 2 to determine increments deltaa and deltab 0.0000 0.0000 0.0000 3.1412 -0.0003 0.0000 0.0000 -0.0010 0.0000 4.4146 5.6827 1.0103 -3.8953 0.4814 15.1734 -1.8751 0.4823 0.2326 4.4146 7.9117 -1.2482 -2.8344 -0.5935 8.0339 1.6822 -0.5928 0.3514 4.4146 13.4047 1.5798 -0.2201 0.7524 0.0484 -0.1656 0.7525 0.5662 2.8858 11.0271 -2.0706 1.1551 -1.2184 1.3342 -1.4074 -1.2187 1.4852 2.8858 13.6143 0.2046 2.6781 0.1209 7.1720 0.3238 0.1203 0.0145 2.8858 17.1555 0.5183 4.7626 0.3056 22.6824 1.4553 0.3044 0.0927 21.9012 68.7961 -0.0057 0.0000 0.0000 54.4443 0.0132 0.0002 2.7427 Iteration 3 to determine fi, gi, and hi 0.0000 0.0000 0.0000 0.0000 0.5788 2.1597 0.0000 0.0000 0.0000 0.0000 1.2000 0.2000 3.4000 0.2000 4.4138 1.2875 0.2296 2.1009 2.7050 0.4823 1.9000 0.2000 4.4000 0.2000 4.4138 1.7924 -0.2822 2.1009 3.7657 -0.5928 2.9000 0.2000 7.2000 0.2000 4.4138 3.0366 0.3582 2.1009 6.3795 0.7525 4.0000 0.2000 8.5000 0.4000 2.8855 3.8212 -0.7174 1.6987 6.4909 -1.2187 4.7000 0.2000 10.8000 0.4000 2.8855 4.7176 0.0708 1.6987 8.0137 0.1203 5.9000 0.2000 13.5000 0.4000 2.8855 5.9447 0.1792 1.6987 10.0980 0.3044 Iteration 3 to determine increments deltaa and deltab 0.0000 0.0000 0.0000 3.1414 -0.0000 0.0000 0.0000 -0.0000 0.0000 4.4138 5.6829 1.0133 -3.8947 0.4823 15.1685 -1.8785 0.4823 0.2327 4.4138 7.9113 -1.2454 -2.8340 -0.5928 8.0315 1.6800 -0.5928 0.3514 4.4138 13.4027 1.5809 -0.2202 0.7525 0.0485 -0.1657 0.7525 0.5662 2.8855 11.0258 -2.0702 1.1548 -1.2187 1.3335 -1.4073 -1.2187 1.4852 2.8855 13.6126 0.2043 2.6776 0.1203 7.1695 0.3220 0.1203 0.0145 2.8855 17.1531 0.5171 4.7619 0.3044 22.6756 1.4496 0.3044 0.0927 21.8977 68.7884 -0.0000 0.0000 0.0000 54.4271 0.0001 0.0000 2.7427
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.5788 Estimate of slope 2.1597
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.4764 Standard uncertainty associated with estimate of slope 0.1355 Covariance associated with estimates of intercept and slope -0.0577
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 2.743 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).