Software to support ISO/TS 28037:2010(E)
Contents
Introduction
This document runs the numerical example of Gauss-Markov regression (GMR) described in Clause 9 (Model for uncertainties and covariances associated with 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.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0]'; m = length(x);
Assign y-values.
y = [1.3 4.1 6.9 7.5 10.2 12.0 14.5 17.1 19.5 21.0]';
Assign covariance matrix associated with y-values.
Uy = [2 1 1 1 1 0 0 0 0 0 1 2 1 1 1 0 0 0 0 0 1 1 2 1 1 0 0 0 0 0 1 1 1 2 1 0 0 0 0 0 1 1 1 1 2 0 0 0 0 0 0 0 0 0 0 5 4 4 4 4 0 0 0 0 0 4 5 4 4 4 0 0 0 0 0 4 4 5 4 4 0 0 0 0 0 4 4 4 5 4 0 0 0 0 0 4 4 4 4 5];
Obtain estimates of the straight line calibration function parameters and associated standard uncertainties and covariance
Solve the Gauss-Markov regression problem to obtain best fit straight-line parameters.
Step 1.
Ly = chol(Uy, 'lower');
Step 2.
e = ones(m, 1); f = Ly\e; g = Ly\x; h = Ly\y;
Step 3.
F2 = sum(f.*f);
Step 4.
g0 = sum(f.*g)/F2; h0 = sum(f.*h)/F2;
Step 5.
gt = g - g0*f; ht = h - h0*f;
Step 6.
Gt2 = sum(gt.*gt);
Step 7.
b = sum(gt.*ht)/Gt2; a = h0 - b*g0;
Step 8.
u2a = 1/F2 + g0^2/Gt2; u2b = 1/Gt2; uab = -g0/Gt2;
Step 9.
r = ht - b*gt;
Step 10.
chi_sq_obs = sum(r.*r); nu = m - 2;
Step 11.
chi_sq = calc_chi_sq_95_percent_quantile(nu);
Display information on screen and generate figures
Measurement model.
fprintf('\nMODEL FOR UNCERTAINTIES AND COVARIANCES ASSOCIATED WITH THE THE YI \n\n') fprintf('ISO/TS 28037:2010(E) CLAUSE 9 \n') fprintf('EXAMPLE \n\n')
MODEL FOR UNCERTAINTIES AND COVARIANCES ASSOCIATED WITH THE THE YI ISO/TS 28037:2010(E) CLAUSE 9 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 10 measurement points 1.0 1.3 2.0 4.1 3.0 6.9 4.0 7.5 5.0 10.2 6.0 12.0 7.0 14.5 8.0 17.1 9.0 19.5 10.0 21.0
Covariance matrix Uy.
fprintf('Covariance matrix Uy \n')
disp(Uy)
Covariance matrix Uy 2 1 1 1 1 0 0 0 0 0 1 2 1 1 1 0 0 0 0 0 1 1 2 1 1 0 0 0 0 0 1 1 1 2 1 0 0 0 0 0 1 1 1 1 2 0 0 0 0 0 0 0 0 0 0 5 4 4 4 4 0 0 0 0 0 4 5 4 4 4 0 0 0 0 0 4 4 5 4 4 0 0 0 0 0 4 4 4 5 4 0 0 0 0 0 4 4 4 4 5
Cholesky factor Ly of Uy.
fprintf('Cholesky factor Ly of Uy \n')
disp(Ly)
Cholesky factor Ly of Uy Columns 1 through 7 1.4142 0 0 0 0 0 0 0.7071 1.2247 0 0 0 0 0 0.7071 0.4082 1.1547 0 0 0 0 0.7071 0.4082 0.2887 1.1180 0 0 0 0.7071 0.4082 0.2887 0.2236 1.0954 0 0 0 0 0 0 0 2.2361 0 0 0 0 0 0 1.7889 1.3416 0 0 0 0 0 1.7889 0.5963 0 0 0 0 0 1.7889 0.5963 0 0 0 0 0 1.7889 0.5963 Columns 8 through 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.2019 0 0 0.3698 1.1435 0 0.3698 0.2691 1.1114
Calculation tableau: see write_gmr_tableaux.m.
write_gmr_tableaux(f, g, h, g0, h0, gt, ht, a, b, r, '%9.4f ');
Initial calculation tableau 0.7071 0.7071 0.9192 0.4082 1.2247 2.8169 0.2887 1.7321 4.4167 0.2236 2.2361 3.9578 0.1826 2.7386 5.6963 0.4472 2.6833 5.3666 0.1491 1.6398 3.6522 0.0925 1.8490 4.4284 0.0673 2.2198 5.3208 0.0529 2.6463 5.5360 Main calculation tableau 0.0000 0.0000 0.0000 4.1111 8.4044 0.0000 0.0000 -0.6456 0.0000 0.5000 0.5000 0.6500 -2.1999 -5.0236 4.8395 11.0514 -0.1809 0.0327 0.1667 0.5000 1.1500 -0.4536 -0.6142 0.2058 0.2786 0.3844 0.1477 0.0833 0.5000 1.2750 0.5453 1.9906 0.2973 1.0854 0.7902 0.6245 0.0500 0.5000 0.8850 1.3168 2.0785 1.7340 2.7370 -0.8202 0.6727 0.0333 0.5000 1.0400 1.9880 4.1619 3.9523 8.2739 -0.2145 0.0460 0.2000 1.2000 2.4000 0.8447 1.6080 0.7136 1.3583 -0.2516 0.0633 0.0222 0.2444 0.5444 1.0269 2.3994 1.0546 2.4640 0.1387 0.0192 0.0085 0.1709 0.4094 1.4689 3.6514 2.1578 5.3636 0.4177 0.1745 0.0045 0.1493 0.3579 1.9433 4.7555 3.7763 9.2412 0.4777 0.2282 0.0028 0.1401 0.2930 2.4287 5.0912 5.8986 12.3650 -0.2552 0.0651 1.0714 4.4048 9.0048 0.0000 0.0000 24.6296 54.2185 2.2014 2.0740
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.6456 Estimate of slope 2.2014
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 1.2726 Standard uncertainty associated with estimate of slope 0.2015 Covariance associated with estimates of intercept and slope -0.1669
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 8 Observed chi-squared value 2.074 95 % quantile of chi-squared distribution with 8 degrees of freedom 15.507 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 errorbar(x, y, sqrt(diag(Uy)), 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 6) plot(x, a + b*x, '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(i)], 'k-') plot(x(i), r(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).