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