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) and Annex C EXAMPLE 1 using the orthogonal factorization approach described in C.2.
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 factor Bx of covariance matrix Ux associated with x-values.
Bx = [0.5 0 0 0 0 0 0 0.5 0 0 0 0.5 0 0 0 0 0 0 1 0 0 0 0.5 0 0 0 0 0.5 1 0 0 0 0 0.5 0 0 0 0 0 1 0 0 0 0 0.5 0 0 0.5 0 1 0 0 0 0 0 0.5 0 0 1 1 0 0 0 0 0 0 0.5 0.5 1 1];
Assign factor By of covariance matrix Uy associated with y-values.
By = [2 0 0 0 0 0 0 1 0 2 0 0 0 0 0 1 0 0 2 0 0 0 0 1 0 0 0 2 0 0 0 1 0 0 0 0 2 0 0 1 0 0 0 0 0 2 0 1 0 0 0 0 0 0 2 1];
Build complete factor for x- and y-values (assuming the correlation associated with each (x, y)-pair is zero).
B = blkdiag(Bx, By);
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(By*By')));
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} = []; Q{1} = []; R10{1} = []; R1{1} = []; T{1} = []; Z{1} = []; ft{1} = []; ft1{1} = []; ft2{1} = []; T11{1} = []; T12{1} = []; T22{1} = []; et2{1} = [];
Steps 2 to 8: see algm_ggmr_orthogonal_steps_2_to_8.m.
ind = 1; [tt, dt, f, J, Q, R1, T, Z, ft1, ft2, T11, T12, T22, et2] = ... algm_ggmr_orthogonal_steps_2_to_8(x, y, B, tt, dt, f, J, Q, R1, T, Z, ... ft1, ft2, T11, T12, T22, et2, ind); while any(abs(dt{ind}) > tol)
Update iteration number.
ind = ind + 1;
Step 9. Repeat steps 2 to 8 until convergence has been achieved: see algm_ggmr_orthogonal_steps_2_to_8.m.
[tt, dt, f, J, Q, R1, T, Z, ft1, ft2, T11, T12, T22, et2] = ... algm_ggmr_orthogonal_steps_2_to_8(x, y, B, tt, dt, f, J, Q, R1, T, Z, ... ft1, ft2, T11, T12, T22, et2, ind);
end
a = tt{ind+1}(m+1);
b = tt{ind+1}(m+2);
Step 10. Solve upper-triangular system and evaluate uncertainties.
Ra = R1{ind}(m+1:m+2, m+1:m+2); Ta = T11{ind}(end-1:end, end-1:end); Ka = Ra\Ta; Ua = Ka*Ka'; u2a = Ua(1, 1); u2b = Ua(2, 2); uab = Ua(1, 2);
Step 11. Form observed chi-squared value and degrees of freedom.
chi_sq_obs = sum(et2{ind}.*et2{ind}); nu = m - 2;
Step 12. 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) ANNEX C \n') fprintf('EXAMPLE 1 \n\n')
MODEL FOR UNCERTAINTIES AND COVARIANCES ASSOCIATED WITH THE XI AND THE YI ISO/TS 28037:2010(E) ANNEX C EXAMPLE 1
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(Bx*Bx')
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
Factor Bx of Ux.
fprintf('Factor Bx of Ux \n')
disp(Bx)
Factor Bx of Ux Columns 1 through 7 0.5000 0 0 0 0 0 0 0 0.5000 0 0 0 0 0 0 0 0.5000 0 0 0 0 0 0 0 0.5000 0 0 0 0 0 0 0 0.5000 0 0 0 0 0 0 0 0.5000 0 0 0 0 0 0 0 0.5000 Columns 8 through 10 0.5000 0 0 0 1.0000 0 0.5000 1.0000 0 0 0 1.0000 0.5000 0 1.0000 0 1.0000 1.0000 0.5000 1.0000 1.0000
Covariance matrix Uy.
fprintf('Covariance matrix Uy \n')
disp(By*By')
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
Factor By of Uy.
fprintf('Factor By of Uy \n')
disp(By)
Factor By of Uy 2 0 0 0 0 0 0 1 0 2 0 0 0 0 0 1 0 0 2 0 0 0 0 1 0 0 0 2 0 0 0 1 0 0 0 0 2 0 0 1 0 0 0 0 0 2 0 1 0 0 0 0 0 0 2 1
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.3018e-08 3.4240e-01 1.0011e+00 1.0962e-04 2.1135e-05 7.6429e-09 3.3717e-10 1.0012e+00
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).