- Email: [email protected]
Components of measurement uncertainty from a measurement model with two stages involving two output quantities
Component of measurement uncertainty from a measurement model with two stages involving two output quantities Dimitrios Theodorou, Ypatia Zannikou, Fanourios Zannikos PII: DOI: Reference:
S0169-7439(15)00145-8 doi: 10.1016/j.chemolab.2015.05.025 CHEMOM 3028
To appear in:
Chemometrics and Intelligent Laboratory Systems
Received date: Revised date: Accepted date:
10 February 2015 13 May 2015 22 May 2015
Please cite this article as: Dimitrios Theodorou, Ypatia Zannikou, Fanourios Zannikos, Component of measurement uncertainty from a measurement model with two stages involving two output quantities, Chemometrics and Intelligent Laboratory Systems (2015), doi: 10.1016/j.chemolab.2015.05.025
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT TITLE:
Component of measurement uncertainty from a measurement model
AUTHORS:
SC R
Dimitrios Theodorou, Ypatia Zannikou, Fanourios Zannikos
IP
T
with two stages involving two output quantities
INSTITUTION:
National Technical University of Athens, School of Chemical Engineering, Laboratory of Fuels
NU
and Lubricants Technology, Iroon Polytechneiou 9, Athens 15780, Greece
MA
Author Contact Details Dimitrios Theodorou
Tel./fax: +30 21 0652 5976.
(CORRESPONDING AUTHOR)
E-mail address: [email protected]
AC
CE P
Ypatia Zannikou
TE
D
Permanent Address: 59 A.Tritsi Street, 15238 Halandri, Athens, Greece.
Fanourios Zannikos
Tel 210 7723143 / Fax 210 7723163 E-mail address: [email protected]
Tel 210 7723143 / Fax 210 7723163 E-mail address: [email protected]
Abstract The GUM uncertainty framework and the Monte Carlo Method can be extended to include multistage measurement models and measurement models having more than one output quantity. A typical multistage measurement model, that also includes a measurement stage
ACCEPTED MANUSCRIPT with more than one output quantity is the calculation of the coefficients of a calibration curve, constructed by means of least squares regression and the subsequent use this calibration curve for the estimation of a quantity value. In the present work the Monte Carlo
T
Method was used to evaluate the component of measurement uncertainty from a
IP
calibration curve used for the determination of the total nitrogen found in a gasoline sample. The slope and the intercept of the calibration curve were treated as a vector output
SC R
quantity characterized by a joint probability distribution (bivariate Gaussian) and two types of coverage regions were estimated (rectangular and ellipsoid). The Monte Carlo Method algorithm employing fixed number of trials (106) gave a predicted value of 2.09 mg L-1, a
NU
standard measurement uncertainty of 0.04 mg L-1 and a 95% symmetrical coverage interval [2.01 – 2.17] mg L-1. The standard measurement uncertainty obtained by MCM agrees well
MA
with the outcome of the equation giving the standard error of the estimate.
TE
D
1. Introduction
Until recently there has been a focus on the evaluation of measurement uncertainty of
CE P
simple measurement models consisting of one stage and having a single output quantity (univariate scalar models). There are measurement models, however, where: (i) the output quantities from previous stages become input quantities for subsequent stages (multistage measurement models) and, (ii) there is more than one output quantity, depending on a
AC
common set of input quantities. The principles for the measurement uncertainty evaluation described in the “Guide to Expression of Uncertainty in Measurement” [1] (GUM uncertainty framework) and its Supplement [2] describing the application of Monte Carlo Method (MCM) can be extended to the evaluation of uncertainties associated with vector estimates of multivariate output quantities. In 2011, JCGM published the GUM supplement 2 [3] which provides guidance concerning the application of a generalized GUM uncertainty framework and MCM in measurement models with any number of output quantities. A typical multistage measurement model, that also includes a measurement stage with more than one output quantity is the calculation of the coefficients of a calibration curve, constructed by means of least squares regression and the subsequent use of this calibration curve for the estimation of a quantity value. In analytical chemistry, linear regression is very often used in the construction of calibration curves required for techniques such as liquid or
ACCEPTED MANUSCRIPT gas chromatography and atomic absorption or fluorescence spectrometry [4-6]. A calibration curve is actually the expression of the relation between indication and corresponding measured quantity value. Calibration establishes a relationship between the values of a
T
standard (reference values), xi and the output quantities (response of the instrument), yi.
IP
Using values of a standard (reference values), xi and instrument responses, yi, a relationship (often assumed to be represented by a straight line) is established and thenthe calibration
SC R
model is used in reverse, that is, to predict a value, xpred from an instrument response, y0 [79]. Fig. 1 illustrates the two-stage measurement model involving construction and use of a calibration function. As with most statistics, calibration curve coefficients (slope and
NU
intercept for linear calibration model) are only estimates based on a finite number of measurements, and therefore, their values are associated with uncertainties. This leads to an measurement uncertainty of the predicted value as well. In many cases the component of
MA
measurement uncertainty from the calibration curve contributes significantly to the
[FIGURE 1]
CE P
TE
D
measurement uncertainty [10] and therefore it has to be reliably evaluated.
The determination of the total nitrogen found in liquid hydrocarbons is conducted using oxidative combustion and chemiluminescence detection method [11], which requires the construction of a linear calibration curve. This determination is important and has to be
AC
reliable because even trace amounts of nitrogenous materials contained in the feedstocks may lead to poisoning of process catalysts used in petroleum and chemical refining. In the present work, the measurement uncertainty of the nitrogen content in liquid fuel products determined using a calibration curve is evaluated using MCM. The coefficients (slope and intercept) of the calibration curve are treated as a vector output quantity characterized by a joint probability distribution [12]. Bivariate coverage regions are specified taking the form an ellipsoid or a rectangle.
ACCEPTED MANUSCRIPT 2. Experimental
An Analytik Jena AG - multi EA® 3100 Element Analyzer equipped with automatic sampler
T
was employed in this work. This analyzer fully complies with ASTM D 4629-12 “Standard Test
IP
Method for Trace Nitrogen in Liquid Petroleum Hydrocarbons by Syringe/Inlet Oxidative
SC R
Combustion and Chemiluminescence Detection” [11]. Table 1 presents the operating conditions of the instrument. Six calibration standards (Analytik Jena nitrogen calibration kit) with concentrations 0.0, 0.1, 0.5, 1.0, 5.0 and 10.0 mg L-1 were analyzed in triplicate. The instrument responses (AU – Area Units) were recorded and a calibration curve with 18
NU
points was constructed. Then a gasoline sample was measured and its total nitrogen concentration in mg L-1 was predicted using the calibration curve constructed. The sample
MA
and the calibration standards were injected using a 40 μL syringe. Data produced are
CE P
TE
D
presented in Table 2.
[TABLE 1]
[TABLE 2]
AC
3. Linear regression by least-squares method
The least-squares method is used for data evaluation problems related to measurement models with an underlying theoretical relationship between an independent variable and a dependent variable. This relationship constitutes the basis of a curve-fitting problem. The task is to estimate the parameters of a calibration function from pairs of measured quantity values and their corresponding indication values [13].
The measured quantity values of an independent variable are typically those of measurement standards. The values of the dependent variables are the indication values returned by the measuring system for the corresponding values of the independent variables. In chemistry, the establishment of a calibration function has the purpose of obtaining a function that allows one to calculate the concentration of the analyte as a
ACCEPTED MANUSCRIPT function of an instrumental signal. A set of standards are prepared containing a known amount of the analyte of interest and the instrument response of each standard is measured. Then a relationship is established between the instrument response and the
T
analyte concentration. The most frequent type of calibration in chemical analyses is the
IP
indirect calibration, where the equipment gives a value (signal or instrument response),
SC R
which has a different quantity from that of the standard (reference value) [7].
Calibration data (pairs of analyte concentrations, xi, and instrument responses, yi ) are used as an input to the least square regression analysis which fits a predetermined measurement
NU
model to these data. The simplest measurement model is the one described by the linear
MA
function:
(1)
D
Y b0 b1 x
TE
where Y is the instrument response (dependent variable), x is the analyte concentration (independent variable) and b0 and b1 are the coefficients of the model known as the
CE P
intercept and slope, respectively. The slope and the intercept are estimated from the calibration data (xi, yi).
It has to be noted that least square regression process can validly be applied under the following conditions [7,8]:
AC
a. The linear model holds for the data (i.e. the response does indeed vary linearly with concentration). b. The errors in the independent variable are insignificant compared with those of the dependent variable. c. These errors are normally distributed and independent of the values of the dependent variable (i.e. the data are homoscedastic). Only when these conditions hold, the application of linear least square regression is expected to realize the best fit to the calibration data.
The quantities b0 and b1 are correlated and the degree of their correlation is characterized by the estimated correlation coefficient r(b0, b1):
ACCEPTED MANUSCRIPT n
x
r (b0 , b1 )
i 1
i
(2)
n
n xi2
T
i 1
IP
where n is the number calibration data (pairs of xi, yi).
SC R
Once the slope and the intercept have been calculated, it becomes possible to calculate the standard error of regression (residual standard deviation) which is required for the calculation of a number of regression statistics. The standard error of regression is a
NU
measure of the deviation of the data points from the regression line and is calculated using
MA
the equation:
yi yˆ i 2
SEregression
(3)
TE
D
n2
where yˆ (predicted y value) is the value of y calculated for a given x value using the linear
CE P
calibration equation.
The standard deviation for the estimate of the intercept, SE(b0), is calculated using the
AC
equation:
SE (b0 ) SEregression
n
x i 1
2 i
(4)
n
n ( xi x )
2
i 1
where x is the average concentration of all standards used. The standard deviation for the estimate of the slope, SE(b1), is calculated using the equation:
SE (b1 )
SEregression n
(x i 1
i
(5)
x)2
Once the least-squares method has been used to estimate the parameters of a calibration function and evaluate the associated standard uncertainties and covariances, the measuring
ACCEPTED MANUSCRIPT system is subsequently used for a measurement. When analyzing unknown test samples the relationship of Eq. (1) is used in reverse in order to predict a value from an instrument
IP
y 0 b0 b1
(6)
SC R
x pred
T
response:
where y0 is the instrument response of the unknown test sample and xpred is the prediction (or estimation) of the measurement model concerning the true value of the measurand (e.g. analyte concentration).
NU
The standard measurement uncertainty associated with the predicted value estimate is evaluated using the standard uncertainties and covariances associated with the calibration
MA
function parameter estimates and the standard measurement uncertainty associated with the indication value.
TE
D
4. Measurement uncertainty evaluation
The result obtained from a calibration curve is actually an estimate of the true value of the
CE P
measurand. Therefore, it is associated with a measurement uncertainty due to lack of fit of the linear model to the data, which leads to uncertainties (standard errors) of the parameters b0 and b1.
AC
The GUM and its supplements [1-3] provide guidance on measurement uncertainty evaluation using as information the best estimates and the standard uncertainties of the input quantities (b0, b1, y0). One approach is to propagate this information, using the law of propagation of uncertainty through a first-order Taylor series approximation to the measurement model to provide the estimate of the output quantity xpred and the standard measurement uncertainty u(xpred) associated with xpred (GUM uncertainty framework) [13]. It has to be noted though that, the GUM uncertainty framework might not be satisfactory, in situations, where: (i)
the measurement function is non-linear,
(ii)
the probability distributions for the input quantities are asymmetric,
(iii)
the measurement uncertainty contributions are not of approximately the same magnitude,
ACCEPTED MANUSCRIPT (iv)
the probability distribution for the output quantity is either asymmetric, or not a Gaussian or a t-distribution.
T
An alternative method for evaluating measurement uncertainty is the Monte Carlo method
IP
(MCM). This method is described in the first supplement of GUM [2] and involves no restrictions for valid application [14,15] except for the fact that it requires prior knowledge
SC R
concerning the type of the probability distribution assigned to each input parameter. The MCM actually combines and propagates distributions rather than propagating uncertainties as in the GUM uncertainty framework. The measurement uncertainty evaluation is based on
NU
a probabilistic approach that combines the whole distribution of the input parameters and is not just based on their best estimates and standard measurement uncertainties.
MA
The MCM consists of simulating draws from the distribution of the output quantity (xpred), based on simulated draws from the distributions of the input quantities (b0, b1, y0). In our case, the quantities b0 and b1 (output quantities for the 1st stage of the measurement model
D
and input quantities for the second stage of the measurement model) are not independent
TE
(non zero covariance) and therefore they should be described by a multivariate (joint) distribution. This distribution is used for the draws during the Monte Carlo Method
CE P
implementation. The number of trials selected for MCM determines the degree of numerical accuracy of the MCM outputs [16]. An adaptive procedure by which the number of trials can be chosen, may be used [17,18]. Nevertheless, the Monte Carlo numerical simulation tends
AC
to require up to 106 trials for calculating a standard measurement uncertainty or a coverage interval which is correct to one or two significant decimals digits.
The estimate of the output quantity xpred is estimated by the average of the M MCM trials which produce M measurement model values (x pred (k) , k = 1,…,M):
x pred
1 M
M
x k 1
(k ) pred
(7)
while the standard measurement uncertainty u(xpred) associated with xpred is evaluated as the standard deviation of the M model values:
ACCEPTED MANUSCRIPT 1 M (k ) ( x pred x pred ) M 1 k 1
(8)
T
u ( xpred )
SC R
IP
5. Results and discussion
5.1 Calibration results
NU
A least square linear regression was carried out for the calibration data (see Table 2) using the LINEST function of Microsoft Excel. The results of the regression analysis are presented
MA
in Table 3. The coefficient of determination for the straight line, r2, is 0.9994, and the regression ANOVA gave a Fisher – Snedecor value F equal to 26821 with a significance (p – value) of 3.1 10-27. This implies a strong linear relationship between concentration and
D
instrument response values. Fig. 2 shows the plot produced by the calibration data and the
[TABLE 3]
[FIGURE 2]
AC
CE P
TE
fitted line.
5.2 Measurement uncertainty of calibration function parameters
The estimation of the slope and the intercept of the calibration curve is based on a measurement model of multiple outputs (two in our case), which involves correlated data. It is more convenient to treat the coeffiecients of the calibration function as a vector quantity E (Eq. 11), associated with a covariance matrix V (Eq. 12)
b E 0 b1
(9)
ACCEPTED MANUSCRIPT u 2 (b0 ) u (b0 , b1 ) SE 2 (b0 ) r (b0 , b1 ) SE (b0 ) SE (b1 ) V 2 SE 2 (b1 ) u (b0 , b1 ) u (b1 ) r (b0 , b1 ) SE (b0 ) SE (b1 )
(10)
where u (b0 , b1 ) = u (b1 , b0 ) is the estimated covariance between b0 and b1. The standard
T
uncertainties u(b0) and u(b1) are provided by the regression statistics as standard error of
SC R
IP
intercept, SE(b0), and slope, SE(b1), respectively.
Once the probability density function (PDF) is established (bivariate Gaussian) a coverage region can be determined, for a coverage probability, p. A coverage region specifies a region
NU
in 2-dimensional space that contains E with probability p. Two types of coverage region can be considered [3]:
rectangle centered coverage region (separately determined coverage intervals for
MA
b1 and b0).
ellipse centered coverage region
D
TE
The rectangle centered coverage region (sides parallel to the axes) is determined by the marginal PDFs of b1 and b0 and the coverage intervals are: b1 ± kq u(b1) and b0 ± kq u(b0). For
CE P
coverage probability p=0.95, a coverage factor kq = 2.24, is used.
AC
The ellipse centered coverage region is determined as: (η – E)T V-1 (η – E) = kp2
or
n0 b0
n1 b1
(11)
1
u 2 (b0 ) u (b0 , b1 ) n0 b0 k p2 2 u (b0 , b1 ) u (b1 ) n1 b1
(12)
where η = (η0, η1)T is a vector variable describing possible values of the output quantity E. For a coverage probability p=0.95, a coverage factor kp = 2.45, is used.
[FIGURE 3]
ACCEPTED MANUSCRIPT Fig. 3 shows the marginal distributions for the slope (b1) and the intercept (b0) produced by the application of MCM (see section 5.3). Fig. 4 shows elliptical and rectangular coverage regions (coverage probability p=0.95) for the bivariate quantity E characterized by a
T
Gaussian probability distribution for which the component quantities b1 and b0 are
IP
correlated. 1000 draws were carried out (instead of 106) for illustration purposes.
SC R
As in our case the two outputs of the 1st stage of the measurement model (b1 and b0) are mutually correlated, the ellipse centered coverage region is more appropriate than the rectangular one. Generally the ellipse is the smallest coverage region for the stipulated
NU
probability, while the rectangle does not reflect at all the correlation between component quantities. Therefore, the rectangular coverage region might be considered inappropriate as a coverage region (the coverage probability for the rectangular region exceeds 0.95). A
MA
rectangle with sides parallel to the axes of the ellipse would have smaller area and might be considered more appropriate, but might be inconvenient since it would be expressed in
D
terms of variables that would be artificial in terms of the application.
TE
Fig. 5 depicts, in the form of contour lines, elliptical coverage regions (for various coverage probabilities p). The parameters kp that were used (Eq. 14) were calculated from the chi-
CE P
squared distribution so that:
p Pr( 22 k p2 )
AC
(13)
where 22 has a chi-squared distribution with two degrees of freedom.
[FIGURE 4]
[FIGURE 5]
5.3 Measurement uncertainty of the predicted value
ACCEPTED MANUSCRIPT
The algorithm of Monte Carlo Method was implemented in MATLAB® [19] using a fixed number of trials (M=106). MATLAB® code is presented in Table 4. Samples for b0 and b1 were
T
drawn from a bivariate (or joint) Gaussian distribution N(E,V) characterized by the
SC R
IP
expectation and the covariance (or uncertainty) matrices, E and V (Eqs 11 and 12).
Sample for y0 was drawn from rectangular (uniform) distribution R(y0-3u(y0), y0+3 u(y0)), where u(y0) is the standard measurement uncertainty of the instrument response that can
NU
be obtained from manufacturer specification which states a tolerance of ± 3%
MA
(u(y0)=tolerance/3).
[TABLE 4]
D
The results produced by the implementation of MCM were an estimate xpred =2.088 mg L-1 with an associated standard measurement uncertainty u(xpred) =0.043 mg L-1. Fig. 6 presents
TE
the PDF of xpred provided by MCM.
CE P
The 95% coverage interval for xpred can be obtained by the 2.5- and 97.5-percentiles of the distribution of MCM results (probabilistically equal coverage interval). The GUM supplement 1 [2] provides guidance for an alternative coverage, the “shortest coverage interval.” This
AC
interval is the shortest among all possible intervals for a distribution of MCM results where each interval has the same coverage probability (95%). The PDF obtained by MCM is bellshaped and symmetrical, which leads to identical shortest coverage interval and the probabilistically equal shortest interval [2.00 – 2.17] mg L-1.
[FIGURE 6] The results of the MCM were compared with the outcome of the equation giving the standard error of the estimate. The estimation of the uncertainty of a result from a linear calibration is very often approximated through the so called error of prediction, s(xpred) which is given by the following equation:
ACCEPTED MANUSCRIPT s( x pred )
SE regression
1 1 n N
b1
y
y
0
x n
2 1
b
i 1
2
x
i
(14)
2
T
where N is the number of repeated measurements made on the unknown test sample and
SC R
IP
y is the average of all measured responses.
The derivation of the Eq. (14) is described in detail by Hibbert [7]. It has to be noted that Eq. (14) actually calculates the combined uncertainty of the regression line and the uncertainty compound due to the repeatability of the response. The latter has a variance which is
MA
2 SEregression
N
(15)
D
Var ( y0 )
NU
estimated as:
TE
If we want to calculate the uncertainty of the regression line only, then the term 1/N in Eq.
SE regression b1
1 n
2 1
b
y
0
y
x n
i 1
i
2
x
(16)
2
AC
s' ( x pred )
CE P
(14) should be omitted. This leads to the expression:
The above expression does not take into account the uncertainty u(y0) of the response y0 which is an input parameter of the measurement model (Eq. (6)). This may lead to results underestimated and not comparable with the results of MCM. The standard measurement uncertainty of the predicted value including the uncertainty of the response is given by the equation:
1 u ( xpred ) s' ( xpred ) u ( y 0 ) b1 2
2
(17)
The application of the Eq. (16) and (17) to our calibration data provides an estimate of the standard measurement uncertainty u(xpred) equal to 0.043 mg L-1 which agrees well with the results of the MCM.
ACCEPTED MANUSCRIPT
T
6 Conclusions
IP
Measurement models, consisting of successive distinct stages (multistage measurement
SC R
models) and having more than one output quantity are receiving increasing attention concerning measurement uncertainty evaluation. The GUM uncertainty framework and the Monte Carlo Method (MCM) can be extended to include such measurement models. The construction and use of a calibration curve in order to determine the value of a measurand
NU
(e.g. analyte concentration) is an integral part of many analytical procedures and is actually a typical multistage measurement model, that also includes a measurement stage with more
MA
than one output quantity.
In the present work the MCM was used to evaluate the component of measurement
D
uncertainty from a calibration curve used for the determination of the total nitrogen found in a gasoline sample. The parameters of the calibration curve (slope and intercept) were
TE
calculated using linear regression. The estimated regression statistics showed a strong linear relationship between concentration and instrument response values and a regression line
CE P
that can reasonably be assumed to pass through zero. The slope and the intercept of the calculated calibration function were treated as a vector output quantity characterized by a joint probability distribution (bivariate Gaussian). Two
AC
types of coverage region were considered in order represent the vector quantity value with a stated probability (95%): a rectangle centered coverage region and an ellipse centered coverage region. Because of the mutual correlation between slope and intercept, the ellipse centered coverage region is more appropriate than the rectangular one.
The vector output quantity of the slope and intercept and the associated covariance matrix were then used as an input to the second stage of the measurement model which is the use of the calibration function for the estimation the quantity value (nitrogen content) corresponding to the unknown gasoline sample indication. The predicted value, its measurement uncertainty and a 95% coverage region were evaluated using MCM employing fixed number of trials (106). The algorithm of the MCM, implemented in MATLAB, gave a predicted value of 2.09 mg L-1, a standard measurement uncertainty of 0.04 mg L-1 and a
ACCEPTED MANUSCRIPT 95% symmetrical coverage interval [2.01 – 2.17] mg L-1. The results of the MCM agreed well
IP
T
with the outcome of the equation giving the standard error of the estimate.
SC R
References
[1] JCGM 100:2008, Guide to the expression of uncertainty in measurement (GUM 1995 with minor corrections), Joint Committee for Guides in Metrology (2009)
NU
[2] JCGM 101:2008, Evaluation of measurement data - Supplement 1 to the “Guide to the expression of uncertainty in measurement" - Propagation of distributions using a
MA
Monte Carlo method, Joint Committee for Guides in Metrology, 2008. [3] JCGM 102:2011, Evaluation of measurement data – Supplement 2 to the "Guide to the expression of uncertainty in measurement" – Extension to any number of output
D
quantities, Joint Committee for Guides in Metrology, 2011.
TE
[4] R. Bro, A. Rinnan, N.M. Faber, Standard error of prediction for multilinear PLS 2. Practical implementation in fluorescence spectroscopy, Chemometrics and
CE P
Intelligent Laboratory Systems 75 (2005) 69–76. [5] E. Saari, P. Perämäki, J. Jalonen, Measurement uncertainty in the determination of total petroleum hydrocarbons (TPH) in soil by GC-FID, Chemometrics and Intelligent
AC
Laboratory Systems 92 (2008) 3–12. [6] Carotenuto, G. Giovinco, B. Viglietti, L. Vanoli, A new procedure for the determination of calibration curves for a gas chromatograph used in natural gas analysis, Chemometrics and Intelligent Laboratory Systems 75 (2005) 209– 217
[7] D.B. Hibbert, The uncertainty of a result from a linear calibration, Analyst 131 (2006) 1273–1278. [8] V.J. Barwick, Preparation of calibration curves—a guide to best practice, VAM, Teddington UK, 2003. [9] L. Cuadros-Rodríguez, L. Gámiz-Gracia, E. Almansa-López, J. Laso-Sánchez Calibration in chemical measurement processes: I. A metrological approach, TrAC - Trends in Analytical Chemistry 20 (2001) 195–206. [10]D. Theodorou, L. Meligotsidou, S. Karavoltsos, A. Burnetas, M. Dassenakis, M. Scoullos, Comparison of ISO-GUM and Monte Carlo methods for the evaluation of
ACCEPTED MANUSCRIPT measurement uncertainty: Application to direct cadmium measurement in water by GFAAS, Talanta 83(2011) 1568-1574 [11]ASTM D4629-12, Standard Test Method for Trace Nitrogen in Liquid Petroleum
T
Hydrocarbons by Syringe/Inlet Oxidative Combustion and Chemiluminescence
IP
Detection, ASTM International, West Conshohocken, PA, 2012, www.astm.org [12]A. Possolo, Copulas for uncertainty analysis, Metrologia 47 (2010) 262-271.
SC R
[13]JCGM 104:2009. Evaluation of measurement data - An introduction to the “Guide to the expression of uncertainty in measurement" and related documents, Joint Committee for Guides in Metrology, 2009.
NU
[14]M.A. Herrador, A.G. Asuero, A.G. González, Estimation of the uncertainty of indirect measurements from the propagation of distributions by using the Monte-Carlo method: An overview, Chemometrics and Intelligent Laboratory Systems 79 (2005)
MA
115–122
[15]P.M. Harris, M.G. Cox, On a Monte Carlo method for measurement uncertainty evaluation and its implementation, Metrologia 51 (2014) S176-S182.
D
[16]W. Bich, M. G. Cox, P.M. Harris Evolution of the 'Guide to the Expression of
TE
Uncertainty in Measurement', Metrologia 43 (2006) S161–S166 [17]D. Theodorou, Y. Zannikou, G. Anastopoulos, F. Zannikos, Coverage interval
CE P
estimation of the measurement of Gross Heat of Combustion of fuel by bomb calorimetry: Comparison of ISO GUM and adaptive Monte Carlo Method, Thermochimica Acta 526 (2011) 122– 129
AC
[18]M. Solaguren-Beascoa Fernández, J.M. Alegre Calderón, P.M. Bravo Díez, Implementation in MATLAB of the adaptive Monte Carlo method for the evaluation of measurement uncertainties, Accreditation and Quality Assurance 14 (2009) 95– 106. [19]MATLAB® and Statistics Toolbox Release 2012b, The MathWorks, Inc., Natick, Massachusetts, United States.
ACCEPTED MANUSCRIPT Quantity values, xi, provided by measurement (calibration) standards and corresponding instrument indications (‘signals’, yi)
1st stage of the measurement model Measurement model with multiple output quantities
IP
T
Estimation of the parameters (e.g slope and intercept) of the calibration function (relation between the indication and the quantity value)
SC R
2nd stage of the measurement model
Analysis of an unknown sample, indication y0
NU
Estimation of the quantity value, xpred, corresponding to the unknown sample, indication, y0
MA
Fig. 1 Two-stage measurement model involving the estimation of the parameters of a calibration function (1st stage of the measurement model, having multiple output quantities) and the use of a calibration function for the estimation of a quantity value (2nd stage of the
AC
CE P
TE
D
measurement model)
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Fig. 2 Calibration graph describing relationship between nitrogen content and instrument
AC
CE P
TE
D
response
ACCEPTED MANUSCRIPT
NU
SC R
IP
T
[20]
AC
CE P
TE
D
MA
Fig. 3 Marginal distributions of the slope (b1) and the intercept (b0)
ACCEPTED MANUSCRIPT
T
5450
IP SC R
slope (L/mg)
5400
5350
NU
5300
-200
-100
0
100
MA
5250
200 300 intercept (AU)
400
500
600
700
800
Fig. 4 Elliptical and rectangular coverage regions (coverage probability p=0.95) of the joint
AC
CE P
TE
D
probability density function of the slope, b1 and the intercept,b0.
ACCEPTED MANUSCRIPT
Data p=0.500 p=0.750 p=0.900 p=0.950 p=0.990 p=0.999
T
5450
IP
5350
SC R
slope (L/mg)
5400
5300
-200
0
200 400 intercept (AU)
MA
-400
NU
5250
600
800
1000
Fig. 5 Elliptical coverage regions (for various coverage probabilities p) of the joint probability
AC
CE P
TE
D
density function of the slope, b1 and the intercept,b0.
ACCEPTED MANUSCRIPT 8
T
7
IP
6
SC R
Probability density
5
4
NU
3
1
0 1,90
2,09 nitrogen content (mg/L)
2,17
2,25
TE
D
2,01
MA
2
Fig. 6 Probability density function of xpred provided by Monte Carlo Method. The vertical
AC
CE P
dashed lines indicate the 95% probabilistically symmetrical coverage interval.
ACCEPTED MANUSCRIPT Table 1 Instrument parameters used for nitrogen determination Value
Volume injected (μL)
40
Syringe drive rate (μL s-1)
1
IP
T
Parameter
Furnace temperature (oC) Inlet oxygen flowmeter setting (mL min-1)
SC R
1050
AC
CE P
TE
D
MA
NU
Inlet carrier (Argon) flowmeter setting (mL min-1)
100 100
ACCEPTED MANUSCRIPT Table 2 Calibration and sample measerement raw data Calibration standard response (AU) replicate
replicate replicate
(mg L-1)
1
2
0.0
557.8
546.5
0.1
966.5
0.5
3057
1.0
5304
5353
5324
5.0
25890
26490
27850
53740
54840
53310
CE P
TE
IP 447.8
SC R
Unknown sample response (AU)
AC
3
895.8
927.6
2963
2868
NU
MA D
10.0
T
Concentration
11440
ACCEPTED MANUSCRIPT Table 3 Results of least square linear regression Regression parameter
Value
Estimate of the intercept, b0
263.99 AU 5352.99 AU mg-1 L
T
Estimate of the slope, b1
149.94 AU
IP
Standard error of the estimate of the intercept, SE(b0)
32.69 AU mg-1 L
SC R
Standard error of the estimate of the slope, SE(b1) Coefficient of determination, r2 standard
error
of
regression
(residual standard deviation), SEregression
value (p-value) Degrees of freedom, df
D
Regression sum of squares, SSregression
MA
The F statistic or Fisher – Snedecor F-
AC
CE P
TE
The residual sum of squares, SSresidual
507.4 AU
NU
The
0.9994
26821 (3.1 10 -27)
16 6.91 109 AU2 4.12 106 AU2
ACCEPTED MANUSCRIPT Table 4 Implementation program of Monte Carlo Method in MATLAB®
Line
Matlab Code function [x_pred]=regression_chem_compact(M)
2
% Calibration Data
3
x1=0; x2=0.1; x3=0.5; x4=1; x5=5; x6=10; x7=0; x8=0.1; x9=0.5; x10=1; x11=5; x12=10; x13=0; x14=0.1; x15=0.5; x16=1; x17=5; x18=10;
4
y1=557.8; y2=966.5; y3=3057.0; y4=5304; y5=25890; y6=53740; y7=546.5; y8=895.8; y9=2963; y10=5353; y11=26490; y12=54840; y13=447.8; y14=927.6; y15=2868; y16=5324; y17=27850; y18=53310;
5
%mean of all concentrations
6
x_av=(x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+x15+x16+x17+x18)/18;
7
% mean of all responses
8
y_av=(y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y11+y12+y13+y14+y15+y16+y17+y18)/18;
9
% Obtain estimates for the slope and the intercept
10
Xi=[x1; x2; x3; x4; x5; x6; x7; x8; x9; x10; x11; x12; x13; x14; x15; x16; x17; x18];
11
Yi=[y1; y2; y3; y4; y5; y6; y7; y8; y9; y10; y11; y12; y13; y14; y15; y16; y17; y18];
12
b1=sum((Xi-x_av).*(Yi-y_av))./sum((Xi-x_av).^2)
13
b0=y_av-b1.*x_av
14
%standard errors of slope and intercept
15
u_b1=32.69;
16
u_b0=149.94;
17
% simulate b1,b0 from joint pdf
18
r_b1b0=-sum(Xi)/sqrt(sumsqr(Xi)*numel(Xi)); % correlation coefficient
19
u_b1b0=u_b1*u_b0*r_b1b0;
20
%covariance matrix
21
V=[u_b1^2 u_b1b0; u_b1b0 u_b0^2];
22
%mean matrix
23
E=[b1; b0];
24
X=mvnrnd(E,V,M); % Mx2 matrix
25
B1=X(:,1); % 1st column of Mx2 matrix
26
B0=X(:,2); % 2nd column of Mx2 matrix
27
% simulate sample of y0
28
[a,b]=unif_param(11440,198.1466124);
AC
CE P
TE
D
MA
NU
SC R
IP
T
1
ACCEPTED MANUSCRIPT Line
Matlab Code y0=unifrnd(a,b,M,1);
30
% obtain x_pred
31
x_pred=(y0-B0)./B1;
32
% Results
33
x_pred=sort(x_pred);
34
Mean_Value=mean(x_pred)
35
Standard_Uncertainty=std(x_pred)
36
figure;
37
lower=min(x_pred);
38
upper=max(x_pred);
39
xc=lower:(upper-lower)/499:upper;
40
n=hist(x_pred,xc);
41
bar(xc,n./(((upper-lower)/499)*M),1);
42
hold on;
43
%ESTIMATE AND PLOT COVERAGE REGIONS FOR b1, b0
44
data=[B0 B1];
45
figure; % 95% rectangular and elliptical coverage region
46
%Plot the original data
47
plot(data(:,1), data(:,2), 'k.');
48
xlim([min(B0)-50, max(B0)+50]);
49
ylim([min(B1)-50, max(B1)+50]);
50
hold on;
51
% Set the axis labels
52
hXLabel = xlabel('b0');
53
hYLabel = ylabel('b1');
54
% Plot 95% rectangular coverage region
55
rectangle('Position',[b0-2.24*u_b0, b1-2.24*u_b1, 2*2.24*u_b0, 2*2.24*u_b1], 'Linestyle', '--')
56
% Plot 95% elliptical coverage region
57
% Calculate the eigenvectors and eigenvalues
58
covariance = cov(data);
59
[eigenvec, eigenval ] = eig(covariance);
60
% Get the index of the largest eigenvector
AC
CE P
TE
D
MA
NU
SC R
IP
T
29
ACCEPTED MANUSCRIPT Line
Matlab Code [largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
62
largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);
63
% Get the largest eigenvalue
64
largest_eigenval = max(max(eigenval));
65
% Get the smallest eigenvector and eigenvalue
66
if(largest_eigenvec_ind_c == 1)
SC R
67
smallest_eigenval = max(eigenval(:,2));
68
smallest_eigenvec = eigenvec(:,2); else
NU
69
IP
T
61
smallest_eigenval = max(eigenval(:,1));
71
smallest_eigenvec = eigenvec(1,:);
MA
70
end
73
% Calculate the angle between the x-axis and the largest eigenvector
74
angle = atan2(largest_eigenvec(2), largest_eigenvec(1));
75
% This angle is between -pi and pi.
76
% Shift it such that the angle is between 0 and 2pi
77
if(angle < 0)
angle = angle + 2*pi;
CE P
78
TE
D
72
end
80
% Get the coordinates of the data mean
81
avg = mean(data);
82
% Get the 95% confidence interval error ellipse
83
chisquare_val = 2.4477; % set chisquare_val = 1.1774 for 50%, 1.6651 for 75%, 2.1460 for 90%, 3.0349 for 99%, 3.7169 for 99.9% confidence levels
84
theta_grid = linspace(0,2*pi);
85
phi = angle;
86
X0=avg(1);
87
Y0=avg(2);
88
a=chisquare_val*sqrt(largest_eigenval);
89
b=chisquare_val*sqrt(smallest_eigenval);
90
% the ellipse in x and y coordinates
91
ellipse_x_r
= a*cos( theta_grid );
92
ellipse_y_r
= b*sin( theta_grid );
AC
79
ACCEPTED MANUSCRIPT Line
Matlab Code %Define a rotation matrix
94
R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];
95
%Rotate the ellipse to some angle phi
96
r_ellipse = [ellipse_x_r;ellipse_y_r]' * R;
97
% Draw the error ellipse
98
plot(r_ellipse(:,1) + X0,r_ellipse(:,2) + Y0,'k-')
99
figure % marginal probabilities of b0, b1
SC R
IP
T
93
100 subplot(1,2,1);
NU
101 B1=sort(B1); 102 Mean_value=mean(B1);
MA
103 lower=min(B1); 104 upper=max(B1);
105 xc=lower:(upper-lower)/199:upper;
D
106 n=hist(B1,xc);
108 title('b1')
CE P
109 subplot(1,2,2);
TE
107 bar(xc,n./(((upper-lower)/199)*M),1);
110 B0=sort(B0);
111 Mean_value=mean(B0); 112 lower=min(B0);
AC
113 upper=max(B0);
114 xc=lower:(upper-lower)/199:upper; 115 n=hist(B0,xc);
116 bar(xc,n./(((upper-lower)/199)*M),1); 117 title('b0') 118 end 119 function [a,b]=unif_param(m,v) 120 % A function to obtain the parameters of the uniform (rectangular) 121 % distribution given its mean and standard deviation 122 b=m+v*sqrt(12)/2; 123 a=2*m-b; 124 end
ACCEPTED MANUSCRIPT Highlights •
The measurement uncertainty of the nitrogen content in gasoline determined using a calibration curve is evaluated This is a typical multistage measurement model with more than one output quantity
•
The coefficients of the calibration curve are characterized by a joint probability
IP
T
•
Monte Carlo Method is used to evaluate the calibration curve uncertainty
CE P
TE
D
MA
NU
component
AC
•
SC R
distribution
S0169-7439(15)00145-8 doi: 10.1016/j.chemolab.2015.05.025 CHEMOM 3028
To appear in:
Chemometrics and Intelligent Laboratory Systems
Received date: Revised date: Accepted date:
10 February 2015 13 May 2015 22 May 2015
Please cite this article as: Dimitrios Theodorou, Ypatia Zannikou, Fanourios Zannikos, Component of measurement uncertainty from a measurement model with two stages involving two output quantities, Chemometrics and Intelligent Laboratory Systems (2015), doi: 10.1016/j.chemolab.2015.05.025
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT TITLE:
Component of measurement uncertainty from a measurement model
AUTHORS:
SC R
Dimitrios Theodorou, Ypatia Zannikou, Fanourios Zannikos
IP
T
with two stages involving two output quantities
INSTITUTION:
National Technical University of Athens, School of Chemical Engineering, Laboratory of Fuels
NU
and Lubricants Technology, Iroon Polytechneiou 9, Athens 15780, Greece
MA
Author Contact Details Dimitrios Theodorou
Tel./fax: +30 21 0652 5976.
(CORRESPONDING AUTHOR)
E-mail address: [email protected]
AC
CE P
Ypatia Zannikou
TE
D
Permanent Address: 59 A.Tritsi Street, 15238 Halandri, Athens, Greece.
Fanourios Zannikos
Tel 210 7723143 / Fax 210 7723163 E-mail address: [email protected]
Tel 210 7723143 / Fax 210 7723163 E-mail address: [email protected]
Abstract The GUM uncertainty framework and the Monte Carlo Method can be extended to include multistage measurement models and measurement models having more than one output quantity. A typical multistage measurement model, that also includes a measurement stage
ACCEPTED MANUSCRIPT with more than one output quantity is the calculation of the coefficients of a calibration curve, constructed by means of least squares regression and the subsequent use this calibration curve for the estimation of a quantity value. In the present work the Monte Carlo
T
Method was used to evaluate the component of measurement uncertainty from a
IP
calibration curve used for the determination of the total nitrogen found in a gasoline sample. The slope and the intercept of the calibration curve were treated as a vector output
SC R
quantity characterized by a joint probability distribution (bivariate Gaussian) and two types of coverage regions were estimated (rectangular and ellipsoid). The Monte Carlo Method algorithm employing fixed number of trials (106) gave a predicted value of 2.09 mg L-1, a
NU
standard measurement uncertainty of 0.04 mg L-1 and a 95% symmetrical coverage interval [2.01 – 2.17] mg L-1. The standard measurement uncertainty obtained by MCM agrees well
MA
with the outcome of the equation giving the standard error of the estimate.
TE
D
1. Introduction
Until recently there has been a focus on the evaluation of measurement uncertainty of
CE P
simple measurement models consisting of one stage and having a single output quantity (univariate scalar models). There are measurement models, however, where: (i) the output quantities from previous stages become input quantities for subsequent stages (multistage measurement models) and, (ii) there is more than one output quantity, depending on a
AC
common set of input quantities. The principles for the measurement uncertainty evaluation described in the “Guide to Expression of Uncertainty in Measurement” [1] (GUM uncertainty framework) and its Supplement [2] describing the application of Monte Carlo Method (MCM) can be extended to the evaluation of uncertainties associated with vector estimates of multivariate output quantities. In 2011, JCGM published the GUM supplement 2 [3] which provides guidance concerning the application of a generalized GUM uncertainty framework and MCM in measurement models with any number of output quantities. A typical multistage measurement model, that also includes a measurement stage with more than one output quantity is the calculation of the coefficients of a calibration curve, constructed by means of least squares regression and the subsequent use of this calibration curve for the estimation of a quantity value. In analytical chemistry, linear regression is very often used in the construction of calibration curves required for techniques such as liquid or
ACCEPTED MANUSCRIPT gas chromatography and atomic absorption or fluorescence spectrometry [4-6]. A calibration curve is actually the expression of the relation between indication and corresponding measured quantity value. Calibration establishes a relationship between the values of a
T
standard (reference values), xi and the output quantities (response of the instrument), yi.
IP
Using values of a standard (reference values), xi and instrument responses, yi, a relationship (often assumed to be represented by a straight line) is established and thenthe calibration
SC R
model is used in reverse, that is, to predict a value, xpred from an instrument response, y0 [79]. Fig. 1 illustrates the two-stage measurement model involving construction and use of a calibration function. As with most statistics, calibration curve coefficients (slope and
NU
intercept for linear calibration model) are only estimates based on a finite number of measurements, and therefore, their values are associated with uncertainties. This leads to an measurement uncertainty of the predicted value as well. In many cases the component of
MA
measurement uncertainty from the calibration curve contributes significantly to the
[FIGURE 1]
CE P
TE
D
measurement uncertainty [10] and therefore it has to be reliably evaluated.
The determination of the total nitrogen found in liquid hydrocarbons is conducted using oxidative combustion and chemiluminescence detection method [11], which requires the construction of a linear calibration curve. This determination is important and has to be
AC
reliable because even trace amounts of nitrogenous materials contained in the feedstocks may lead to poisoning of process catalysts used in petroleum and chemical refining. In the present work, the measurement uncertainty of the nitrogen content in liquid fuel products determined using a calibration curve is evaluated using MCM. The coefficients (slope and intercept) of the calibration curve are treated as a vector output quantity characterized by a joint probability distribution [12]. Bivariate coverage regions are specified taking the form an ellipsoid or a rectangle.
ACCEPTED MANUSCRIPT 2. Experimental
An Analytik Jena AG - multi EA® 3100 Element Analyzer equipped with automatic sampler
T
was employed in this work. This analyzer fully complies with ASTM D 4629-12 “Standard Test
IP
Method for Trace Nitrogen in Liquid Petroleum Hydrocarbons by Syringe/Inlet Oxidative
SC R
Combustion and Chemiluminescence Detection” [11]. Table 1 presents the operating conditions of the instrument. Six calibration standards (Analytik Jena nitrogen calibration kit) with concentrations 0.0, 0.1, 0.5, 1.0, 5.0 and 10.0 mg L-1 were analyzed in triplicate. The instrument responses (AU – Area Units) were recorded and a calibration curve with 18
NU
points was constructed. Then a gasoline sample was measured and its total nitrogen concentration in mg L-1 was predicted using the calibration curve constructed. The sample
MA
and the calibration standards were injected using a 40 μL syringe. Data produced are
CE P
TE
D
presented in Table 2.
[TABLE 1]
[TABLE 2]
AC
3. Linear regression by least-squares method
The least-squares method is used for data evaluation problems related to measurement models with an underlying theoretical relationship between an independent variable and a dependent variable. This relationship constitutes the basis of a curve-fitting problem. The task is to estimate the parameters of a calibration function from pairs of measured quantity values and their corresponding indication values [13].
The measured quantity values of an independent variable are typically those of measurement standards. The values of the dependent variables are the indication values returned by the measuring system for the corresponding values of the independent variables. In chemistry, the establishment of a calibration function has the purpose of obtaining a function that allows one to calculate the concentration of the analyte as a
ACCEPTED MANUSCRIPT function of an instrumental signal. A set of standards are prepared containing a known amount of the analyte of interest and the instrument response of each standard is measured. Then a relationship is established between the instrument response and the
T
analyte concentration. The most frequent type of calibration in chemical analyses is the
IP
indirect calibration, where the equipment gives a value (signal or instrument response),
SC R
which has a different quantity from that of the standard (reference value) [7].
Calibration data (pairs of analyte concentrations, xi, and instrument responses, yi ) are used as an input to the least square regression analysis which fits a predetermined measurement
NU
model to these data. The simplest measurement model is the one described by the linear
MA
function:
(1)
D
Y b0 b1 x
TE
where Y is the instrument response (dependent variable), x is the analyte concentration (independent variable) and b0 and b1 are the coefficients of the model known as the
CE P
intercept and slope, respectively. The slope and the intercept are estimated from the calibration data (xi, yi).
It has to be noted that least square regression process can validly be applied under the following conditions [7,8]:
AC
a. The linear model holds for the data (i.e. the response does indeed vary linearly with concentration). b. The errors in the independent variable are insignificant compared with those of the dependent variable. c. These errors are normally distributed and independent of the values of the dependent variable (i.e. the data are homoscedastic). Only when these conditions hold, the application of linear least square regression is expected to realize the best fit to the calibration data.
The quantities b0 and b1 are correlated and the degree of their correlation is characterized by the estimated correlation coefficient r(b0, b1):
ACCEPTED MANUSCRIPT n
x
r (b0 , b1 )
i 1
i
(2)
n
n xi2
T
i 1
IP
where n is the number calibration data (pairs of xi, yi).
SC R
Once the slope and the intercept have been calculated, it becomes possible to calculate the standard error of regression (residual standard deviation) which is required for the calculation of a number of regression statistics. The standard error of regression is a
NU
measure of the deviation of the data points from the regression line and is calculated using
MA
the equation:
yi yˆ i 2
SEregression
(3)
TE
D
n2
where yˆ (predicted y value) is the value of y calculated for a given x value using the linear
CE P
calibration equation.
The standard deviation for the estimate of the intercept, SE(b0), is calculated using the
AC
equation:
SE (b0 ) SEregression
n
x i 1
2 i
(4)
n
n ( xi x )
2
i 1
where x is the average concentration of all standards used. The standard deviation for the estimate of the slope, SE(b1), is calculated using the equation:
SE (b1 )
SEregression n
(x i 1
i
(5)
x)2
Once the least-squares method has been used to estimate the parameters of a calibration function and evaluate the associated standard uncertainties and covariances, the measuring
ACCEPTED MANUSCRIPT system is subsequently used for a measurement. When analyzing unknown test samples the relationship of Eq. (1) is used in reverse in order to predict a value from an instrument
IP
y 0 b0 b1
(6)
SC R
x pred
T
response:
where y0 is the instrument response of the unknown test sample and xpred is the prediction (or estimation) of the measurement model concerning the true value of the measurand (e.g. analyte concentration).
NU
The standard measurement uncertainty associated with the predicted value estimate is evaluated using the standard uncertainties and covariances associated with the calibration
MA
function parameter estimates and the standard measurement uncertainty associated with the indication value.
TE
D
4. Measurement uncertainty evaluation
The result obtained from a calibration curve is actually an estimate of the true value of the
CE P
measurand. Therefore, it is associated with a measurement uncertainty due to lack of fit of the linear model to the data, which leads to uncertainties (standard errors) of the parameters b0 and b1.
AC
The GUM and its supplements [1-3] provide guidance on measurement uncertainty evaluation using as information the best estimates and the standard uncertainties of the input quantities (b0, b1, y0). One approach is to propagate this information, using the law of propagation of uncertainty through a first-order Taylor series approximation to the measurement model to provide the estimate of the output quantity xpred and the standard measurement uncertainty u(xpred) associated with xpred (GUM uncertainty framework) [13]. It has to be noted though that, the GUM uncertainty framework might not be satisfactory, in situations, where: (i)
the measurement function is non-linear,
(ii)
the probability distributions for the input quantities are asymmetric,
(iii)
the measurement uncertainty contributions are not of approximately the same magnitude,
ACCEPTED MANUSCRIPT (iv)
the probability distribution for the output quantity is either asymmetric, or not a Gaussian or a t-distribution.
T
An alternative method for evaluating measurement uncertainty is the Monte Carlo method
IP
(MCM). This method is described in the first supplement of GUM [2] and involves no restrictions for valid application [14,15] except for the fact that it requires prior knowledge
SC R
concerning the type of the probability distribution assigned to each input parameter. The MCM actually combines and propagates distributions rather than propagating uncertainties as in the GUM uncertainty framework. The measurement uncertainty evaluation is based on
NU
a probabilistic approach that combines the whole distribution of the input parameters and is not just based on their best estimates and standard measurement uncertainties.
MA
The MCM consists of simulating draws from the distribution of the output quantity (xpred), based on simulated draws from the distributions of the input quantities (b0, b1, y0). In our case, the quantities b0 and b1 (output quantities for the 1st stage of the measurement model
D
and input quantities for the second stage of the measurement model) are not independent
TE
(non zero covariance) and therefore they should be described by a multivariate (joint) distribution. This distribution is used for the draws during the Monte Carlo Method
CE P
implementation. The number of trials selected for MCM determines the degree of numerical accuracy of the MCM outputs [16]. An adaptive procedure by which the number of trials can be chosen, may be used [17,18]. Nevertheless, the Monte Carlo numerical simulation tends
AC
to require up to 106 trials for calculating a standard measurement uncertainty or a coverage interval which is correct to one or two significant decimals digits.
The estimate of the output quantity xpred is estimated by the average of the M MCM trials which produce M measurement model values (x pred (k) , k = 1,…,M):
x pred
1 M
M
x k 1
(k ) pred
(7)
while the standard measurement uncertainty u(xpred) associated with xpred is evaluated as the standard deviation of the M model values:
ACCEPTED MANUSCRIPT 1 M (k ) ( x pred x pred ) M 1 k 1
(8)
T
u ( xpred )
SC R
IP
5. Results and discussion
5.1 Calibration results
NU
A least square linear regression was carried out for the calibration data (see Table 2) using the LINEST function of Microsoft Excel. The results of the regression analysis are presented
MA
in Table 3. The coefficient of determination for the straight line, r2, is 0.9994, and the regression ANOVA gave a Fisher – Snedecor value F equal to 26821 with a significance (p – value) of 3.1 10-27. This implies a strong linear relationship between concentration and
D
instrument response values. Fig. 2 shows the plot produced by the calibration data and the
[TABLE 3]
[FIGURE 2]
AC
CE P
TE
fitted line.
5.2 Measurement uncertainty of calibration function parameters
The estimation of the slope and the intercept of the calibration curve is based on a measurement model of multiple outputs (two in our case), which involves correlated data. It is more convenient to treat the coeffiecients of the calibration function as a vector quantity E (Eq. 11), associated with a covariance matrix V (Eq. 12)
b E 0 b1
(9)
ACCEPTED MANUSCRIPT u 2 (b0 ) u (b0 , b1 ) SE 2 (b0 ) r (b0 , b1 ) SE (b0 ) SE (b1 ) V 2 SE 2 (b1 ) u (b0 , b1 ) u (b1 ) r (b0 , b1 ) SE (b0 ) SE (b1 )
(10)
where u (b0 , b1 ) = u (b1 , b0 ) is the estimated covariance between b0 and b1. The standard
T
uncertainties u(b0) and u(b1) are provided by the regression statistics as standard error of
SC R
IP
intercept, SE(b0), and slope, SE(b1), respectively.
Once the probability density function (PDF) is established (bivariate Gaussian) a coverage region can be determined, for a coverage probability, p. A coverage region specifies a region
NU
in 2-dimensional space that contains E with probability p. Two types of coverage region can be considered [3]:
rectangle centered coverage region (separately determined coverage intervals for
MA
b1 and b0).
ellipse centered coverage region
D
TE
The rectangle centered coverage region (sides parallel to the axes) is determined by the marginal PDFs of b1 and b0 and the coverage intervals are: b1 ± kq u(b1) and b0 ± kq u(b0). For
CE P
coverage probability p=0.95, a coverage factor kq = 2.24, is used.
AC
The ellipse centered coverage region is determined as: (η – E)T V-1 (η – E) = kp2
or
n0 b0
n1 b1
(11)
1
u 2 (b0 ) u (b0 , b1 ) n0 b0 k p2 2 u (b0 , b1 ) u (b1 ) n1 b1
(12)
where η = (η0, η1)T is a vector variable describing possible values of the output quantity E. For a coverage probability p=0.95, a coverage factor kp = 2.45, is used.
[FIGURE 3]
ACCEPTED MANUSCRIPT Fig. 3 shows the marginal distributions for the slope (b1) and the intercept (b0) produced by the application of MCM (see section 5.3). Fig. 4 shows elliptical and rectangular coverage regions (coverage probability p=0.95) for the bivariate quantity E characterized by a
T
Gaussian probability distribution for which the component quantities b1 and b0 are
IP
correlated. 1000 draws were carried out (instead of 106) for illustration purposes.
SC R
As in our case the two outputs of the 1st stage of the measurement model (b1 and b0) are mutually correlated, the ellipse centered coverage region is more appropriate than the rectangular one. Generally the ellipse is the smallest coverage region for the stipulated
NU
probability, while the rectangle does not reflect at all the correlation between component quantities. Therefore, the rectangular coverage region might be considered inappropriate as a coverage region (the coverage probability for the rectangular region exceeds 0.95). A
MA
rectangle with sides parallel to the axes of the ellipse would have smaller area and might be considered more appropriate, but might be inconvenient since it would be expressed in
D
terms of variables that would be artificial in terms of the application.
TE
Fig. 5 depicts, in the form of contour lines, elliptical coverage regions (for various coverage probabilities p). The parameters kp that were used (Eq. 14) were calculated from the chi-
CE P
squared distribution so that:
p Pr( 22 k p2 )
AC
(13)
where 22 has a chi-squared distribution with two degrees of freedom.
[FIGURE 4]
[FIGURE 5]
5.3 Measurement uncertainty of the predicted value
ACCEPTED MANUSCRIPT
The algorithm of Monte Carlo Method was implemented in MATLAB® [19] using a fixed number of trials (M=106). MATLAB® code is presented in Table 4. Samples for b0 and b1 were
T
drawn from a bivariate (or joint) Gaussian distribution N(E,V) characterized by the
SC R
IP
expectation and the covariance (or uncertainty) matrices, E and V (Eqs 11 and 12).
Sample for y0 was drawn from rectangular (uniform) distribution R(y0-3u(y0), y0+3 u(y0)), where u(y0) is the standard measurement uncertainty of the instrument response that can
NU
be obtained from manufacturer specification which states a tolerance of ± 3%
MA
(u(y0)=tolerance/3).
[TABLE 4]
D
The results produced by the implementation of MCM were an estimate xpred =2.088 mg L-1 with an associated standard measurement uncertainty u(xpred) =0.043 mg L-1. Fig. 6 presents
TE
the PDF of xpred provided by MCM.
CE P
The 95% coverage interval for xpred can be obtained by the 2.5- and 97.5-percentiles of the distribution of MCM results (probabilistically equal coverage interval). The GUM supplement 1 [2] provides guidance for an alternative coverage, the “shortest coverage interval.” This
AC
interval is the shortest among all possible intervals for a distribution of MCM results where each interval has the same coverage probability (95%). The PDF obtained by MCM is bellshaped and symmetrical, which leads to identical shortest coverage interval and the probabilistically equal shortest interval [2.00 – 2.17] mg L-1.
[FIGURE 6] The results of the MCM were compared with the outcome of the equation giving the standard error of the estimate. The estimation of the uncertainty of a result from a linear calibration is very often approximated through the so called error of prediction, s(xpred) which is given by the following equation:
ACCEPTED MANUSCRIPT s( x pred )
SE regression
1 1 n N
b1
y
y
0
x n
2 1
b
i 1
2
x
i
(14)
2
T
where N is the number of repeated measurements made on the unknown test sample and
SC R
IP
y is the average of all measured responses.
The derivation of the Eq. (14) is described in detail by Hibbert [7]. It has to be noted that Eq. (14) actually calculates the combined uncertainty of the regression line and the uncertainty compound due to the repeatability of the response. The latter has a variance which is
MA
2 SEregression
N
(15)
D
Var ( y0 )
NU
estimated as:
TE
If we want to calculate the uncertainty of the regression line only, then the term 1/N in Eq.
SE regression b1
1 n
2 1
b
y
0
y
x n
i 1
i
2
x
(16)
2
AC
s' ( x pred )
CE P
(14) should be omitted. This leads to the expression:
The above expression does not take into account the uncertainty u(y0) of the response y0 which is an input parameter of the measurement model (Eq. (6)). This may lead to results underestimated and not comparable with the results of MCM. The standard measurement uncertainty of the predicted value including the uncertainty of the response is given by the equation:
1 u ( xpred ) s' ( xpred ) u ( y 0 ) b1 2
2
(17)
The application of the Eq. (16) and (17) to our calibration data provides an estimate of the standard measurement uncertainty u(xpred) equal to 0.043 mg L-1 which agrees well with the results of the MCM.
ACCEPTED MANUSCRIPT
T
6 Conclusions
IP
Measurement models, consisting of successive distinct stages (multistage measurement
SC R
models) and having more than one output quantity are receiving increasing attention concerning measurement uncertainty evaluation. The GUM uncertainty framework and the Monte Carlo Method (MCM) can be extended to include such measurement models. The construction and use of a calibration curve in order to determine the value of a measurand
NU
(e.g. analyte concentration) is an integral part of many analytical procedures and is actually a typical multistage measurement model, that also includes a measurement stage with more
MA
than one output quantity.
In the present work the MCM was used to evaluate the component of measurement
D
uncertainty from a calibration curve used for the determination of the total nitrogen found in a gasoline sample. The parameters of the calibration curve (slope and intercept) were
TE
calculated using linear regression. The estimated regression statistics showed a strong linear relationship between concentration and instrument response values and a regression line
CE P
that can reasonably be assumed to pass through zero. The slope and the intercept of the calculated calibration function were treated as a vector output quantity characterized by a joint probability distribution (bivariate Gaussian). Two
AC
types of coverage region were considered in order represent the vector quantity value with a stated probability (95%): a rectangle centered coverage region and an ellipse centered coverage region. Because of the mutual correlation between slope and intercept, the ellipse centered coverage region is more appropriate than the rectangular one.
The vector output quantity of the slope and intercept and the associated covariance matrix were then used as an input to the second stage of the measurement model which is the use of the calibration function for the estimation the quantity value (nitrogen content) corresponding to the unknown gasoline sample indication. The predicted value, its measurement uncertainty and a 95% coverage region were evaluated using MCM employing fixed number of trials (106). The algorithm of the MCM, implemented in MATLAB, gave a predicted value of 2.09 mg L-1, a standard measurement uncertainty of 0.04 mg L-1 and a
ACCEPTED MANUSCRIPT 95% symmetrical coverage interval [2.01 – 2.17] mg L-1. The results of the MCM agreed well
IP
T
with the outcome of the equation giving the standard error of the estimate.
SC R
References
[1] JCGM 100:2008, Guide to the expression of uncertainty in measurement (GUM 1995 with minor corrections), Joint Committee for Guides in Metrology (2009)
NU
[2] JCGM 101:2008, Evaluation of measurement data - Supplement 1 to the “Guide to the expression of uncertainty in measurement" - Propagation of distributions using a
MA
Monte Carlo method, Joint Committee for Guides in Metrology, 2008. [3] JCGM 102:2011, Evaluation of measurement data – Supplement 2 to the "Guide to the expression of uncertainty in measurement" – Extension to any number of output
D
quantities, Joint Committee for Guides in Metrology, 2011.
TE
[4] R. Bro, A. Rinnan, N.M. Faber, Standard error of prediction for multilinear PLS 2. Practical implementation in fluorescence spectroscopy, Chemometrics and
CE P
Intelligent Laboratory Systems 75 (2005) 69–76. [5] E. Saari, P. Perämäki, J. Jalonen, Measurement uncertainty in the determination of total petroleum hydrocarbons (TPH) in soil by GC-FID, Chemometrics and Intelligent
AC
Laboratory Systems 92 (2008) 3–12. [6] Carotenuto, G. Giovinco, B. Viglietti, L. Vanoli, A new procedure for the determination of calibration curves for a gas chromatograph used in natural gas analysis, Chemometrics and Intelligent Laboratory Systems 75 (2005) 209– 217
[7] D.B. Hibbert, The uncertainty of a result from a linear calibration, Analyst 131 (2006) 1273–1278. [8] V.J. Barwick, Preparation of calibration curves—a guide to best practice, VAM, Teddington UK, 2003. [9] L. Cuadros-Rodríguez, L. Gámiz-Gracia, E. Almansa-López, J. Laso-Sánchez Calibration in chemical measurement processes: I. A metrological approach, TrAC - Trends in Analytical Chemistry 20 (2001) 195–206. [10]D. Theodorou, L. Meligotsidou, S. Karavoltsos, A. Burnetas, M. Dassenakis, M. Scoullos, Comparison of ISO-GUM and Monte Carlo methods for the evaluation of
ACCEPTED MANUSCRIPT measurement uncertainty: Application to direct cadmium measurement in water by GFAAS, Talanta 83(2011) 1568-1574 [11]ASTM D4629-12, Standard Test Method for Trace Nitrogen in Liquid Petroleum
T
Hydrocarbons by Syringe/Inlet Oxidative Combustion and Chemiluminescence
IP
Detection, ASTM International, West Conshohocken, PA, 2012, www.astm.org [12]A. Possolo, Copulas for uncertainty analysis, Metrologia 47 (2010) 262-271.
SC R
[13]JCGM 104:2009. Evaluation of measurement data - An introduction to the “Guide to the expression of uncertainty in measurement" and related documents, Joint Committee for Guides in Metrology, 2009.
NU
[14]M.A. Herrador, A.G. Asuero, A.G. González, Estimation of the uncertainty of indirect measurements from the propagation of distributions by using the Monte-Carlo method: An overview, Chemometrics and Intelligent Laboratory Systems 79 (2005)
MA
115–122
[15]P.M. Harris, M.G. Cox, On a Monte Carlo method for measurement uncertainty evaluation and its implementation, Metrologia 51 (2014) S176-S182.
D
[16]W. Bich, M. G. Cox, P.M. Harris Evolution of the 'Guide to the Expression of
TE
Uncertainty in Measurement', Metrologia 43 (2006) S161–S166 [17]D. Theodorou, Y. Zannikou, G. Anastopoulos, F. Zannikos, Coverage interval
CE P
estimation of the measurement of Gross Heat of Combustion of fuel by bomb calorimetry: Comparison of ISO GUM and adaptive Monte Carlo Method, Thermochimica Acta 526 (2011) 122– 129
AC
[18]M. Solaguren-Beascoa Fernández, J.M. Alegre Calderón, P.M. Bravo Díez, Implementation in MATLAB of the adaptive Monte Carlo method for the evaluation of measurement uncertainties, Accreditation and Quality Assurance 14 (2009) 95– 106. [19]MATLAB® and Statistics Toolbox Release 2012b, The MathWorks, Inc., Natick, Massachusetts, United States.
ACCEPTED MANUSCRIPT Quantity values, xi, provided by measurement (calibration) standards and corresponding instrument indications (‘signals’, yi)
1st stage of the measurement model Measurement model with multiple output quantities
IP
T
Estimation of the parameters (e.g slope and intercept) of the calibration function (relation between the indication and the quantity value)
SC R
2nd stage of the measurement model
Analysis of an unknown sample, indication y0
NU
Estimation of the quantity value, xpred, corresponding to the unknown sample, indication, y0
MA
Fig. 1 Two-stage measurement model involving the estimation of the parameters of a calibration function (1st stage of the measurement model, having multiple output quantities) and the use of a calibration function for the estimation of a quantity value (2nd stage of the
AC
CE P
TE
D
measurement model)
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Fig. 2 Calibration graph describing relationship between nitrogen content and instrument
AC
CE P
TE
D
response
ACCEPTED MANUSCRIPT
NU
SC R
IP
T
[20]
AC
CE P
TE
D
MA
Fig. 3 Marginal distributions of the slope (b1) and the intercept (b0)
ACCEPTED MANUSCRIPT
T
5450
IP SC R
slope (L/mg)
5400
5350
NU
5300
-200
-100
0
100
MA
5250
200 300 intercept (AU)
400
500
600
700
800
Fig. 4 Elliptical and rectangular coverage regions (coverage probability p=0.95) of the joint
AC
CE P
TE
D
probability density function of the slope, b1 and the intercept,b0.
ACCEPTED MANUSCRIPT
Data p=0.500 p=0.750 p=0.900 p=0.950 p=0.990 p=0.999
T
5450
IP
5350
SC R
slope (L/mg)
5400
5300
-200
0
200 400 intercept (AU)
MA
-400
NU
5250
600
800
1000
Fig. 5 Elliptical coverage regions (for various coverage probabilities p) of the joint probability
AC
CE P
TE
D
density function of the slope, b1 and the intercept,b0.
ACCEPTED MANUSCRIPT 8
T
7
IP
6
SC R
Probability density
5
4
NU
3
1
0 1,90
2,09 nitrogen content (mg/L)
2,17
2,25
TE
D
2,01
MA
2
Fig. 6 Probability density function of xpred provided by Monte Carlo Method. The vertical
AC
CE P
dashed lines indicate the 95% probabilistically symmetrical coverage interval.
ACCEPTED MANUSCRIPT Table 1 Instrument parameters used for nitrogen determination Value
Volume injected (μL)
40
Syringe drive rate (μL s-1)
1
IP
T
Parameter
Furnace temperature (oC) Inlet oxygen flowmeter setting (mL min-1)
SC R
1050
AC
CE P
TE
D
MA
NU
Inlet carrier (Argon) flowmeter setting (mL min-1)
100 100
ACCEPTED MANUSCRIPT Table 2 Calibration and sample measerement raw data Calibration standard response (AU) replicate
replicate replicate
(mg L-1)
1
2
0.0
557.8
546.5
0.1
966.5
0.5
3057
1.0
5304
5353
5324
5.0
25890
26490
27850
53740
54840
53310
CE P
TE
IP 447.8
SC R
Unknown sample response (AU)
AC
3
895.8
927.6
2963
2868
NU
MA D
10.0
T
Concentration
11440
ACCEPTED MANUSCRIPT Table 3 Results of least square linear regression Regression parameter
Value
Estimate of the intercept, b0
263.99 AU 5352.99 AU mg-1 L
T
Estimate of the slope, b1
149.94 AU
IP
Standard error of the estimate of the intercept, SE(b0)
32.69 AU mg-1 L
SC R
Standard error of the estimate of the slope, SE(b1) Coefficient of determination, r2 standard
error
of
regression
(residual standard deviation), SEregression
value (p-value) Degrees of freedom, df
D
Regression sum of squares, SSregression
MA
The F statistic or Fisher – Snedecor F-
AC
CE P
TE
The residual sum of squares, SSresidual
507.4 AU
NU
The
0.9994
26821 (3.1 10 -27)
16 6.91 109 AU2 4.12 106 AU2
ACCEPTED MANUSCRIPT Table 4 Implementation program of Monte Carlo Method in MATLAB®
Line
Matlab Code function [x_pred]=regression_chem_compact(M)
2
% Calibration Data
3
x1=0; x2=0.1; x3=0.5; x4=1; x5=5; x6=10; x7=0; x8=0.1; x9=0.5; x10=1; x11=5; x12=10; x13=0; x14=0.1; x15=0.5; x16=1; x17=5; x18=10;
4
y1=557.8; y2=966.5; y3=3057.0; y4=5304; y5=25890; y6=53740; y7=546.5; y8=895.8; y9=2963; y10=5353; y11=26490; y12=54840; y13=447.8; y14=927.6; y15=2868; y16=5324; y17=27850; y18=53310;
5
%mean of all concentrations
6
x_av=(x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+x15+x16+x17+x18)/18;
7
% mean of all responses
8
y_av=(y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y11+y12+y13+y14+y15+y16+y17+y18)/18;
9
% Obtain estimates for the slope and the intercept
10
Xi=[x1; x2; x3; x4; x5; x6; x7; x8; x9; x10; x11; x12; x13; x14; x15; x16; x17; x18];
11
Yi=[y1; y2; y3; y4; y5; y6; y7; y8; y9; y10; y11; y12; y13; y14; y15; y16; y17; y18];
12
b1=sum((Xi-x_av).*(Yi-y_av))./sum((Xi-x_av).^2)
13
b0=y_av-b1.*x_av
14
%standard errors of slope and intercept
15
u_b1=32.69;
16
u_b0=149.94;
17
% simulate b1,b0 from joint pdf
18
r_b1b0=-sum(Xi)/sqrt(sumsqr(Xi)*numel(Xi)); % correlation coefficient
19
u_b1b0=u_b1*u_b0*r_b1b0;
20
%covariance matrix
21
V=[u_b1^2 u_b1b0; u_b1b0 u_b0^2];
22
%mean matrix
23
E=[b1; b0];
24
X=mvnrnd(E,V,M); % Mx2 matrix
25
B1=X(:,1); % 1st column of Mx2 matrix
26
B0=X(:,2); % 2nd column of Mx2 matrix
27
% simulate sample of y0
28
[a,b]=unif_param(11440,198.1466124);
AC
CE P
TE
D
MA
NU
SC R
IP
T
1
ACCEPTED MANUSCRIPT Line
Matlab Code y0=unifrnd(a,b,M,1);
30
% obtain x_pred
31
x_pred=(y0-B0)./B1;
32
% Results
33
x_pred=sort(x_pred);
34
Mean_Value=mean(x_pred)
35
Standard_Uncertainty=std(x_pred)
36
figure;
37
lower=min(x_pred);
38
upper=max(x_pred);
39
xc=lower:(upper-lower)/499:upper;
40
n=hist(x_pred,xc);
41
bar(xc,n./(((upper-lower)/499)*M),1);
42
hold on;
43
%ESTIMATE AND PLOT COVERAGE REGIONS FOR b1, b0
44
data=[B0 B1];
45
figure; % 95% rectangular and elliptical coverage region
46
%Plot the original data
47
plot(data(:,1), data(:,2), 'k.');
48
xlim([min(B0)-50, max(B0)+50]);
49
ylim([min(B1)-50, max(B1)+50]);
50
hold on;
51
% Set the axis labels
52
hXLabel = xlabel('b0');
53
hYLabel = ylabel('b1');
54
% Plot 95% rectangular coverage region
55
rectangle('Position',[b0-2.24*u_b0, b1-2.24*u_b1, 2*2.24*u_b0, 2*2.24*u_b1], 'Linestyle', '--')
56
% Plot 95% elliptical coverage region
57
% Calculate the eigenvectors and eigenvalues
58
covariance = cov(data);
59
[eigenvec, eigenval ] = eig(covariance);
60
% Get the index of the largest eigenvector
AC
CE P
TE
D
MA
NU
SC R
IP
T
29
ACCEPTED MANUSCRIPT Line
Matlab Code [largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
62
largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);
63
% Get the largest eigenvalue
64
largest_eigenval = max(max(eigenval));
65
% Get the smallest eigenvector and eigenvalue
66
if(largest_eigenvec_ind_c == 1)
SC R
67
smallest_eigenval = max(eigenval(:,2));
68
smallest_eigenvec = eigenvec(:,2); else
NU
69
IP
T
61
smallest_eigenval = max(eigenval(:,1));
71
smallest_eigenvec = eigenvec(1,:);
MA
70
end
73
% Calculate the angle between the x-axis and the largest eigenvector
74
angle = atan2(largest_eigenvec(2), largest_eigenvec(1));
75
% This angle is between -pi and pi.
76
% Shift it such that the angle is between 0 and 2pi
77
if(angle < 0)
angle = angle + 2*pi;
CE P
78
TE
D
72
end
80
% Get the coordinates of the data mean
81
avg = mean(data);
82
% Get the 95% confidence interval error ellipse
83
chisquare_val = 2.4477; % set chisquare_val = 1.1774 for 50%, 1.6651 for 75%, 2.1460 for 90%, 3.0349 for 99%, 3.7169 for 99.9% confidence levels
84
theta_grid = linspace(0,2*pi);
85
phi = angle;
86
X0=avg(1);
87
Y0=avg(2);
88
a=chisquare_val*sqrt(largest_eigenval);
89
b=chisquare_val*sqrt(smallest_eigenval);
90
% the ellipse in x and y coordinates
91
ellipse_x_r
= a*cos( theta_grid );
92
ellipse_y_r
= b*sin( theta_grid );
AC
79
ACCEPTED MANUSCRIPT Line
Matlab Code %Define a rotation matrix
94
R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];
95
%Rotate the ellipse to some angle phi
96
r_ellipse = [ellipse_x_r;ellipse_y_r]' * R;
97
% Draw the error ellipse
98
plot(r_ellipse(:,1) + X0,r_ellipse(:,2) + Y0,'k-')
99
figure % marginal probabilities of b0, b1
SC R
IP
T
93
100 subplot(1,2,1);
NU
101 B1=sort(B1); 102 Mean_value=mean(B1);
MA
103 lower=min(B1); 104 upper=max(B1);
105 xc=lower:(upper-lower)/199:upper;
D
106 n=hist(B1,xc);
108 title('b1')
CE P
109 subplot(1,2,2);
TE
107 bar(xc,n./(((upper-lower)/199)*M),1);
110 B0=sort(B0);
111 Mean_value=mean(B0); 112 lower=min(B0);
AC
113 upper=max(B0);
114 xc=lower:(upper-lower)/199:upper; 115 n=hist(B0,xc);
116 bar(xc,n./(((upper-lower)/199)*M),1); 117 title('b0') 118 end 119 function [a,b]=unif_param(m,v) 120 % A function to obtain the parameters of the uniform (rectangular) 121 % distribution given its mean and standard deviation 122 b=m+v*sqrt(12)/2; 123 a=2*m-b; 124 end
ACCEPTED MANUSCRIPT Highlights •
The measurement uncertainty of the nitrogen content in gasoline determined using a calibration curve is evaluated This is a typical multistage measurement model with more than one output quantity
•
The coefficients of the calibration curve are characterized by a joint probability
IP
T
•
Monte Carlo Method is used to evaluate the calibration curve uncertainty
CE P
TE
D
MA
NU
component
AC
•
SC R
distribution