Issue 
Int. J. Metrol. Qual. Eng.
Volume 13, 2022



Article Number  9  
Number of page(s)  15  
DOI  https://doi.org/10.1051/ijmqe/2022008  
Published online  02 August 2022 
Research Article
Determining the covariance matrix for a nonlinear implicit multivariate measurement equation uncertainty analysis
Department of Mechanical Engineering, University of South Africa, Private Bag X6, Florida 1710, South Africa
^{*} Corresponding author: ramnav@unisa.ac.za
Received:
27
February
2022
Accepted:
7
June
2022
The application of the Guide to the Expression of Uncertainty in Measurement (GUM) for multivariate measurand equations requires an expected vector value and a corresponding covariance matrix in order to accurately calculate measurement uncertainties for models that involve correlation effects. Typically in scientific metrology applications the covariance matrix is estimated from Monte Carlo numerical simulations with the assumption of a Gaussian joint probability density function, however this procedure is often times considered too complex or cumbersome for many practicing metrologists in industrial metrology calibration laboratories, and as a result a problem which occurs is that correlation effects are frequently omitted so that uncertainties are approximated through a simple rootsumsquare of uncertainties which leads to inaccuracies of measurement uncertainties. In this paper, a general purpose deterministic approach is developed using a computer algebra system (CAS) approach that avoids the need for Monte Carlo simulations in order to analytically construct the covariance matrix for arbitrary nonlinear implicit multivariate measurement models. An illustrative example for a multivariate SakumaHattori pyrometer equation with the proposed method is demonstrated with explanations of underlying Python code.
Key words: Covariance matrix / Guide to the Expression of Uncertainty in Measurement (GUM) / MonteCarlo simulation (MCS) / multivariate measurement uncertainty / pyrometry
© V. Ramnath, Published by EDP Sciences, 2022
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
1.1 Research motivation
The original Guide to the Expression of Uncertainty in Measurement that is commonly simply abbreviated as the GUM [1] has largely replaced the earlier Kline & McClintock uncertainty analysis technique [2] that has historically been utilized in many engineering research and applications work. It achieves this by focusing on explicit univariate measurand models of the form Y = f (X _{1}, … , X _{ N }), or equivalently Y = f ( X ) with X = [X _{1}, … , X _{ N }]^{T} ∈ ℝ ^{ N }, where X _{1}, … , X _{ N } are known scalar inputs with corresponding uncertainty information inclusive of possible covariances Cov(X _{ i } X _{ j }) for i ≠ j and i, j = 1, … , N, and where Y is a single scalar output with an unknown uncertainty. This uncertainty method is utilized in order to determine the combined standard uncertainty u _{ c }(Y) and corresponding expanded uncertainty U(Y) = k _{ p } ⋅ u(Y) for Y under certain limiting assumptions where k _{ p } is a suitable coverage factor for a specified confidence level p.
Initially the validity assumptions associated with theGUM were mainly restricted to three conditions, namely (i) the necessity for specifying a linearisation of f ( X ) in a small local neighbourhood of X ∈ ℝ ^{ N } that would specify a domain for which u(Y) is valid, (ii) the assumption of Gaussian uncertainties u(X _{ i }) for the inputs such that where μ _{ i } and are corresponding equivalent expected and variance values for X _{ i } for i = 1, … , N using standard statistical transformations to convert a nonGaussian probability density function (PDF) such as rectangular or triangular PDFs to equivalent Gaussian PDFs and where x _{ i } is a random variable corresponding to X _{ i }, and (iii) the validity of the WelchSatterthwaite formula for estimating an effective degreesoffreedom v _{ eff} in order to obtain a suitable coverage factor k _{ p } to calculate expanded uncertainties and corresponding confidence intervals for specified confidence levels in cases whereacorrelation in the inputs were present such that Cov(X _{ i }, X _{ j }) ≠ 0 for i, j = 1, … , N.
For the original GUM method the PDF for Y is g _{ Y }(y) where y is a random variable associated with Y such that y ∼ g _{ Y }(y) and this PDF is obtained by using the WelchSatterthwaite formula such that g _{ Y }(y) is approximated as a Student's tdistribution using v _{eff}. As a result, the key simplicity for performing uncertainty calculations for u(Y) with an explicit univariate model Y = f ( X ) is that the uncertainty may in principle be fully accomplished analytically without any need for advanced numerical simulations.
Following the adoption of the GUM a GUM Supplement 1 [3] was introduced that extended the original validity assumptions to models for (a) implicit univariate equations of the form h(Y, X ) = 0, and (b) models in which the input PDFs did not follow a normal distribution. With this new development general and possibly nonGaussian joint PDFs for the input PDFs x ∼ g _{ x }(ξ) and output PDF y ∼ g _{ Y }(η) were now possible to model correlations in the inputs and nonGaussian behaviour for both the inputs and outputs. A key achievement of the GUM Supplement 1 is that it officially introduced the metrology community to the use of theMonte Carlo method as an advanced uncertainty analysis technique where as previous earlier uncertainty analysis work was almost wholly analytical based. This then allowed metrologists to investigate underlying measurement systems without any limiting validity assumptions. Consequently newer theoretical techniques for reporting univariate measurement uncertainties that were nonGaussian were then developed by Willink [4] using quantile functions, and this was later applied by Ramnath [5] in practical engineering applications.
Multivariate uncertainty analysis frameworks were then subsequently published in the GUM Supplement 2 [6] which further extended the underlying GUM uncertainty framework approach for models for both explicit equations as well as implicit equations, and again allowed for the possibility of nonGaussian PDFs for both inputs as well as outputs.
An explicit multivariate equation of the form(1)has an known input X = [X _{1}, … , X _{ N }]^{T} of nominal values with an associated joint PDF g _{ X }(ξ) to model the input uncertainty that is introduced into the model in order to work out the output Y = [y _{1}, … , y _{ m }]^{T} and its associated joint PDF g _{ Y }(η). Although the output uncertainty may or may not be a multivariate Gaussian PDF within the GUM framework the unknown output uncertainty is specifically modelled as a multivariate Gaussian PDF for which an output covariance matrix U _{ Y } ensures mathematical closure. The output uncertainty u( Y ) is formally defined in terms of the matrix equation(2)
For the above equations the matrices are specified as(3) (4) (5)
Referring to the above system of equations it may be seen that the sensitivity matrix C _{ x } may be conveniently and explicitly obtained by calculating the partial derivatives for i = 1, … , m and for j = 1, … , N respectively and evaluating the expressions at X = x for a specified value of the random variable x , and that the covariance matrix for Y i.e. the uncertainty u( Y ) may then be directly calculated through an explicit matrix equation by a simple matrix multiplication once U _{ x } and C _{ x } are determined in order to give the corresponding uncertainty for y that satisfies Y = f ( X ). As an example if calculations are performed in Matlab or Gnu Octave as discussed by Hansen [7] the matrix equation for U _{ y } may be simply calculated as a one line of code as
As a result it may be observed that in the case of a general multivariate explicit model Y = f ( X ) that the uncertainty analysis approach is conceptually straightforward, and that there is no need for any further advanced numerical simulation in order to obtain the covariance matrix U _{ y }.
Special techniques for determining the covariance matrix of a linear multivariate explicit model in the special case of a straight line y = αx + b with inputs x _{ i } and outputs y _{ i } for datapoints i = 1, 2, … , N modelled as a multivariate regression equation in order to determine the uncertainties of the parameters has been reported earlier by Ramnath [8]. In this work the research objective was to calculate the uncertainty of y at a specified value of x with the formula .
Considering on the other hand an implicit multivariate measurement equation, with an illustrative example in the Appendix to outline and explain the mathematical formulation, of the form(6)the uncertainty u( Y ) may now be obtained from the covariance matrix U _{ y } that is specified from the matrix equation(7)
In the above equation the input sensitivity matrix C _{ x } is of dimension m × N with matrix elements C _{ ij } defined as where i = 1, … , m and j = 1, … , N, whilst the output sensitivity matrix C _{ y } is of dimension m × m with matrix elements C _{ℓk } defined as where l = 1, … , m and k = 1, … , m, such that(8) (9)
Referring to the above uncertainty analysis equation it may be observed that there is no simple or direct analytical approach to solve for the unknown covariance matrix U _{ y }, and this observation provides the research motivation to investigate and develop a method that may be used to analytically construct the covariance matrix without recourse to more complex stochastic numerical techniques such as Monte Carlo simulations which are computationally expensive and require specialist techniques to postprocess as reported by Ramnath [9]. A key simplification that is achieved if a covariance matrix is calculated is that a multivariate measurement uncertainty is completely encapsulated within this matrix and the distribution function may simply be taken as a multivariate Gaussian distribution.
In scientific metrology applications the Monte Carlo method as a stochastic technique is generally preferred to solve for the output covariance matrix U _{ Y } due to the lack of a convenient and quick method for the solution of the matrix equation . This is formally achieved as an application of the GUM Supplement 2 by solving the Markov formula that specifies a convolution integral of the form(10)
In the above equation δ(⋅) is the Dirac delta function as discussed by Cox & Siebert [10] and the Markov formula may be used in the most general case in order to numerically estimate U _{ y } by postprocessing the Monte Carlo data Ω _{ Y }.
Implementing this theoretical stochastic approach in scientific metrology applications such as applied engineering problems as reported by Ramnath [11] whilst feasible, is nevertheless an approach that demands the use of advanced software engineering considerations for sequentially estimating the covariance matrix and confidence region as discussed Harris and Cox [12], and is thus frequently too complex an undertaking for many metrologists within industry that do not have access to the requisite specialist mathematical and statistical training.
Due to these technical observations in the absence of full Monte Carlo simulations, many metrologists often do not correctly quantitatively account for input and output correlation effects in multivariate implicit model uncertainty analysis problems.
1.2 Research objective
Motivated by the above problem formulation the research objective in this paper is to develop an approach to solve for the unknown covariance matrix U _{ Y } in a nonlinear implicit multivariate measurement equation h ( Y , X ) = 0 by analytical calculations to reduce and minimize the need for unnecessarily complicated numerical techniques that may not be accessible to metrology practitioners working in calibration laboratories within the industrial metrology field.
In this paper, to accomplish this research goal a general purpose algebraic approach is investigated in order to develop a method that avoids Monte Carlo simulations by instead directly analytically constructing the output covariance matrix for arbitrary nonlinear implicit multivariate measurement models. Results of the proposed method will be utilized to validate and verify the proposed method's feasibility by using a threedimensional vector measurement equation model based on a SakumaHattori pyrometry equation where the temperaturesignal has the form(11)
The objective of this paper is to determine the covariance matrix for the parameters A, B and C in the equation(12)
When the above covariance matrix is determined it may then conveniently be used in the corresponding inverse signal equation for the temperature(13)to conveniently calculate the temperature uncertainty u(T) when the uncertainty in the signal u(S) is specified.
The practical utility of determining a calibrated pyrometer's SakumaHattori equation parameter expected values and covariances is that many national measurement laboratories that maintain their respective country's national high temperature radiation thermometry scale from 600°C to about 2000°C use this particular equation to completely calibrate and characterize their respective pyrometer standards. When such a pyrometer laboratory standard is calibrated by determining the underlying covariance matrix for the instrument all the relevant uncertainty information is then fully quantified, so that the pyrometer may then be used to uniquely determine measured high temperatures at the highest possible accuracy levels.
In practical terms such a pyrometer standard when characterized with a covariance matrix may then be used to measure the true temperature T _{ BB } of a Variable Temperature High Temperature Blackbody (VTHTBB). The resulting blackbody temperature source may then be used as a transfer medium to in turn calibrate a client pyrometer from industry to give a Unit Under Test (UUT) measured value T _{ UUT }. The true temperature of the VTHTBB is subtracted from the measured client pyrometer temperature to give a UUT error e _{ UUT } = T _{ UUT } − T _{ BB } which may then be used to compensate the UUT measurements to produce a true temperature. Consequently, in this measurement traceability scheme every client pyrometer within a particular country that is used in industries such as steel fabrication, manufacturing, energy production and materials processing amongst others that is calibrated with pyrometry measurements traceable to a particular country's national measurement high temperature scale, is ultimately dependent on the determined covariance matrix of that country's pyrometer standard.
This particular traceability scheme is also applicable in other measurement instruments such as pressure balances, mass transducers, and flow rate meters amongst other physical laboratory instruments in national measurement laboratories and commercial industrial calibration laboratories. The result is that all measurement instrument results and accuracies within industry and the engineering sector are ultimately characterized by a national measurement laboratory's specific measurement instrument covariance matrix.
2 Literature review
2.1 Overview of existing statistical approaches
Earlier work by Warsza and Puchalski [13] to study the effect of correlations in multivariate measurement models only considered an explicit equation of the form Y = F(X) where X = [X _{1}, X _{2}, … , X _{ N }]^{T} and Y = [Y _{1}, Y _{2}, … , Y _{ m }]^{T}. An example of a multivariate measurement equation is the case of a multiparameter digital multimeter which is a single instrument that can measure voltage V _{ meas } = Y _{1}, current I _{ meas } = Y _{2} and resistance R _{ meas } = Y _{3} based on input standards of resistance R _{ source } = X _{1} and voltage V _{ source } = X _{2}. The inputs X _{1} and X _{2} by necessity include both nominal expected values of these quantities as well as any correlation Cov(X _{1}, X _{2}) in these inputs for statistical completeness. When implemented in a model Y = F(X) the underlying input uncertainty U _{ X } will be logically carried through into the measurement mathematical model to produce both the nominal expected values Y _{1}, Y _{2} and Y _{3} as well as their corresponding uncertainties. Covariances conveniently model the output uncertainties and are defined by the terms Cov(Y _{1}, Y _{2}), Cov(Y _{1}, Y _{3}) and Cov(Y _{2}, Y _{3}) on the offdiagonal elements of U _{ Y }, whilst the variances Var(Y _{1}, Y _{1}), Var(Y _{2}, Y _{2}) and Var(Y _{3}, Y _{3}) are located on the diagonal elements of the U _{ Y } matrix. In this work the multivariate uncertainty of the output was approximated as where the elements of the sensitivity matrix of dimension m × N propagate the output uncertainty as U _{ Y } = S ⋅ U _{ X } ⋅ S ^{ T } and the input and output matrices are
The key observation is that in the special case of an explicit multivariate model that the sensitivity matrix S may be analytically determined through simple partial derivatives, and if the output Y is further processed through an equation Z = G(Y) then the covariance will again be carried through such that the covariance matrix for Z is then . The disadvantage of this approach is that it cannot be readily applied to implicit multivariate measurement equations for which no algebraic solution exists and this presents a research gap in determining the covariance matrix.
A practical measurement traceability chain that illustrates how the effect of covariances are propagated is shown in Figure 1 where a national metrology institute (NMI) for a country has a set of apex measurement standards that are used to calibrate a national laboratory standard that is then in turn used to calibrate an instrument from industry and in this way the information in the covariance matrix is propagated to every instrument in industry within that country using a twostage process as later discussed by Forbes [14].
This measurement traceability scheme by international measurement conventions also applies in pyrometry calibration practice where a national laboratory typically has a set of three reference temperature measurements, such as fixed point reference blackbodies that provide known thermal radiance sources for the silver, copper and gold freezing temperatures, which are used to calibrate a national laboratory's pyrometer standard. This national laboratory pyrometry standard in turn is completely characterized by a measurement equation such as a SakumaHattori equation set of parameter values and parameter covariances that completely and uniquely defines the laboratory standards measurement performance. When the national laboratory pyrometer standard for a particular country is then used to calibrate client pyrometers from industry every radiation thermometry measurement in that country be it in the manufacturing, energy production or materials processing sectors of that country that uses high temperature measurements are ultimately traceable to and affected by the covariance matrix for the national laboratory standard which must be determined and quantified in order to provide credible and accurate measurement traceability.
As reported by van der Heen and Cox [15] sometimes national measurement laboratories may omit or ignore correlations because it is either too difficult to compute or for other reasons and this omission can lead to poor measurement decisions or logical absurdities when determining measurement results and equivalences such as in key comparison reference values (KCRVs) which are used to produce national laboratory standards.
The above scientific metrology approaches are fundamentally based on the GUM and GUM supplements, and although they may appear similar to uncertainty analysis in other areas of physics and engineering certain differences of approach are present. As an example in the area of nuclear physics as recently discussed by Kornilov et al. [16] the covariance matrix is generally produced through simulations by not considering the underlying measurement uncertainty through a χ ^{2} merit function goodness of fit as originally communicated by Press et al. [17] which thus differs from the recommended approach of the GUM Supplement 2 for multivariate uncertainties. When discrepancies arise which cannot be explained by the uncertainty of the measurement data then additional Monte Carlo simulations are performed in an attempt to generate a Systematic Distortion Function (SDF) which may be thought of as a “correction function” to eliminate systematic bias in the data. In the case where a SDF cannot adequately account for underlying discrepancies then an additional Unrecognized Source of Uncertainties (USU) functional is generated which links to the SDF in an attempt to modify the underlying model to more closely align to the data even though this may render the model as deviating from relevant physical principles or laws. Alternatives to SDF/USU schemes also include the method of determining a covariance matrix by analysing the noise in a model as discussed by Chhabra et al. [18] i.e. where a signal y is composed as a superposition of an input x and a separate noise term ε so that y = x + ϵ. This approach is generally inconsistent with the conventional metrology approach which considers the probability density function of a signal y as a random variable g(y) which continuously varies as analysed from a wholly Bayesian statistics framework where all the sources of uncertainty are statistical aleatoric uncertainties, and which does not have a known constant value or is composed of separate systematic epistemic uncertainties and random aleatoric uncertainties. These alternative uncertainty approaches in other fields of study some of which mix aleatoric and epistemic uncertainties in “grey box models” as discussed by Brastein et al. [19] in some branches of engineering and physics, are therefore considered to be inconsistent with the accepted guidelines in scientific metrology work where any mathematical technique must be underpinned by an appropriate physical law. A fundamental difference in approach particularly for scientific metrology and industrial metrology calibration work is therefore that all uncertainty information must be fully encapsulated completely within the measurement equation through reporting of a covariance matrix as discussed by Smith et al. [20]. Consequently, data analysis postprocessing approaches such as SDF/USU functions amongst other techniques are considered inappropriate and inconsistent for national metrology institutes that maintain and disseminate national measurement standards for various countries and their industries.
Fig. 1 Illustration of measurement traceability scheme showing how information from a covariance matrix of a national laboratory standard is transmitted to a client instrument measurement standard in industry to propagate uncertainty analysis results. 
2.2 Pyrometry theory and measurement modelsc
Saunders [21] gives a general expression for the signal S measured by a pyrometer as(14)where s(λ) is the pyrometer spectral responsivity for wavelengths in the range λ _{min} λ λ _{max}, K is an instrument calibration constant that includes geometrical, optical and electrical factors for the pyrometer, and L _{ b }(λ, T) is the blackbody spectral radiance using Planck's law such that(15)
For the above equation c _{1} and c _{2} are the first and second radiation constants, λ is the wavelength of the optical radiation in the medium in which the pyrometer is immersed which may be assumed to be air, and n is the refractive index of the medium. Using the above Planck version for the spectral radiance a general form of the Sakumahattori equation that is suitable for both narrow band and wide band wavelengths may take the form as previously discussed along with the corresponding inverse temperature equation of the SakumaHattori equation. If the covariance matrix is known then the uncertainty in temperature from the pyrometer signal may then be calculated as(16)
In general, for wide band and broad band spectral responsivities there are no explicit formulae to analytically calculate the parameter values for A, B, C as this is a nonlinear regression analysis and a trialanderror approach must be adopted to obtain suitable starting values for a subsequent LevenbergMarquardt optimization as discussed by Press et al. [17].
Nevertheless the instrument parameters in the special case narrow band spectral responsivities where the relative bandwidth of the spectral responsivity is small as discussed by Saunders & White [22] the parameters A, B, C may then be approximated as(17) (18) (19)
In the above equation λ _{0} is the mean wavelength for λ _{min} λ λ _{max} and σ is the standard deviation of the spectral responsivity using standard statistical formulae.
A commonly accepted experimental approach to determine the spectral responsivity of pyrometers in high temperature measurements was reported by Briaudeau et al. [23] for both absolute spectral responsivity R(λ) as well as relative spectral responsivity s(λ) characterizations of pyrometers. In this work the overall purpose is to cancel out any speckle noise caused by laser interferences inside an integrating sphere, so that a known absolute spectral radiance may be measured by the pyrometer.
Use is made of a laser as a tunable monochromatic quasiLambertian source which first produces a laser beam of a known constant spectral radiance at a wavelength λ, that then secondly passes through a multimode optical fibre housed inside an ultrasonic bath, and which then thirdly passes through an optical integrating sphere, before finally entering the pyrometer. The pyrometer which is used as a radiation thermometry device is composed of various filters and detector components that convert the input optical radiation in units of W ⋅ sr^{−1} ⋅ m^{−3} into a measured photocurrent signal in units of A. For this measurement system the pyrometer photocurrent I _{pyro}(λ, r = 0)/[A] is specified as(20)where the absolute spectral responsivity is is a nonlinearity correction factor for pyrometer responsivity, SSE is the sizeofsource effect for the pyrometer, and ε is constant blackbody emissivity.
In the special case the above calibration terms , SSE and ε may all be assumed to not have any spectral dependence for simplicity e.g. ϵ ≠ ϵ(λ) if a narrow band of radiation from a tunable laser is considered to be effectively monochromatic.
Later work by Yoon et al. [24] reported further technical details on the experimental aspects to consider for the measurement of spectral responsivities and extended the use of the SakumaHattori over a wider temperature from 400K to 1300K at NIST. These very high accuracy results were achieved by using a custom designed pyrometer with absolute spectral responsivities capable of absolute thermodynamic temperature measurements and used currenttovoltage amplifiers that were traceable to the quantum resistance and quantum voltage standards at NIST, and were able to achieve accuracies of around ±0.05 mK.
Manoi et al. [25] investigated radiation thermometry using two fixed points, which are typically taken as silver, copper or gold fixed points, corresponding to a n = 2 scheme for calibrating a temperaturesignal relationship for a pyrometer. In this approach for a pyrometer with a sufficiently narrow bandwidth it may be shown that equation parameters are independent of the shape of the spectral responsivity and may be quantified solely in terms of the mean wavelength λ _{0} of the spectral responsivity and standard deviation σ of the spectral responsivity. This approach typically produces errors smaller than ± 3 mK for temperature from 600°C to 3000°C when the relative bandwidth is . With this constraint for a 650 nm pyrometer a full width half maximum (FWHM) as shown in Figure 2 must be less than about 20 nm for the special case to work of the form(21)
When these special conditions apply the constants in a SakumaHattori equation may be obtained using c _{2} = 0.014388 m ⋅ K as C = c ^{′}, , and which are three unknowns to be determined. If the spectral responsivity is available and σ is estimated from the FWHM, then there are two remaining free parameters λ _{0} and c ^{′} and these can then be simultaneously determined from the pyrometer signals at the two known temperature fixed points using the (T _{1}, S _{1}) and (T _{2}, S _{2}) pairs of experimentally measured signal points where T _{1} and T _{2} are known fixed points by definition.
An earlier similar concept by Saunders [26] using two temperatures T _{0} and T _{1} and the temperaturesignal relationship , where R(λ) is the spectral responsivity that includes any constants due to geometrical and electrical parameters, allows for the temperature to be calculated to the ratio of the signal at a temperature relative to the signal at a reference temperature. Using this approach the integration may be removed by defining a temperaturedependent mean effective wavelength λ _{ m } so that the ratio of signals may be expressed as(22)
In the above equation the mean effective wavelength is calculated in terms of the limiting effective wavelength λ _{ T } that is a function of the single temperature which is calculated in terms of the known spectral responsivity such that(23)from which λ _{ m } is finally obtained as(24)where λ _{ T } is considered as an explicit function of the temperature.
Calibration of pyrometers using the concept of a mean effective wavelength λ _{ m } thus requires knowledge of the pyrometer detector linearity and spectral responsivity R(λ) which is a difficult and timeconsuming calibration that must be performed in a spectral radiometry laboratory with expensive specialist equipment.
Saunders [27] considered a univariate nonlinear regression for pairs of data (x _{1}, y _{1}), … , (X _{ N }, y _{ N }) where the model took the form(25)and obtained expressions with M data points in the absence of correlations to fit a a SakumaHattori equation where(26)where the sensitivity terms and were defined in terms of sums of products of matrices. A limitation of this earlier approach is that it is not easily generalizable to incorporate correlation effects in the parameters of the SakumaHattori equation and is more amenable to a Monte Carlo simulation.
Fig. 2 Illustration of the full width half maximum of a curve [Graphic Source: Arne Nordmann (norro), CC BYSA 3.0 http://creativecommons.org/licenses/bysa/3.0/, via Wikimedia Commons]. 
3 Mathematical modelling
In order to solve the matrix equation for the unknown output covariance matrix U _{ y } the relationship between the input sensitivity matrix C _{ x } and the output sensitivity C _{ y } as reported by the GUM Supplement 2 may be formally expressed as(27) (28)
Although there is an ostensible theoretical formula to compute the covariance matrix U _{ y } the key technical challenge is that the inverse matrix computation renders the formula as numerically unstable, and a further challenge is that a full symbolic calculation of the matrix inverse whilst in principle theoretically possible is not generally feasible on most desktop and laptop computers due to the excessive number of operations.
These challenges may be overcome by using a hybrid symbolic/numerical approach to first symbolically derive algebraic expressions for C _{ x } and C _{ y } and then evaluating their numerical values. Thereafter once the numerical values of the matrices are known a linear algebra decomposition algorithm provided for in the GUM Supplement 2 as summarized in Figure 3 may then be numerically implemented using readily available nonspecialist routines to calculate U _{ y }.
Considering the general form of the SakumaHattori equation and three pairs of experimental data points (T _{1}, S _{1}), (T _{2}, S _{2}) and (T _{3}, S _{3}) the equation to determine the parameters may be constructed as(29) (30) (31) (32)
The approach adopted in this paper is to use the sympy symbolic Python package reported by Meurer et al. [28] to utilize the computer algebra system to automatically calculate the sensitivity matrices.
When the above symbolic expressions are processed in sympy the following meritfunction expression results:(33)
The system of three nonlinear equations is then obtained by partial differentiation as:(34) (35) (36)
The input sensitivity matrix C _{ X } with elements specified by the partial derivatives for i = 1, 2, 3 and j = 1, 2, 3, 4, 5, 6 and similarly the output sensitivity matrix C _{ Y } with elements specified by the partial derivatives for i = 1, 2, 3 and j = 1, 2, 3 are omitted as these symbolic matrices are too complex and unwieldy to print. Nevertheless these symbolic matrices may easily be constructed in sympy and numerically evaluated using the numerical values of X _{1} = T _{1}, X _{2} = T _{2}, X _{3} = T _{3}, X _{4} = S _{1}, X _{5} = S _{2}, X _{6} = S _{3} and Y _{1} = A, Y _{2} = B, Y _{3} = C where the numerical values of Y _{1}, Y _{2} and Y _{3} are obtained from the numerical solution of h _{1}( Y , X ) = 0, h _{2}( Y , X ) = 0, and h _{1}( Y , X ) = 0.
Fig. 3 Numerically stable linear algebra algorithm for determining covariance matrix. 
4 Numerical simulations
Using the accepted International Temperature Scale of 1990 (ITS90) temperature data as reported by PrestonThomas [29] let the silver, gold and copper fixed point freezing temperatures be T _{1} = T _{ Ag } = 1234.93 K, T _{2} = T _{ Au } = 1337.33 K and T _{3} = T _{ Cu } = 1357.77 K and assume medium accuracy standard uncertainties of u(T _{1}) = ± 0.3 K, u(T _{2}) = ± 0.3 K and u(T _{3}) = ± 0.3 K based on representative fixed point blackbody furnace sources as reported in Sakuma & Hattori [30]. Then using the medium accuracy representative pyrometer relative spectral responsivity data by Saunders [21] as shown in Figure 4 and the NIST CODATA 2018 recommended values for physical constants available at https://physics.nist.gov/cuu/Constants/ set the first and second radiation constants to c _{1} = 3.741771852 × 10^{−16} W ⋅ m^{2} and c _{2} = 1.438776877 × 10^{−2} m ⋅ K.
From this data the pyrometer signal values may be calculated using the formula to yield the following underlying temperaturesignal data to fit the parameters in a SakumaHattori equation such that(37)and(38)
For simplicity assume that the standard uncertainty of the signals are all relatively large at ±2% of the nominal signal values so that u(S _{1}) = ± 1.191152160860418×10^{2} u(S _{2}) = ± 4.693939769213445×10^{2} and u(S _{3}) = ± 6.020984643486366×10^{2}.
Setting the input as X = [T _{1}, T _{2}, T _{3}, S _{1}, S _{2}, S _{3}]^{T} and the output as Y = [A, B, C]^{T} it immediately follows that N = 6 and m = 3 according to the previous notation. As a result the input sensitivity matrix C _{ x } will be of dimension 3 × 6 and the output sensitivity matrix C _{ y } will be of dimension 3 × 3. The sympy package in Python may conveniently be installed and then imported into Python to first symbolically construct the matrices C _{ x } and C _{ y } and then evaluate their numerical values as shown in Figure 5.
The final step is to solve the matrix equation using the algorithm in Figure 3 and this may conveniently be achieved with the Python code in Figure 6.
In order to check the validity of the results from the proposed method the supplied specified signaltemperature values (S _{1}, T _{1}), (S _{2}, T _{2}), (S _{3}, T _{3}) are compared to the corresponding predictions for the temperatures as obtained from the SakumaHattori equation using the inverse temperature equation and the optimized fitted parameter values as shown in Figure 7. Ideally if there was a perfect parity the predicted points would lie exactly on the diagonal line, however it is observed that this is not the case due to the imperfection of the pyrometer's spectral responsivity which has a finite spectral bandpass and is not a Diracdelta point function. In general pyrometry practice there is usually a balance or tradeoff that is involved with a selection of a spectral responsivity bandpass. A narrower bandpass in for example a filter radiometer will result in increased accuracies but over a smaller temperature range, whilst a broader bandpass will result in a greater sensitivity over a broader temperature range but at the expense of a lower temperature accuracy. Nevertheless when the expanded temperature uncertainties as calculated from equation (16) that incorporates the covariance information it is seen that the error bars are within limits of the parity thus indicating experimental consistency in uncertainty analysis predictions. These results may also be directly visualized with the signaltemperature curves in Figure 8 from which it is observed that for the same signal values that the predicted temperature values taking into account the uncertainties from the covariances are reasonably consistent. It may be observed that the slight inaccuracies are due to a combination of the the minimum number of points with n = 3 that are used for the nonlinear regression and the “leverage effect” where a single point is located “far away” from the main cluster of points. These effects in radiation thermometry studies using high accuracy pyrometers are largely mitigated by including a larger number of fixed points in the regression scheme as discussed by Yoon et al. [24] who used six fixed points, Wooliams et al. [31] who used fourteen fixed points, and with more recent work by Todd et al. [32] that seeks to extend and refine the number of high temperature fixed points.
Consequently once the covariance matrix U _{ Y } is determined for a national measurement laboratory's pyrometry standard it may in turn be used to perform high accuracy radiation thermometry measurements of VTHTBB's and in turn calibrate client pyrometers from industry, using a similar measurement traceability scheme as previously illustrated in Figure 1.
Fig. 5 Python code for symbolically and numerically constructing the matrices C _{ x }, C _{ y } and U _{ x }. 
Fig. 6 Python code for linear algebra implementation to postprocess C _{ x }, C _{ y } and U _{ x } to calculate U _{ y }. 
Fig. 7 Parity line curve to validate the accuracy of the model fit for the SakumaHattori equation. 
Fig. 8 Comparison of signaltemperature results from specified and predicted data for the SakumaHattori equation. 
5 Conclusions
Based on the research reported in this paper the following conclusions were determined:

A new method using a computer algebra system has been developed that demonstrated the utility of symbolically directly calculating the input and output sensitivity matrices for an arbitrary multivariate implicit measurement equation for a real practical problem in metrology.

An implementation method using standard linear algebra routines has been developed that allows for the output covariance matrix to be easily calculated without the need for lengthy and complex Monte Carlo simulations.
6 Influences and implications
Based on the research reported in this paper the following influences and implications are:

Metrologists working in industry are no longer forced to attempt to perform advanced Monte Carlo simulations for a category of uncertainty problems which may be symbolically solved using computer algebra systems.

Computer programs written in the open source language Python are now freely and publicly available for metrologists to freely modify and adapt to assist them in uncertainty analysis problems in other measurement fields.
Acknowledgments
This work was performed with funds provided by the Department of Higher Education, Science and Technology (DHEST) on behalf of the South African government for research by public universities.
Appendix A. Mathematical formulation of h(Y, X) = 0
In order to illustrate and explain the mathematical formalism of an implicit multivariate equation h ( Y , X ) = 0 consider a pressure balance experiment that is performed with a deadweight tester in which a mass m _{ i } is stacked on top of a piston within a cylinder that generates a known weight force in a downwards direction that is counterbalanced by a hydrostatic pressure P _{ i } that pushes the piston in an upwards direction. By a simple freebody diagram and equilibrium of forces it follows that the governing equation is(A.1)
In the above equation let i denote a particular index for an applied pressure, say i = 5 so that the sequence of applied pressures are P _{1}, P _{2}, P _{3}, P _{4}, P _{5} which must be solved for. Corresponding to these different pressures are known masses m _{1}, m _{2}, m _{3}, m _{4}, m _{5} and known mass densities ρ _{1}, ρ _{2}, ρ _{3}, ρ _{4}, ρ _{5} for each of the stacked weights. For simplicity assume that the atmospheric air density ρ _{ atm } and local gravitational acceleration are both constant, while A _{0} and λ are known equation parameters that were previously experimentally determined.
Rearrange the known nonlinear equation as a system of five equations(A.2)
In the above system of equations the input X will be a 12 × 1 vector(A.3)and the output Y will be a 5 × 1 vector(A.4)
Each of the five scalar components of the nonlinear equations h _{ j }( Y , X ) = 0 from the vector equation h ( Y , X ) = 0 would have to be simultaneously solved for specified values of X in order to determine Y .
In an actual laboratory the mass values m _{ i } would always be correlated to the mass density values ρ _{ i } since both measured masses and measured mass densities are ultimately traceable to a country's national kilogram standard, and from pressure balance theory the zeropressure area A _{0} and distortion coefficient λ are also always statistically correlated. This means in practical terms that the corresponding covariance matrix U _{ X } which has elements of covariances like Cov(m _{2}, ρ _{2}) and Cov(A _{0}, λ) will have various nonzero elements.
Since the measurement equation h ( Y , X ) = 0 may contain many inputs and outputs, it may be observed that the use of a symbolic algebra calculation of the various sensitivity coefficients produces considerable benefits for constructing the matrices C _{ X } and C _{ Y } in order to solve for the covariance matrix U _{ Y }.
The previous Python code for the input sensitivity matrix C _{ X } may simply be edited as follows:
Similarly the rows of the matrix C _{ X } may simply be calculated as
and then assembled as
Similarly the Python code for the output sensitivity matrix C _{ Y } may in turn simply be edited as follows:
and assembled as
The above symbolic matrices C _{ X } and C _{ Y } may then be numerically evaluated using the supplied information, and the final uncertainty matrix U _{ Y } calculated using the linear algebra algorithm from Figure 6. The final result in the illustrative example of a pressure balance will be the solutions P _{1} from h _{1}( Y , X ) through to P _{5} from h _{5}( Y , X ) along with a covariance matrix of pressures that indicates the pressure uncertainties u ^{2}(P _{1}) through to u ^{2}(P _{5}) along with the covariances Cov(P _{ i }, P _{ j }) of the form
In the particular case of a pressure calibration the pressure balance determined known values of pressures Y = [P _{1}, … , P _{5}]^{ T } would then be used to calibrate either a digital or analogue pressure transducer or other simlar pressure measurement UUT device. Depending on the linearity/nonlinearity of the pressure UUT being calibrated with UUT measured pressure values Z = [Z _{1}, … , Z _{5}]^{ T } the above covariance matrix may either be used to fit a straight line calibration curve e.g. P _{ UUT } = aP _{ true } + b or possibly a nonlinear calibration curve e.g. with a model Z = G(Y) and corresponding UUT covariance matrix U _{ Z }.
References
 BIPM, IEC, IFCC, ILAC, ISO, IUPAP, and OIML, Evaluation of measurement data − Guide to the expression of uncertainty in measurement, tech. rep., JCGM/WG1 GUM (2008). Revised 1st edition − https://www.bipm.org/en/publications/guides/ [Google Scholar]
 S. Kline, F. Mcclintock, Describing uncertainties in singlesample experiments, Mech. Eng. 75 , 3–8 (1953) [Google Scholar]
 BIPM, IEC, IFCC, ILAC, ISO, IUPAP, and OIML, Evaluation of measurement data − Supplement 1 to the “Guide to the expression of uncertainty in measurement” − Propogation of distributions using a Monte Carlo method, tech. rep., JCGM/WG1 GUM Supplement 1, 2008. 1st edition −https://www.bipm.org/en/publications/guides/ [Google Scholar]
 R. Willink, Representating Monte Carlo output distributions?for transfereability in uncertainty analysis: modelling with quantile functions, Metrologia 46 , 154–166 (2009) [CrossRef] [Google Scholar]
 V. Ramnath, Application of quantile functions for the analysis and comparison of gas pressure balance uncertainties, Int. J. Metrol. Qual. Eng. 8 , 4 (2017) [CrossRef] [EDP Sciences] [Google Scholar]
 BIPM, IEC, IFCC, ILAC, ISO, IUPAP, and OIML, Evaluation of measurement data − Supplement 2 to the “Guide to the expression of uncertainty in measurement” − Propogation of distributions using a Monte Carlo method, tech. rep., JCGM/WG1 GUM Supplement 2 (2011). 1st edition −https://www.bipm.org/en/publications/guides/ [Google Scholar]
 J.S. Hansen, GNU Octave Beginner’s Guide (PACKT Publishing, 2011) [Google Scholar]
 V. Ramnath, Comparison of straight line curve fit approaches for determining variances and covariances, Int. J. Metrol. Qual. Eng. 11 , 16 (2020) [CrossRef] [EDP Sciences] [Google Scholar]
 V. Ramnath, Analysis and comparison of hyperellipsoidal and smallest coverage regions for multivariate Monte Carlo measurement uncertainty analysis simulation datasets, MAPAN J. Metrol. Soc. India 2019 , 1–16 (2019) [Google Scholar]
 M.G. Cox, B.R.L. Siebert, The use of a Monte Carlo method for evaluating uncertainty and expanded uncertainty, Metrologia 43 , S178–S188 (2006) [Google Scholar]
 V. Ramnath, Numerical analysis of the accuracy of bivariate quantile distributions utilizing copulas compared to the GUM supplement 2 for oil pressure balance uncertainties, Int. J. Metrol. Qual. Eng. 8 , 29 (2017) [CrossRef] [EDP Sciences] [Google Scholar]
 P.M. Harris, M.G. Cox, On a Monte Carlo method for measurement uncertainty evaluation and its implementation, Metrologia 51 , S176–S182 (2014) [Google Scholar]
 Z.L. Warsza, J. Puchalski, Estimation of uncertainties in indirect multivariable measurements: Part 2. influence of the processing function accuracy, in Automation 2020 Towards Industry of the Future, edited by R. Szewczyk, C. Zielinski, M. Kaliczynska (Springer, Warsaw, Poland, 2020), pp. 326–344 [Google Scholar]
 A. Forbes, Approximate models of CMM behavior and point cloud uncertainties, Measurement: Sensors 18 , 100304 (2021) [CrossRef] [Google Scholar]
 M.G. Cox, A.M.H. van der Veen, Understanding and treating correlated quantities in measurement uncertainty evaluation, in Good practice in evaluating measurement uncertainty: Compendium of examples, edited A.M.H. van der Veen, M.G. Cox (EURAMET/NPL, 2021), pp. 29–44 [Google Scholar]
 N.V. Kornilov, V.G. Pronyaev, S.M. Grimes, Systematic distortion factor and unrecognized source of uncertanties in nuclear data measurements and evaluations, Metrology 2 , 1–18 (2021) [CrossRef] [Google Scholar]
 W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 2007), 3rd edn. [Google Scholar]
 A. Chhabra, J.R. Venepally, D. Kim, Measurement noise covarianceadapting Kalman filters for varying sensor noise situations, Sensors 21 , 8304 (2021) [CrossRef] [PubMed] [Google Scholar]
 O.M. Brastein, A. Ghaderi, C.F. Pfeiffer, N.O. Skeie, Analysing uncertainty in parameter estimation and prediction for greybox building thermal behavior models, Energy Build. 224 , 110236 (2020) [CrossRef] [Google Scholar]
 I. Smith, Y. Luo, D. Hutzschenreuter, The storage within digital calibration certificates of uncertainty information obtained using a Monte Carlo method, Metrology 2 , 33–45 (2022) [CrossRef] [Google Scholar]
 P. Saunders, Uncertainty propogation through integrated quantities for radiation thermometry, Metrologia 55 , 863–871 (2018) [CrossRef] [Google Scholar]
 P. Saunders, D.R. White, Physical basis of interpolation equations for radiation thermometry, Metrologia 40 , 195–203 (2003) [CrossRef] [Google Scholar]
 S. Briaudeau, M. Sadli, F. Bourson, B. Rougi, A. Rihan, J.J. Zondy, Primary radiometry for the miseenpratique: the laserbased radiance method applied to a pyrometer, Int. J. Thermophys. 32 , 2183–2196 (2011) [CrossRef] [Google Scholar]
 H.W. Yoon, V.B. Khromchenko, G.P. Eppeldauer, C.E. Gibson, J.T. Woodward, P.S. Shaw, K.R. Lykke, Towards highaccuracy primary spectral radiometry from 400K to 1300 K, Philos. Trans. Roy. Soc. A 374 , 1–12 (2016) [Google Scholar]
 A. Manoi, P. Wongnut, X. Lu, P. Bloembergen, P. Saunders, Calibration of standard radiation thermometers using two fixed points, Metrologia 57 , 014002 (2020) [CrossRef] [Google Scholar]
 P. Saunders, General interpolation equations for the calibration of radiation thermometers, Metrologia 34 , 201–210 (1997) [CrossRef] [Google Scholar]
 P. Saunders, Propagation of uncertainty for nonlinear calibration equations with an application in radiation thermometry, Metrologia 40 , 93–101 (2003) [CrossRef] [Google Scholar]
 A. Meurer, C.P. Smith, M. Paprocki, O. Čertík, S.B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J.K. Moore, S. Singh, T. Rathnayake, S. Vig, B.E. Granger, R.P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F. Pedregosa, M.J. Curry, A.R. Terrel, S. Rouficka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman, A. Scopatz, Sympy: symbolic computing in python, PeerJ Comput. Sci. 3 , e103 (2017) [CrossRef] [Google Scholar]
 H. PrestonThomas, The international temperature scale of 1990 (ITS90), Metrologia 27 , 3–10 (1990) [CrossRef] [Google Scholar]
 F. Sakuma, S. Hattori, A practical type fixed point blackbody furnance, in Temperature: Its Measurement and Control in Science and Industry , edited by J.F. Schooley (AIP, New York, 1982), vol. 5 , pp. 535–539 [Google Scholar]
 E.R. Wooliams, G. Machin, D.H. Lowe, R. Winkler, Metal (carbide)carbon eutectics for thermometry and radiometry: a review of the first seven years, Metrologia 43 , R11–R25 (2006) [CrossRef] [Google Scholar]
 A.D.W. Todd, K. Anhalt, P. Bloembergen, B.B. Khlevnoy, D.H. Lowe, G. Machin, N. Sasajima, P. Saunders, On the uncertainties in the realization of the kelvin based on thermodynamic temperatures of hightemperature fixedpoint cells, Metrologia 58 , 035007 (2021) [CrossRef] [Google Scholar]
Cite this article as: Vishal Ramnath, Determining the covariance matrix for a nonlinear implicit multivariate measurement equation uncertainty analysis, Int. J. Metrol. Qual. Eng. 13, 9 (2022)
All Figures
Fig. 1 Illustration of measurement traceability scheme showing how information from a covariance matrix of a national laboratory standard is transmitted to a client instrument measurement standard in industry to propagate uncertainty analysis results. 

In the text 
Fig. 2 Illustration of the full width half maximum of a curve [Graphic Source: Arne Nordmann (norro), CC BYSA 3.0 http://creativecommons.org/licenses/bysa/3.0/, via Wikimedia Commons]. 

In the text 
Fig. 3 Numerically stable linear algebra algorithm for determining covariance matrix. 

In the text 
Fig. 4 Representative pyrometer relative spectral responsitivity data reported by Saunders [21]. 

In the text 
Fig. 5 Python code for symbolically and numerically constructing the matrices C _{ x }, C _{ y } and U _{ x }. 

In the text 
Fig. 6 Python code for linear algebra implementation to postprocess C _{ x }, C _{ y } and U _{ x } to calculate U _{ y }. 

In the text 
Fig. 7 Parity line curve to validate the accuracy of the model fit for the SakumaHattori equation. 

In the text 
Fig. 8 Comparison of signaltemperature results from specified and predicted data for the SakumaHattori equation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.