Contents
- Function for validation and display of results
- Step 1: Display the results obtained from the GUM uncertainty framework and a Monte Carlo method
- Step 2: Validate the results obtained from the GUM uncertainty framework against those from a Monte Carlo method
- Step 3: Display the probability distributions obtained from the GUM uncertainty framework and a Monte Carlo method
- Step 4: Display the coverage intervals obtained from the GUM uncertainty framework and a Monte Carlo method
- Function for plotting a Gaussian or t-distribution
- Function for plotting an approximation to the probability density function
- Program identifier
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