Contents

Function for validation and display of results

This function compares the results obtained from applications of the GUM uncertainty framework and a Monte Carlo method as implementations of the propagation of distributions. The values of the estimates of the output quantity, the associated standard uncertainties and the coverage intervals for the output quantity returned by the two approaches are displayed. The probability distributions provided by the two approaches to characterize the output quantity are displayed graphically. The results obtained using a Monte Carlo method are used to validate those obtained using the GUM uncertainty framework.

Section 6 of M G Cox, P M Harris and I M Smith, Software specifications for uncertainty evaluation, 2010, describes the procedure for validating the results obtained from the GUM uncertainty framework against those obtained from a Monte Carlo method.

function GUM_validation(yGUM, uyGUM, IyGUM, pdfGUM, ...
                        yMCM, uyMCM, IyMCM, pdfMCM, conv, Nbins, nval)

Input

yGUM          Estimate of the output quantity obtained from the
              GUM uncertainty framework
uyGUM         Standard uncertainty associated with the estimate yGUM
IyGUM  1 x 2  Lower and upper endpoints of a coverage interval obtained
              from the GUM uncertainty framework
pdfGUM 1 x 4  Probability density function for the output quantity
              obtained from the GUM uncertainty framework
yMCM          Estimate of the output quantity obtained from a
              Monte Carlo method
uyMCM         Standard uncertainty associated with the estimate yMCM
IyMCM  1 x 2  Lower and upper endpoints of a coverage interval obtained
              from a Monte Carlo method
pdfMCM 2 x Nb Approximation to the probability density function for the
              output quantity obtained from a Monte Carlo method
conv          Indicates whether the results of the Monte Carlo
              calculation have stabilized (1) or not (0)
Nbins         Number of bins used to display the approximation pdfMCM
nval          Number of significant digits for validating the results
              obtained from the GUM uncertainty framework

Step 1: Display the results obtained from the GUM uncertainty framework and a Monte Carlo method

Display the estimates, the standard uncertainties associated with the estimates, and the endpoints of the coverage intervals for the output quantity returned by the GUM uncertainty framework and a Monte Carlo method.

   disp('Results from GUM uncertainty framework:')
   fprintf(1, 'y    = %g\n', yGUM)
   fprintf(1, 'u(y) = %g\n', uyGUM)
   fprintf(1, 'I(y) = %g, %g\n\n', IyGUM(1), IyGUM(2))
   disp('Results from a Monte Carlo method:')
   fprintf(1, 'y    = %g\n', yMCM)
   fprintf(1, 'u(y) = %g\n', uyMCM)
   fprintf(1, 'I(y) = %g, %g\n\n', IyMCM(1), IyMCM(2))

Step 2: Validate the results obtained from the GUM uncertainty framework against those from a Monte Carlo method

The results obtained from the GUM uncertainty framework are only validated against those from a Monte Carlo method in the case that the results obtained from the Monte Carlo method have stabilized. The validation is based on comparing the absolute differences between the endpoints of the coverage intervals provided by the two approaches with a numerical tolerance. The numerical tolerance is set according to the number of significant decimal digits that are regarded as meaningful in the standard uncertainty associated with an estimate of the output quantity.

   if conv
      disp('MCM calculation has stabilized')
      if uyMCM == 0
         tol = 0;
      else
         r   = -(floor(log10(uyMCM)) - (nval - 1));
         tol = 0.5*10^(-r);
      end
      if abs(IyGUM(1) - IyMCM(1)) <= tol && abs(IyGUM(2) - IyMCM(2)) <= tol
         disp('GUF is validated against MCM')
      else
         disp('GUF is NOT validated against MCM')
      end
   else
      disp('MCM calculation has NOT stabilized')
   end

Step 3: Display the probability distributions obtained from the GUM uncertainty framework and a Monte Carlo method

The approximation to the probability distribution obtained from a Monte Carlo method is displayed as a scaled frequency distribution (histogram). That obtained from the GUM uncertainty framework is a Gaussian or t-distribution, and is displayed as a red continuous curve. (A probability distribution is not displayed in the case that a standard uncertainty is zero.)

   if uyMCM > 0
      axis_MCM = plot_pdfe(pdfMCM, Nbins);
      hold on
      if uyGUM > 0
         axis_GUM = plot_pdf(pdfGUM);
      end
   end

Set the axis limits so that both probability distributions are displayed.

   if uyMCM > 0
      axis_pdf = axis_MCM;
      if uyGUM > 0
         axis_pdf(1) = min([axis_pdf(1), axis_GUM(1)]);
         axis_pdf(2) = max([axis_pdf(2), axis_GUM(2)]);
         axis_pdf(3) = min([axis_pdf(3), axis_GUM(3)]);
         axis_pdf(4) = max([axis_pdf(4), axis_GUM(4)]);
      end
      title('Probability distribution for the output quantity')
      xlabel('Value')
      ylabel('Probability density')
      axis(axis_pdf)
   end

Step 4: Display the coverage intervals obtained from the GUM uncertainty framework and a Monte Carlo method

Add vertical lines to show the endpoints of the coverage intervals obtained from the two approaches. The endpoints obtained from a Monte Carlo method are displayed as blue continuous lines. Those obtained from the GUM uncertainty framework are displayed as red broken lines. (A coverage interval is not displayed in the case that a standard uncertainty is zero.)

   if uyMCM > 0
      plot([IyMCM(1) IyMCM(1)],  [axis_pdf(3) axis_pdf(4)], 'b-')
      plot([IyMCM(2) IyMCM(2)],  [axis_pdf(3) axis_pdf(4)], 'b-')
      if uyGUM > 0
         plot([IyGUM(1) IyGUM(1)],  [axis_pdf(3) axis_pdf(4)], 'r--')
         plot([IyGUM(2) IyGUM(2)],  [axis_pdf(3) axis_pdf(4)], 'r--')
      end
      hold off
   end

Function for plotting a Gaussian or t-distribution

This function displays a Gaussian or t-distribution. Such distributions are provided by the GUM uncertainty framework to characterize the output quantity of a measurement model.

function axlim = plot_pdf(pdf)

Input

pdf    1 x 4 Probability density function for the output quantity:
             [1 or 2, yGUM, uyGUM, veff ]

Output

axlim  1 x 4 Axis limits for the graph

Display a Gaussian distribution (assuming the standard deviation is positive).

   if pdf(1) == 1
      c = pdf(2);
      s = pdf(3);
      x = linspace(c-5*s, c+5*s, 1001)';
      y = (1/(s*sqrt(2*pi)))*exp(-(x - c).^2/(2*s^2));
      plot(x, y, 'r')
   end

Display a t-distribution (assuming the scale parameter is positive).

   if pdf(1) == 2
      c  = pdf(2);
      s  = pdf(3);
      nu = pdf(4);
      x = linspace(c-5*s, c+5*s, 1001)';
      k = gamma((nu+1)/2)/(s*sqrt(nu*pi)*gamma(nu/2));
      y = k*(1+((x-c)/s).^2/nu).^(-(nu+1)/2);
      plot(x, y, 'r')
   end

Set the axis limits for the graph.

   axlim = [min(x), max(x), -0.05*max(y), 1.05*max(y)];
   axis(axlim)

Function for plotting an approximation to the probability density function

This function displays an approximation to the probability distribution for a quantity such as may be provided by a Monte Carlo calculation.

function axlim = plot_pdfe(pdf, Nbins)

Input

pdf   2 x Nb Approximation to the probability density function for
             the output quantity (assuming Nb > 1)
Nbins        Number of bins to be used for displaying the approximation

Output

axlim 1 x 4  Axis limits for the plot

Start by adding bins at each end, with zero probability density, so that the total number of bins is an integer multiple of Nbins.

   gx = pdf(1,:);
   gy = pdf(2,:);
   bw = gx(2) - gx(1);
   nb = Nbins - rem(length(gx), Nbins);
   if nb > 0 && nb < Nbins
      nbl = round(nb/2);
      gxl = gx(1)-nbl*bw:bw:gx(1)-bw;
      gyl = zeros(1,nbl);
      nbr = nb - nbl;
      gxr = gx(end)+bw:bw:gx(end)+nbr*bw;
      gyr = zeros(1,nbr);
      gx  = [gxl, gx, gxr];
      gy  = [gyl, gy, gyr];
   end

Combine consecutive bins so that the number of bins is Nbins.

   R   = length(gx)/Nbins;
   gxp = mean(reshape(gx, R, Nbins));
   gyp = mean(reshape(gy, R, Nbins));

Use MATLAB's bar function to display the approximation to the probability density function.

   bar(gxp, gyp, 1)
   h = findobj(gca, 'Type', 'patch');
   set(h, 'FaceColor', 'w', 'EdgeColor', 'b')

Set the axis limits for the graph.

   axlen = (max(gxp) - min(gxp))/10;
   axlim = [min(gxp)-axlen, max(gxp)+axlen, -0.05*max(gyp), 1.05*max(gyp)];
   axis(axlim)

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