© V. Ramnath, Published by EDP Sciences, 2022

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 ₁, … , X _N ), or equivalently Y = f ( X ) with X  = [X ₁, … , X _N ]^T ∈ ℝ ^N , where X ₁, … , 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 $x_i\sim N({\mu }_i,{\sigma }_i^2)$ where μ _i and ${\sigma }_i^2$ are corresponding equivalent expected and variance values for X _i for i = 1, … , N using standard statistical transformations to convert a non-Gaussian 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 Welch-Satterthwaite formula for estimating an effective degrees-of-freedom 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 Welch-Satterthwaite formula such that g _Y (y) is approximated as a Student's t-distribution 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 non-Gaussian 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 non-Gaussian 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 non-Gaussian 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 non-Gaussian PDFs for both inputs as well as outputs.

An explicit multivariate equation of the form$$Y=f(X)$$(1)has an known input X  = [X ₁, … , 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 ₁, … , 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$${\text{U}}_y={\text{C}}_x{\text{U}}_x{\text{C}}_x^{\text{T}}$$(2)

For the above equations the matrices are specified as$$U_x=\left[\begin{array}{lll}u(X_1,X_1)\hfill & \cdots \hfill & u(X_1,X_N)\hfill \\ ⋮\hfill & \ddots \hfill & ⋮\hfill \\ u(X_N,X_1)\hfill & \cdots \hfill & u(X_N,X_N)\hfill \end{array}\right]$$(3) $$C_x=\left[\begin{array}{lll}\frac{\partial f_1}{\partial X_1}\hfill & \cdots \hfill & \frac{\partial f_1}{\partial X_N}\hfill \\ \hfill ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \frac{\partial f_m}{\partial X_1}\hfill & \cdots \hfill & \frac{\partial f_m}{\partial X_N}\hfill \end{array}\right]$$(4) $${\text{U}}_y=\left[\begin{array}{lll}u(Y_1,Y_1)\hfill & \cdots \hfill & u(Y_1,Y_m)\hfill \\ ⋮\hfill & \ddots \hfill & ⋮\hfill \\ u(Y_m,Y_1)\hfill & \cdots \hfill & u(Y_m,Y_m)\hfill \end{array}\right]$$(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 $\frac{\partial f_i}{\partial x_j}$ 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$${\text{U}}_y={\text{C}}_x*{\text{U}}_x*{\text{C}}_x\prime$$

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 ${\text{C}}_Y=\left[\begin{array}{ll}{u}^2(a)\hfill & Cov(a,b)\hfill \\ Cov(a,b)\hfill & u^2(b)\hfill \end{array}\right]$ 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 data-points 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 $u^2(y)={\left(\frac{\partial y}{\partial x}\right)}^2{u}^2(x)+{\left(\frac{\partial y}{\partial a}\right)}^2{u}^2(a)+{\left(\frac{\partial y}{\partial b}\right)}^2{u}^2(b)+2\left(\frac{\partial y}{\partial a}\right)\left(\frac{\partial y}{\partial b}\right)Cov(a,b)$.

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$$h(Y,X)=0$$(6)the uncertainty u( Y ) may now be obtained from the covariance matrix U _y that is specified from the matrix equation$${\text{C}}_y{\text{U}}_y{\text{C}}_{\text{y}}^{\text{T}}={\text{C}}_x{\text{U}}_x{C}_x^{\text{T}}$$(7)

In the above equation the input sensitivity matrix C _x is of dimension m × N with matrix elements C _ij defined as $C_{ij}=\frac{\partial h_i}{\partial X_j}$ 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 $C_{\ell k}=\frac{\partial h_{\ell }}{\partial Y_k}$ where l = 1, … , m and k = 1, … , m, such that$${\text{C}}_x\left[\begin{array}{lll}\frac{\partial h_1}{\partial X_1}\hfill & \cdots \hfill & \frac{\partial h_1}{\partial X_N}\hfill \\ ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \frac{\partial h_m}{\partial X_1}\hfill & \cdots \hfill & \frac{\partial h_m}{\partial X_N}\hfill \end{array}\right]$$(8) $${\text{C}}_y\left[\begin{array}{lll}\frac{\partial h_1}{\partial Y_1}\hfill & \cdots \hfill & \frac{\partial h_1}{\partial Y_m}\hfill \\ ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \frac{\partial h_m}{\partial Y_1}\hfill & \cdots \hfill & \frac{\partial h_m}{\partial Y_N}\hfill \end{array}\right]$$(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 post-process 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 ${\text{C}}_y{\text{U}}_y{\text{C}}_y^{\text{T}}={\text{C}}_x{\text{U}}_x{\text{C}}_x^{\text{T}}$. 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$$\begin{array}{l}{g}_Y(\mathit{\eta })={\int }_{-\infty }^{\infty }{g}_X(\mathbf{\xi })\hfill \\ \times \delta (\mathit{\eta }-f(\mathbf{\xi }))d{\xi }_1\cdots d{\xi }_N\hfill \end{array}$$(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 post-processing 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 three-dimensional vector measurement equation model based on a Sakuma-Hattori pyrometry equation where the temperature-signal has the form$$S(T)=\frac{C}{\exp \left(\frac{c_2}{AT+B}\right)-1}$$(11)

The objective of this paper is to determine the covariance matrix for the parameters A, B and C in the equation$$\begin{array}{lll}{u}^2(A)\hfill & Cov(A,B)\hfill & Cov(A,C)\hfill \\ Cov(A,B)\hfill & u^2(B)\hfill & Cov(B,C)\hfill \\ Cov(A,C)\hfill & Cov(B,C)\hfill & u^2(C)\hfill \end{array}$$(12)

When the above covariance matrix is determined it may then conveniently be used in the corresponding inverse signal equation for the temperature$$T=\frac{c_2}{A\ln \left(\frac{C}{S}+1\right)}-\frac{B}{A}$$(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 Sakuma-Hattori 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 ₁, X ₂, … , X _N ]^T and Y = [Y ₁, Y ₂, … , Y _m ]^T. An example of a multivariate measurement equation is the case of a multi-parameter digital multimeter which is a single instrument that can measure voltage V _meas  = Y ₁, current I _meas  = Y ₂ and resistance R _meas  = Y ₃ based on input standards of resistance R _source  = X ₁ and voltage V _source  = X ₂. The inputs X ₁ and X ₂ by necessity include both nominal expected values of these quantities as well as any correlation Cov(X ₁, X ₂) 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 ₁, Y ₂ and Y ₃ as well as their corresponding uncertainties. Covariances conveniently model the output uncertainties and are defined by the terms Cov(Y ₁, Y ₂), Cov(Y ₁, Y ₃) and Cov(Y ₂, Y ₃) on the off-diagonal elements of U _Y , whilst the variances Var(Y ₁, Y ₁), Var(Y ₂, Y ₂) and Var(Y ₃, Y ₃) are located on the diagonal elements of the U _Y matrix. In this work the multivariate uncertainty of the output was approximated as $\Delta \text{Y}\approx \frac{\partial \text{F}(\mathbf{\text{X}})}{\partial \mathbf{\text{X}}}\equiv \text{S}(\mathbf{\text{X}})\cdot \Delta \text{X}(\mathbf{\text{X}})$ where the elements of the sensitivity matrix $S_{ji}=\frac{\partial y_i}{\partial x_i}$ of dimension m × N propagate the output uncertainty as U _Y  = S ⋅ U _X  ⋅ S ^T and the input and output matrices are$${\text{U}}_{\text{X}}=\left[\begin{array}{lll}{u}^2(x_1)\hfill & \dots \hfill & \rho (x_1,x_N)u(x_1)u(x_N)\hfill \\ ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \rho (x_1,x_N)u(x_1)u(x_N)\hfill & \dots \hfill & u^2(x_N)\hfill \end{array}\right]$$ $${\text{U}}_{\text{Y}}=\left[\begin{array}{lll}\hfill u^2(y_1)\hfill & \dots \hfill & \rho (y_1,y_m)u(y_1)u(y_m)\hfill \\ ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \rho (y_1,y_m)u(y_1)u(y_m)\hfill & \dots \hfill & u^2(y_m)\hfill \end{array}\right]$$

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 ${\mathbf{\text{U}}}_{\mathbf{\text{Z}}}={\mathbf{S}}_{\mathbf{\text{G}}}\cdot {\mathbf{\text{U}}}_{\mathbf{\text{Y}}}\cdot {\mathbf{\text{U}}}_{\mathbf{\text{G}}}^{\mathbf{\text{T}}}$. 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 two-stage 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 Sakuma-Hattori 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 χ ² 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 post-processing 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$$S=K{\int }_0^{\infty }S(\lambda )L_b(\lambda ,T)d\lambda$$(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$$L_b(\lambda ,T)=\frac{c_1}{n^2{\lambda }^5\left[\exp \left(\frac{c_2}{n\lambda T}\right)-1\right]}$$(15)

For the above equation c ₁ and c ₂ 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 Sakuma-hattori 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 Sakuma-Hattori equation. If the covariance matrix is known then the uncertainty in temperature from the pyrometer signal may then be calculated as$$\begin{array}{l}{u}^2(T)={\left(\frac{\partial T}{\partial S}\right)}^2{u}^2(S)+{\left(\frac{\partial T}{\partial A}\right)}^2{u}^2(A)\hfill \\ +{\left(\frac{\partial T}{\partial B}\right)}^2{u}^2(B)+{\left(\frac{\partial T}{\partial C}\right)}^2{u}^2(C)\hfill \\ +2\frac{\partial T}{\partial A}\frac{\partial T}{\partial B}Cov(A,B)+2\frac{\partial T}{\partial A}\frac{\partial T}{\partial C}Cov(A,C)\hfill \\ +2\frac{\partial T}{\partial B}\frac{\partial T}{\partial C}Cov(B,C)\hfill \end{array}$$(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 non-linear regression analysis and a trial-and-error approach must be adopted to obtain suitable starting values for a subsequent Levenberg-Marquardt 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 $r=\frac{\sigma }{{\lambda }_0}<<1$ is small as discussed by Saunders & White [22] the parameters A, B, C may then be approximated as$$A=n{\lambda }_0\left[1-6{\left(\frac{\sigma }{{\lambda }_0}\right)}^2\right]$$(17) $$B=\frac{c_2}{2}{\left(\frac{\sigma }{{\lambda }_0}\right)}^2$$(18) $$C=K\frac{c_1}{n^2{\lambda }_0^5}\left[1+15{\left(\frac{\sigma }{{\lambda }_0}\right)}^2\right]{\int }_0^{\infty }s(\lambda )d\lambda$$(19)

In the above equation λ ₀ 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 quasi-Lambertian 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⁻¹ ⋅ m⁻³ into a measured photo-current signal in units of A. For this measurement system the pyrometer photo-current I _pyro(λ, r = 0)/[A] is specified as$$\begin{array}{l}{I}_{\text{pyro}}(T)=\left(1+C_{\text{pyro}}^{\text{NL}}\right)(1+SSE)ϵ\\ \times {\int }_0^{\infty }{L}_{bb}(\lambda ,T)S_{\text{pyro}}(\lambda )d\lambda \end{array}$$(20)where the absolute spectral responsivity is $S_{\text{pyro}}(\lambda )/[\text{A}\cdot {\text{W}}^{-1}\cdot \text{s}{\text{r}}^{-1}\cdot {\text{m}}^3,C_{\text{pyro}}^{\text{NL}}$ is a non-linearity correction factor for pyrometer responsivity, SSE is the size-of-source effect for the pyrometer, and ε is constant blackbody emissivity.

In the special case the above calibration terms $C_{\text{pyro}}^{\text{NL}}$, 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 Sakuma-Hattori 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 current-to-voltage 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 temperature-signal 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 λ ₀ 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 $\frac{\sigma }{{\lambda }_0}<0.01$. 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$$S=c^{\prime }{\left[\exp \left(\frac{c_2}{{\lambda }_0\left(1-6\frac{{\sigma }^2}{{\lambda }_0^2}\right)T+\frac{c_2{\sigma }^2}{2{\lambda }_0^2}}\right)-1\right]}^{-1}$$(21)

When these special conditions apply the constants in a Sakuma-Hattori equation may be obtained using c ₂ = 0.014388 m ⋅ K as C = c ^′, $A={\lambda }_0\left(1-6\frac{{\sigma }^2}{{\lambda }_0^2}\right)$, and $B=\frac{c_2}{2}\frac{{\sigma }^2}{{\lambda }_0^2}$ 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 λ ₀ and c ^′ and these can then be simultaneously determined from the pyrometer signals at the two known temperature fixed points using the (T ₁, S ₁) and (T ₂, S ₂) pairs of experimentally measured signal points where T ₁ and T ₂ are known fixed points by definition.

An earlier similar concept by Saunders [26] using two temperatures T ₀ and T ₁ and the temperature-signal relationship $S(T)={\int }_0^{\infty }R(\lambda )L_b(\lambda ,T)d\lambda$, 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 temperature-dependent mean effective wavelength λ _m so that the ratio of signals may be expressed as$$Q=\frac{S(T_1)}{S(T_0)}=\frac{L_b({\lambda }_m,T_1)}{L_b({\lambda }_m,T_0)}$$(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$$\frac{1}{{\lambda }_T}=\frac{{\int }_0^{\infty }\frac{1}{\lambda }R(\lambda )L_b(\lambda ,T)d\lambda }{{\int }_0^{\infty }R(\lambda )L_b(\lambda ,T)d\lambda }$$(23)from which λ _m is finally obtained as$$\frac{1}{{\lambda }_m}=\frac{1}{\frac{1}{T_1}\frac{1}{T_0}}{\int }_{\frac{1}{T_0}}^{\frac{1}{T_1}}\frac{1}{{\lambda }_T}d\left(\frac{1}{T}\right)$$(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 time-consuming calibration that must be performed in a spectral radiometry laboratory with expensive specialist equipment.

Saunders [27] considered a univariate non-linear regression for pairs of data (x ₁, y ₁), … , (X _N , y _N ) where the model took the form$$y=y\left(x;x_1,\dots ,x_N,y_1,\dots ,y_N\right)$$(25)and obtained expressions with M data points in the absence of correlations to fit a a Sakuma-Hattori equation where$$u_y^2=\sum _{i=1}^{\mathrm{M}}{\left(\frac{\partial y}{\partial x_i}\right)}^2{u}_{x_i}^2+\sum _{i=1}^{\mathrm{M}}{\left(\frac{\partial y}{\partial y_i}\right)}^2{u}_{y_i}^2+{\left(\frac{\partial y}{\partial x}\right)}^2{u}_x^2$$(26)where the sensitivity terms $\frac{\partial y}{\partial x_i}$ and $\frac{\partial y}{\partial y_i}$ 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 Sakuma-Hattori 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 BY-SA 3.0 http://creativecommons.org/licenses/by-sa/3.0/, via Wikimedia Commons].

3 Mathematical modelling

In order to solve the matrix equation ${\mathbf{\text{C}}}_{\mathbf{y}}{\mathbf{\text{U}}}_{\mathbf{y}}{\mathbf{\text{C}}}_{\mathbf{y}}^{\mathbf{\text{T}}}={\mathbf{\text{C}}}_{\mathbf{x}}{\mathbf{\text{U}}}_{\mathbf{x}}{\mathbf{\text{C}}}_{\mathbf{y}}^{\mathbf{\text{T}}}$ 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$${\mathbf{\text{U}}}_{\mathbf{y}}=\mathbf{C}{\mathbf{U}}_{\mathbf{x}}{\mathbf{C}}^{\mathbf{\text{T}}}$$(27) $$\mathbf{\text{C}}={\mathbf{\text{C}}}_y^{-1}{\mathbf{\text{C}}}_{\mathbf{\text{x}}}$$(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 non-specialist routines to calculate U _y .

Considering the general form of the Sakuma-Hattori equation and three pairs of experimental data points (T ₁, S ₁), (T ₂, S ₂) and (T ₃, S ₃) the equation to determine the parameters may be constructed as$${\chi }^2={\displaystyle \sum _{i=1}^3}{[S(T_i)-S_i]}^2$$(29) $$X={[T_1,T_2,T_3,S_1,S_2,S_3]}^{\text{T}}$$(30) $$Y={[A,B,C]}^{\text{T}}$$(31) $$\mathbf{h}(\mathbf{Y},X)=\left[\begin{array}{c}\hfill \frac{\partial {\chi }^2}{\partial A}\hfill \\ \hfill \frac{\partial {\chi }^2}{\partial B}\hfill \\ \hfill \frac{\partial {\chi }^2}{\partial C}\hfill \end{array}\right]$$(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 merit-function expression results:$$\begin{array}{l}{\chi }^2={\left(-X_4+\frac{Y_3}{\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)-1}\right)}^2\hfill \\ +{\left(-X_5+\frac{Y_3}{\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)-1}\right)}^2\hfill \\ +{\left(-X_6+\frac{Y_3}{\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)-1}\right)}^2\hfill \end{array}$$(33)

The system of three non-linear equations is then obtained by partial differentiation as:$$\begin{array}{l}{h}_1(Y_1,Y_2,Y_3,X_1,X_2,X_3,X_4,X_5,X_6)\hfill \\ =\frac{2X_1{Y}_3{c}_2\left(-X_4+\frac{Y_3}{\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)-1}\right)\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)}{{(X_1{Y}_1+Y_2)}^2{\left(\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)-1\right)}^2}\hfill \\ +\frac{2X_2{Y}_3{c}_2\left(-X_5+\frac{Y_3}{\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)-1}\right)\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)}{{(X_2{Y}_1+Y_2)}^2{\left(\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)-1\right)}^2}\hfill \\ +\frac{2X_3{Y}_3{c}_2\left(-X_6+\frac{Y_3}{\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)-1}\right)\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)}{{(X_3{Y}_1+Y_2)}^2{\left(\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)-1\right)}^2}\hfill \end{array}$$(34) $$\begin{array}{l}{h}_2(Y_1,Y_2,Y_3,X_1,X_2,X_3,X_4,X_5,X_6)\hfill \\ =\frac{2Y_3{c}_2\left(-X_4+\frac{Y_3}{\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)-1}\right)\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)}{{(X_1{Y}_1+Y_2)}^2{\left(\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)-1\right)}^2}\hfill \\ +\frac{2Y_3{c}_2\left(-X_5+\frac{Y_3}{\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)-1}\right)\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)}{{(X_2{Y}_1+Y_2)}^2{\left(\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)-1\right)}^2}\hfill \\ +\frac{2Y_3{c}_2\left(-X_6+\frac{Y_3}{\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)-1}\right)\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)}{{(X_3{Y}_1+Y_2)}^2{\left(\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)-1\right)}^2}\hfill \end{array}$$(35) $$\begin{array}{l}{h}_3(Y_1,Y_2,Y_3,X_1,X_2,X_3,X_4,X_5,X_6)\hfill \\ =\frac{2\left(-X_4\right)+\frac{Y_3}{\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)-1}}{\exp \left(\frac{c_2}{X_1{Y}_1+Y_2}\right)-1}\hfill \\ +\frac{2\left(-X_5\right)+\frac{Y_3}{\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)-1}}{\exp \left(\frac{c_2}{X_2{Y}_1+Y_2}\right)-1}\hfill \\ +\frac{2\left(-X_6\right)+\frac{Y_3}{\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)-1}}{\exp \left(\frac{c_2}{X_3{Y}_1+Y_2}\right)-1}\hfill \end{array}$$(36)

The input sensitivity matrix C _X with elements specified by the partial derivatives $\frac{\partial h_i}{\partial X_j}$ 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 $\frac{\partial h_i}{\partial Y_j}$ 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 ₁ = T ₁,  X ₂ = T ₂,  X ₃ = T ₃,  X ₄ = S ₁,  X ₅ = S ₂,  X ₆ = S ₃ and Y ₁ = A,  Y ₂ = B,  Y ₃ = C where the numerical values of Y ₁, Y ₂ and Y ₃ are obtained from the numerical solution of h ₁( Y , X ) = 0,  h ₂( Y , X ) = 0,  and h ₁( 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 (ITS-90) temperature data as reported by Preston-Thomas [29] let the silver, gold and copper fixed point freezing temperatures be T ₁ = T _Ag  = 1234.93 K,  T ₂ = T _Au  = 1337.33 K and T ₃ = T _Cu  = 1357.77 K and assume medium accuracy standard uncertainties of u(T ₁) = ± 0.3 K,  u(T ₂) = ± 0.3 K and u(T ₃) = ± 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 ₁ = 3.741771852 × 10⁻¹⁶ W ⋅ m² and c ₂ = 1.438776877 × 10⁻² m ⋅ K.

From this data the pyrometer signal values may be calculated using the formula $S(T)={\int }_0^{\infty }R(\lambda )L_b(\lambda ,T)d\lambda$ to yield the following underlying temperature-signal data to fit the parameters in a Sakuma-Hattori equation such that$$\left[\begin{array}{c}\hfill T_1\hfill \\ \hfill T_2\hfill \\ \hfill T_3\hfill \end{array}\right]=\left[\begin{array}{c}\hfill \mathrm{1234.93\; 0000000000}\hfill \\ \hfill \mathrm{1337.33\; 0000000000}\hfill \\ \hfill \mathrm{1357.77\; 0000000000}\hfill \end{array}\right]$$(37)and$$\left[\begin{array}{c}\hfill S_1\hfill \\ \hfill S_2\hfill \\ \hfill S_3\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 5.955760804302088e-01\hfill \\ \hfill 2.346969884606723e+00\hfill \\ \hfill 3.010492321743183e+00\hfill \end{array}\right]$$(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.191152160860418×10^-2  u(S ₂) = ± 4.693939769213445×10^-2 and u(S ₃) = ± 6.020984643486366×10^-2.

Setting the input as X  = [T ₁, T ₂, T ₃, S ₁, S ₂, S ₃]^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 ${\mathbf{\text{C}}}_{\mathbf{y}}{\mathbf{\text{U}}}_{\mathbf{y}}{\mathbf{\text{C}}}_{\mathbf{y}}^{\mathbf{\text{T}}}={\mathbf{\text{C}}}_{\mathbf{x}}{\mathbf{\text{U}}}_{\mathbf{x}}{\mathbf{\text{C}}}_{\mathbf{x}}^{\mathbf{\text{T}}}$ 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 signal-temperature values (S ₁, T ₁),  (S ₂, T ₂),  (S ₃, T ₃) are compared to the corresponding predictions for the temperatures as obtained from the Sakuma-Hattori 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 Dirac-delta point function. In general pyrometry practice there is usually a balance or trade-off 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 signal-temperature 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. 4 Representative pyrometer relative spectral responsitivity data reported by Saunders [21].

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 post-process 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 Sakuma-Hattori equation.

Fig. 8 Comparison of signal-temperature results from specified and predicted data for the Sakuma-Hattori 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 dead-weight 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 counter-balanced by a hydrostatic pressure P _i that pushes the piston in an upwards direction. By a simple free-body diagram and equilibrium of forces it follows that the governing equation is$$P_i=\frac{m_i\left(1-\frac{{\rho }_{atm}}{{\rho }_i}\right)g}{A_0(1+\lambda P_i)}$$(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 ₁, P ₂, P ₃, P ₄, P ₅ which must be solved for. Corresponding to these different pressures are known masses m ₁, m ₂, m ₃, m ₄, m ₅ and known mass densities ρ ₁, ρ ₂, ρ ₃, ρ ₄, ρ ₅ for each of the stacked weights. For simplicity assume that the atmospheric air density ρ _atm and local gravitational acceleration are both constant, while A ₀ and λ are known equation parameters that were previously experimentally determined.

Rearrange the known non-linear equation as a system of five equations$$\begin{array}{ccc}\hfill h_1\hfill & \hfill =\hfill & \hfill P_1-\frac{m_i\left(1-\frac{{\rho }_{atm}}{{\rho }_1}\right)g}{A_0(1+\lambda P_1)}\hfill \\ \hfill \hfill & \hfill ⋮\hfill & \hfill \hfill \\ \hfill h_5\hfill & \hfill =\hfill & \hfill P_5-\frac{m_i\left(1-\frac{{\rho }_{atm}}{{\rho }_5}\right)g}{A_0(1+\lambda P_5)}\hfill \end{array}$$(A.2)

In the above system of equations the input X will be a 12 × 1 vector$$X=\left[\begin{array}{l}{m}_1\hfill \\ ⋮\hfill \\ m_5\hfill \\ {\rho }_1\hfill \\ ⋮\hfill \\ {\rho }_5\hfill \\ A_0\hfill \\ \lambda \hfill \end{array}\right]=\left[\begin{array}{l}{X}_1\hfill \\ ⋮\hfill \\ X_5\hfill \\ X_6\hfill \\ ⋮\hfill \\ X_{10}\hfill \\ X_{11}\hfill \\ X_{12}\hfill \end{array}\right]$$(A.3)and the output Y will be a 5 × 1 vector$$Y=\left[\begin{array}{c}\hfill P_1\hfill \\ \hfill ⋮\hfill \\ \hfill P_5\hfill \end{array}\right]=\left[\begin{array}{c}\hfill Y_1\hfill \\ \hfill ⋮\hfill \\ \hfill Y_5\hfill \end{array}\right]$$(A.4)

Each of the five scalar components of the non-linear 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 zero-pressure area A ₀ 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 ₂, ρ ₂) and Cov(A ₀, λ) will have various non-zero 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 ₁ from h ₁( Y , X ) through to P ₅ from h ₅( Y , X ) along with a covariance matrix of pressures that indicates the pressure uncertainties u ²(P ₁) through to u ²(P ₅) along with the covariances Cov(P _i , P _j ) of the form

$${\text{U}}_Y=\left[\begin{array}{cccc}\hfill \mathrm{Cov}(P_{1}P2)\hfill & \hfill Cov(P_1,P_2)\hfill & \hfill \dots \hfill & \hfill Cov(P_1,P_5)\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill Cov(P_1,P_5)\hfill & \hfill Cov(P_2,P_5)\hfill & \hfill \dots \hfill & \hfill Cov(P_5,P_5)\hfill \end{array}\right]$$(A.5)

In the particular case of a pressure calibration the pressure balance determined known values of pressures Y  = [P ₁, … , P ₅] ^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/non-linearity of the pressure UUT being calibrated with UUT measured pressure values Z  = [Z ₁, … , Z ₅] ^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 non-linear calibration curve e.g. $P_{UUT}=\alpha P_{true}^2+\beta P_{true}+\gamma$ with a model Z = G(Y) and corresponding UUT covariance matrix U _Z .

References

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 BY-SA 3.0 http://creativecommons.org/licenses/by-sa/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 post-process 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 Sakuma-Hattori equation.
In the text

	Fig. 8 Comparison of signal-temperature results from specified and predicted data for the Sakuma-Hattori equation.
In the text