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