Software to support ISO/TS 28037:2010(E)
Contents
Introduction
This document runs the numerical example of weighted least squares (WLS) with known equal weights described in Clause 6 (Model for uncertainties associated with the y_i), and performs the prediction described in 11.1 EXAMPLE 1 and forward evaluation described in 11.2 EXAMPLE.
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]'; m = length(x);
Assign y-values.
y = [3.3 5.6 7.1 9.3 10.7 12.1]';
Assign uncertainties associated with y-values.
uy = [0.5 0.5 0.5 0.5 0.5 0.5]';
Obtain estimates of the straight line calibration function parameters and associated standard uncertainties and covariance
Solve the weighted least squares problem to obtain best fit straight-line parameters.
Step 1.
w = ones(m, 1)./uy; F2 = sum(w.*w);
Step 2.
g0 = (sum(w.*w.*x))/F2; h0 = (sum(w.*w.*y))/F2;
Step 3.
g = w.*(x - g0); h = w.*(y - h0);
Step 4.
G2 = sum(g.*g);
Step 5.
b = (sum(g.*h))/G2; a = h0 - b*g0;
Step 6.
u2a = 1/F2 + g0^2/G2; u2b = 1/G2; uab = -g0/G2;
Step 7.
r = w.*(y - a - b*x);
Step 8.
chi_sq_obs = sum(r.*r); nu = m - 2;
Step 9.
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 YI \n\n') fprintf('ISO/TS 28037:2010(E) CLAUSE 6 \n') fprintf('EXAMPLE (EQUAL WEIGHTS) \n\n')
MODEL FOR UNCERTAINTIES ASSOCIATED WITH THE YI ISO/TS 28037:2010(E) CLAUSE 6 EXAMPLE (EQUAL WEIGHTS)
Measurement data.
fprintf('FITTING \n') fprintf(['Data representing ', num2str(m),' measurement points, equal weights \n']) for i = 1:m fprintf('%8.1f %8.1f %8.1f \n', [x(i), y(i), uy(i)]) end fprintf('\n')
FITTING Data representing 6 measurement points, equal weights 1.0 3.3 0.5 2.0 5.6 0.5 3.0 7.1 0.5 4.0 9.3 0.5 5.0 10.7 0.5 6.0 12.1 0.5
Calculation tableau: see write_wls_tableau.m.
write_wls_tableau(x, y, w, g0, h0, g, h, a, b, r, '%8.3f ');
Calculation tableau 0.000 0.000 0.000 0.000 3.500 8.017 0.000 0.000 1.867 0.000 2.000 4.000 4.000 13.200 -5.000 -9.433 25.000 47.167 -0.648 0.419 2.000 4.000 8.000 22.400 -3.000 -4.833 9.000 14.500 0.438 0.192 2.000 4.000 12.000 28.400 -1.000 -1.833 1.000 1.833 -0.076 0.006 2.000 4.000 16.000 37.200 1.000 2.567 1.000 2.567 0.810 0.655 2.000 4.000 20.000 42.800 3.000 5.367 9.000 16.100 0.095 0.009 2.000 4.000 24.000 48.400 5.000 8.167 25.000 40.833 -0.619 0.383 0.000 24.000 84.000 192.400 0.000 0.000 70.000 123.000 1.757 1.665
Solution estimates.
fprintf('Estimate of intercept \n'), fprintf('%8.3f \n\n', a) fprintf('Estimate of slope \n'), fprintf('%8.3f \n\n', b)
Estimate of intercept 1.867 Estimate of slope 1.757
Standard uncertainties associated with solution estimates.
fprintf('Standard uncertainty associated with estimate of intercept \n'), fprintf('%8.3f \n\n', sqrt(u2a)) fprintf('Standard uncertainty associated with estimate of slope \n'), fprintf('%8.3f \n\n', sqrt(u2b)) fprintf('Covariance associated with estimates of intercept and slope \n'), fprintf('%8.3f \n\n', uab)
Standard uncertainty associated with estimate of intercept 0.465 Standard uncertainty associated with estimate of slope 0.120 Covariance associated with estimates of intercept and slope -0.050
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 1.665 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 errorbar(x, y, 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)


Prediction
yp = 10.5; uyp = 0.5; xp = (yp - a)/b; ca = -1/b; cb = -(yp - a)/(b^2); cy = 1/b; u2xp = ca^2*u2a + cb^2*u2b + 2*ca*cb*uab + cy^2*uyp^2; fprintf('ISO/TS 28037:2010(E) 11.1 \n') fprintf('EXAMPLE 1 \n\n') fprintf('PREDICTION \n') fprintf('Measured y-value \n'), fprintf('%8.3f \n\n', yp) fprintf('Standard uncertainty associated with measured y-value \n'), fprintf('%8.3f \n\n', uyp) fprintf('Estimate of x \n'), fprintf('%8.3f \n\n', xp) fprintf('Sensitivity coefficient wrt a \n'), fprintf('%8.3f \n\n', ca) fprintf('Sensitivity coefficient wrt b \n'), fprintf('%8.3f \n\n', cb) fprintf('Sensitivity coefficient wrt y \n'), fprintf('%8.3f \n\n', cy) fprintf('Standard uncertainty associated with estimate of x \n'), fprintf('%8.3f \n\n', sqrt(u2xp))
ISO/TS 28037:2010(E) 11.1 EXAMPLE 1 PREDICTION Measured y-value 10.500 Standard uncertainty associated with measured y-value 0.500 Estimate of x 4.913 Sensitivity coefficient wrt a -0.569 Sensitivity coefficient wrt b -2.796 Sensitivity coefficient wrt y 0.569 Standard uncertainty associated with estimate of x 0.322
Forward evaluation
xf = 3.5; uxf = 0.2; yf = a + b*xf; ca = 1; cb = xf; cx = b; u2yf = ca^2*u2a + cb^2*u2b + 2*ca*cb*uab + cx^2*uxf^2; fprintf('ISO/TS 28037:2010(E) 11.2 \n') fprintf('EXAMPLE \n\n') fprintf('FORWARD EVALUATION \n') fprintf('Measured x-value \n'), fprintf('%8.3f \n\n', xf) fprintf('Standard uncertainty associated with measured x-value \n'), fprintf('%8.3f \n\n', uxf) fprintf('Estimate of y \n'), fprintf('%8.3f \n\n', yf) fprintf('Sensitivity coefficient wrt a \n'), fprintf('%8.3f \n\n', ca) fprintf('Sensitivity coefficient wrt b \n'), fprintf('%8.3f \n\n', cb) fprintf('Sensitivity coefficient wrt x \n'), fprintf('%8.3f \n\n', cx) fprintf('Standard uncertainty associated with estimate of y \n'), fprintf('%8.3f \n\n', sqrt(u2yf))
ISO/TS 28037:2010(E) 11.2 EXAMPLE FORWARD EVALUATION Measured x-value 3.500 Standard uncertainty associated with measured x-value 0.200 Estimate of y 8.017 Sensitivity coefficient wrt a 1.000 Sensitivity coefficient wrt b 3.500 Sensitivity coefficient wrt x 1.757 Standard uncertainty associated with estimate of y 0.406
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).