Contents

Function implementing the GUM uncertainty framework

This function implements the GUM uncertainty framework for uncertainty evaluation for a measurement model

$$Y = f(X_1, X_2, ..., X_N),$$

with a scalar output quantity and a general number of uncorrelated input quantities.

The inputs to the function are a specification model of the measurement model, a specification pdfin of the probability distributions assigned to the input quantities in the measurement model, and a coverage probability p. The function returns an estimate yGUM of the output quantity, the standard uncertainty uyGUM associated with the estimate, the endpoints IyGUM of the probabilistically symmetric coverage interval for the output quantity corresponding to coverage probability p, and a specification pdfGUM of the probability distribution used to characterize the output quantity.

The estimate is obtained by evaluating the model at the expectations of the input quantities and the associated standard uncertainty by applying the law of propagation of uncertainty. The output quantity is characterized by a Gaussian distribution or a scaled and shifted t-distribution, and the distribution is used as the basis of evaluating the coverage interval for the output quantity.

model is a string containing the name of a function for evaluating the measurement model for values of the input quantities and calculating sensitivity coefficients. For an example of such a function, see model_additive.

Each input quantity in the measurement model is assigned a probability distribution, described by a probability density function, which may be (1) a Gaussian (or normal) distribution, (2) a scaled and shifted t-distribution, (3) a rectangular (or uniform) distribution, (4) a curvilinear trapezoidal distribution, or (5) a U-shaped distribution. (This list includes all the univariate probability distributions encountered in the examples of GUM Supplement 1. There is no consideration of joint probability distributions, such as the multivariate Gaussian distribution.)

pdfin is a matrix of N rows and 4 columns, the ith row of which defines the probability distribution for the ith input quantity. For a measurement model with a single input quantity, pdfin can take the following forms:

pdfin = [1, mu, sigma, inf] defines a Gaussian distribution with expectation mu and standard deviation sigma (the fourth element of pdfin is not used);

pdfin = [2, mu, sigma, nu] defines a t-distribution with shift parameter mu, scale parameter sigma and degrees of freedom nu;

pdfin = [3, a, b, inf] defines a rectangular distribution with lower limit a and upper limit b (the fourth element of pdfin is not used);

pdfin = [4, a, b, r] defines a curvilinear trapezoidal distribution with lower limit a, upper limit b and (fractional) reliability r for the semi-width;

pdfin = [5, a, b, inf] defines a U-shaped distribution with lower limit a and upper limit b (the fourth element of pdfin is not used).

Guidance on the assignment, in some common circumstances, of probability density functions for the input quantities in a measurement model is given in GUM Supplement 1 [3, section 6].

The coverage probability p takes one of the values 0.90, 0.95 or 0.99.

Section 4.1 of M G Cox, P M Harris and I M Smith, Software specifications for uncertainty evaluation, 2010, describes the procedure implemented in the GUM uncertainty framework.

Section 3.1.1 of M G Cox, P M Harris and I M Smith, Software for GUM Supplement 1: User Manual, 2011, describes the specifications of the different probability distributions. (Note that this document describes other probability distributions that are not considered here.)

function [yGUM, uyGUM, IyGUM, pdfGUM] = GUM_calculation(model, pdfin, p)

Input

model  string Function defining the measurement model
pdfin  N x 4  Probability density functions for the input quantities in
              the measurement model
p             Coverage probability (0.90, 0.95 or 0.99)

Output

yGUM          Estimate of the output quantity
uyGUM         Standard uncertainty associated with the estimate of the
              output quantity
IyGUM  1 x 2  Lower and upper endpoints of a coverage interval for the
              output quantity
pdfGUM 1 x 4  Probability density function for the output quantity:
              [1 or 2, yGUM, uyGUM, veff]

Step 1: Summarize the probability distributions for the input quantities

Each input quantity is summarized by an estimate (expectation) x, associated standard uncertainty (standard deviation) ux, and associated degrees of freedom nu, obtained from the probability density function assigned to the quantity.

  N = size(pdfin,1);
  for i = 1:N
     [x(i,1), ux(i,1), nu(i,1)] = convert_pdf(pdfin(i,:));
  end

Step 2: Apply the law of propagation of uncertainty

The model is evaluated for the estimates of the input quantities to obtain an estimate yGUM of the output quantity. The first order partial derivatives of the model with respect to the input quantities are evaluated for the estimates of the input quantities to obtain the sensitivity coefficients c. The law of propagation of uncertainty is applied to obtain the standard uncertainty uyGUM associated with yGUM.

See sections 4.2 and 4.3 of M G Cox, P M Harris and I M Smith, Software specifications for uncertainty evaluation, 2010.

  [yGUM, c] = eval([model, '(x)']);
  uyGUM     = sqrt(c*diag(ux.^2)*c');

Step 3: Characterize the output quantity by a probability distribution

The output quantity is characterized by a Gaussian distribution (for an effective degrees of freedom exceeding 100) or a scaled and shifted t-distribution otherwise. The Welch-Satterthwaite formula is used to evaluate the effective degrees of freedom veff.

See section 4.4 of M G Cox, P M Harris and I M Smith, Software specifications for uncertainty evaluation, 2010.

  denom_sum = sum(((c'.*ux).^4)./nu);
  if denom_sum == 0
     veff   = inf;
     pdfGUM = [1, yGUM, uyGUM, veff];
  else
     veff   = fix((uyGUM^4)/denom_sum);
     if veff >  100
        veff   = inf;
        pdfGUM = [1, yGUM, uyGUM, veff];
     else
        pdfGUM = [2, yGUM, uyGUM, veff];
     end
  end

Step 4: Evaluate the coverage factor

The coverage factor kp is obtained as a percentage point of the standard Gaussian distribution (for an effective degrees of freedom exceeding 100) or a t-distribution otherwise.

See section 4.5 of M G Cox, P M Harris and I M Smith, Software specifications for uncertainty evaluation, 2010.

  kp = coverage_factor(p, veff);

Step 5: Evaluate the expanded uncertainty and coverage interval

The expanded uncertainty Uy is the product of the coverage factor kp and the standard uncertainty uyGUM. The probabilistically symmetric coverage interval (also the shortest coverage interval) is an interval centred on the estimate yGUM with half-width equal to the expanded uncertainty.

See section 4.5 of M G Cox, P M Harris and I M Smith, Software specifications for uncertainty evaluation, 2010.

  Uy    = kp*uyGUM;
  IyGUM = [(yGUM - Uy) (yGUM + Uy)];

Function for summarizing the probability distribution for an input quantity

This function summarizes an input quantity in a measurement model by an estimate (expectation), associated standard uncertainty (standard deviation), and associated degrees of freedom, obtained from the probability distribution assigned to the quantity.

function [x, ux, nu] = convert_pdf(pdf)

Input

pdf   1 x 4 Probability density function for an input quantity in a
            measurement model

Output

x           Estimate of the input quantity
ux          Standard uncertainty associated with the estimate of the
            input quantity
nu          Degrees of freedom associated with the standard uncertainty

The calculation depends on the nature of the probability distribution assigned to the input quantity.

   switch pdf(1)

Gaussian distribution.

      case {1}
         x  = pdf(2);
         ux = pdf(3);
         nu = inf;
%%%
% t-distribution.
      case {2}
         x  = pdf(2);
         ux = pdf(3);
         nu = pdf(4);
%%%
% Rectangular distribution.
      case {3}
         a  = pdf(2);
         b  = pdf(3);
         x  = (b + a)/2;
         ux = (b - a)/(2*sqrt(3));
         nu = inf;
%%%
% Curvilinear trapezoidal distribution.
      case {4}
         a  = pdf(2);
         b  = pdf(3);
         r  = pdf(4);
         x  = (b + a)/2;
         ux = (b - a)/(2*sqrt(3));
         nu = (1/2)*(1/r)^2;
%%%
% U-shaped distribution.
      case {5}
         a  = pdf(2);
         b  = pdf(3);
         x  = (b + a)/2;
         ux = (b - a)/(2*sqrt(2));
         nu = inf;
   end

Function for evaluating a coverage factor

In the GUM uncertainty framework, the random variable

$$Z = (Y - y)/u(y)$$

is characterized by the standard Gaussian distribution N(0, 1) in the case that the effective degrees of freedom veff associated with u(y) is infinite, or by a t-distribution with veff degrees of freedom otherwise. This function evaluates the coverage factor kp corresponding to coverage probability p such that

$${\rm Prob}(|Z| < kp) = p.$$

The coverage factor is evaluated as a percentage point of the distribution (standard Gaussian or t-distribution) used to characterize Z.

See section 4.5 of M G Cox, P M Harris and I M Smith, Software specifications for uncertainty evaluation, 2010.

function kp = coverage_factor(p, veff)

Input

p     Coverage probability (0.90, 0.95 or 0.99)
veff  Effective degrees of freedom (inf for a Gaussian distribution,
      finite for a t-distribution)

Output

kp    Coverage factor

Distinguish between a degrees of freedom that is infinite or finite.

   if isinf(veff)
      if p == 0.90
         kp = 1.645;
      elseif p == 0.95
         kp = 1.960;
      elseif p == 0.99
         kp = 2.576;
      end
   else
      kp = t_dist_coverage(p, veff);
   end

Function for evaluating a coverage factor for the t-distribution

This function evaluates the coverage factor kp corresponding to coverage probability p such that

$${\rm Prob}(|T| < kp) = p,$$

where T is characterized by a t-distribution with nu degrees of freedom.

See Table G.2 of BIPM, IEC, IFCC, ISO, IUPAC, IUPAP and OIML, Guide to the expression of uncertainty in measurement (ISBN 92-67-10188-9, corrected and reprinted, 1995).

function kp = t_dist_coverage(p, nu)

Input

p     Coverage probability (0.90, 0.95 or 0.99)
nu    Degrees of freedom

Output

kp    Coverage factor for the t-distribution

Look-up table TABLE for the coverage factor.

   TABLE = [1   6.31  12.71  63.66
            2   2.92   4.30   9.92
            3   2.35   3.18   5.84
            4   2.13   2.78   4.60
            5   2.02   2.57   4.03
            6   1.94   2.45   3.71
            7   1.89   2.36   3.50
            8   1.86   2.31   3.36
            9   1.83   2.26   3.25
            10  1.81   2.23   3.17
            11  1.80   2.20   3.11
            12  1.78   2.18   3.05
            13  1.77   2.16   3.01
            14  1.76   2.14   2.98
            15  1.75   2.13   2.95
            16  1.75   2.12   2.92
            17  1.74   2.11   2.90
            18  1.73   2.10   2.88
            19  1.73   2.09   2.86
            20  1.72   2.09   2.85
            25  1.71   2.06   2.79
            30  1.70   2.04   2.75
            35  1.70   2.03   2.72
            40  1.68   2.02   2.70
            45  1.68   2.01   2.69
            50  1.68   2.01   2.68
            100 1.660  1.984  2.626];

Identify which column of TABLE to use (column 2 for 90 %, column 3 for 95 %, column 4 for 99 %).

   if p == 0.90
      col = 2;
   elseif p == 0.95
      col = 3;
   elseif p == 0.99
      col = 4;
   end

Use MATLAB's interp1 function to linearly interpolate the values in the identified column of TABLE. (For degrees of freedom greater than 100, take nu = 100.)

   if nu > 100
      kp = TABLE(27, col);
   else
      kp = interp1(TABLE(:, 1), TABLE(:, col), nu);
   end

Program identifier

NPLUnc_101 Open-Source Software (Release 1.1): Software for measurement uncertainty evaluation

M G Cox, P M Harris and I M Smith, Mathematics and Scientific Computing Group, National Physical Laboratory, UK

Crown copyright 2011


Published with MATLAB® 7.11