© V. Ramnath, published by EDP Sciences, 2020

1 Introduction

1.1 Research motivation

The use of least squares based statistical regression analysis is a widely used tool to analyse experimental data when estimating and fitting model parameter values and uncertainties. In the particular case of straight line equations of the form y = ax + b known input data x_i and output data y_i for data-points i = 1, … , n inclusive of possible available uncertainty information are used in the regression analysis where standard uncertainty analysis theory from the GUM [1] specifies the model uncertainty as(1) (2)

In the above formula since x is a statistically independent random variable it follows that cov(x, a)= 0, cov(x, b) = 0, and in general cov(a, b) ≠ 0. As a result, the probability density function (PDF) of the output g_y(η_y) may be determined where η_y is a univariate random variable of the model output y, g_x(ξ_x) is a univariate random variable of the model input x, and g_a,b(ξ_a, ξ_b) is a bivariate random variable of the joint PDF of the model parameters a and b that incorporates the coupling effect between random variables ξ_a and ξ_b for the model parameters a and b that are modelled in terms of cov(a, b).

The research problem that occurs when OLS, WLS and GLS regressions are used to estimate the parameters a and b to varying degrees of accuracies, is that additional Monte Carlo simulations must then be performed in order to determine the variances u²(a) and u²(b) along with the covariance cov(a, b), and that this information must then be further processed in order to infer u²(A₀), u²(λ) and cov(A₀, λ). Consequently, in the absence of covariance information the model uncertainty will be systematically under-estimated if the covariance is unknown or unavailable and only the variances of the model parameters are utilized in uncertainty calculations.

1.2 Research focus

This paper focuses on demonstrating the conceptual ease and functionality of the WTLSC approach and how it is able to incorporate measurement variance information u²(x_i) and u²(y_i) along with covariance information cov(x_i, y_i) of the data when available, in order to determine both the variances u²(a) and u²(b) along with the covariance information cov(a, b) of the model parameters for statistical completeness. A new simple analytical calculation approach to directly infer both the variances u²(A₀) and u²(λ) as well as covariance term cov(A₀, λ) from knowledge of u²(a), u²(b) and cov(a, b) as obtained from the WTLSC algorithm implemented with a Gnu Octave/Matlab program is presented that does not rely on additional Monte Carlo simulations or statistical post-processing.

2 Literature review

2.1 Physical theory of pressure balances

Pressure balance theory from Dadson [2] for a cross-floated pressure balance in equilibrium may be used to define the applied pressure P in terms of a generalized force F and an effective area S as where F is a generalized force term to account for buoyancy effects, and S is a cross-sectional effective area that accounts for pressure induced deformation and temperature expansion/contraction effects such that where S = A₀(1 + λP)φ, A₀ is a effective area when extrapolated to a zero pressure, λ is a distortion coefficient based on linear elasticity theory, and φ is an appropriate temperature compensation function that adjusts the area to an appropriate reference temperature.

In general since there is a pressure dependence within the effective area term S the equation is usually rearranged as a homogeneous equation and then solved for specified inputs to obtain a value of P such that h = 0. In the event that the pressure P is known and the area S must be determined the equation is rearranged as and again rearranged as another homogeneous equation in order to determine the value of S that solves the homogeneous equation for specified inputs.

In the above equation m is the mass on top of the piston, ρ_a is the ambient fluid density in contact with the piston to account for buoyancy forces acting on the piston, ρ_m is the mass density, g is the gravitational acceleration constant, A₀ is the zero-pressure area, λ is the distortion coefficient, P is the applied pressure, and φ = 1 + α(t − t₀) is a temperature compensation function to adjust the cross-sectional area to a reference temperature t₀ where α = α_p + α_c is the linear thermal expansion coefficient of the combination of the piston and cylinder where α_p and α_c are the corresponding linear thermal expansion coefficients of the piston and cylinder respectively.

When the above system of equations is implemented the equation for the cross-floated area S may then be simply modelled as(3)

In the above equation the inputs are the expected values and standard uncertainties of the applied pressures (μ(P_i) ± u(P_i)), the cross-floated areas (μ(S_i) ± u(S_i)), and estimates of the covariances cov(P_i, S_i) for a sequence of measurements i = 1, … , n where n is the number of measurements, whilst the unknowns are considered to be the expected values of the zero-pressure area A₀, the distortion coefficient λ, and the covariance cov(A₀, λ). This modelling approach is equivalent to fitting a straight line of the form y = ax + b from experimental data, and utilizing a least squares regression to quantify the values and uncertainties of the parameters a and b.

2.2 Historical and mathematical context

The fitting of straight lines is a common numerical task to fields such as science, engineering, and economics that require experimental or computational data to be modelled and fitted. Usually a least-squares algorithm is implemented to minimize a merit function as discussed by Saunders [3] where in general a non-linear calibration equation y = y(x ; a₁, … , a_N) is desired to be fitted from measurement pairs of controlled values x_i and observed values y_i. In the case where the number of measurements is M and the number of parameters to determine is N then a least squares solution is mathematically essential in the general case where M > N. For this situation a merit function is constructed as(4)and optimized in terms of the parameters a₁, … , a_N by simultaneously solving ∂χ²/∂ a_i = 0 for i = 1, … , N where w_i is a suitable and appropriate weighting function to account for the underlying measurement accuracy and available information.

The simplest case involves an optimization where all of the data points have an assumed constant uncertainty/accuracy such that u(x_i) = u(y_i) = σ for i = 1, … , M measurements, which corresponds to the situation where no detailed or specific measurement uncertainty information is available, and as a result the optimization of the merit function does not require any specific knowledge of σ as this can be conveniently factored out when optimizing χ². This particular problem thus corresponds to the well known historical total least squares problem as investigated by Van Huffel and Vandervalle [4].

At a next conceptual level the problem of fitting a straight line with errors in both coordinates was originally formulated by York [5] and subsequently involved either attempts to directly solve the original problem posed by York or by various approximation methods in order to infer the variances of the fitting parameters by researchers working in the area of regression analysis of experimental measurements. For this problem the approach of Press et al. [6] has achieved the most notable recognition amongst many non-mathematicians for the fitting of a straight line where both x and y coordinates possess some underlying known statistical uncertainty.

In this particular case the model y(x) = a + bx involves a χ² merit function optimization of the form(5)where u(x_i) and u(y_i) are the standard deviations of the x_i and y_i data points that are allowed to vary. Although conceptually straightforward the main challenge with directly solving this problem is that the presence of the parameter b in the denominator renders the partial derivative equation ∂χ²/∂ b = 0 a non-linear equation, and as a result this problem is considerably more difficult to solve than the earlier total least squares problem which is amenable to an analytical solution.

Consequently, the optimization problem for straight lines with errors in both coordinates is typically solved either with classical optimization techniques such as the Levenberg-Marquardt method or with the Monte Carlo bootstrap technique as outlined by Press et al. [6], and as a result this solution approach is not generally directly accessible for many metrologists working in industrial laboratories due to the mathematical complexity of the method.

It should be noted that a conventional numerical solution of this non-linear optimization problem using either classical techniques such as Newton's, steepest descent and homotopy methods as discussed by Burden and Faires [7] or newer methods such as the particle swarm optimization (PSO) approach would still suffer from two limitations, namely that it would firstly not be able to compute the parameter variances, and secondly that it would not be able to provide any insights into the parameter variances and covariances.

The statistical effects from the choice of OLS, WLS and GLS regressions was subsequently investigated by Duer et al. [8] who determined that whilst the specific choice of OLS, WLS or GLS regressions did not generally produce inconsistent estimates for the intercept b and slope a parameter values, that the different regression schemes nevertheless did in fact produce considerable variations in the corresponding estimates of the uncertainties u(a) and u(b) of the parameters a and b. As a result, Duer et al. concluded that a simplified OLS regression could only be safely utilized in the very specific and special case of where the uncertainties u(x_k) in the x_k data points was negligible.

In the particular case of scientific metrology work an inconsistent or incorrect estimate of a parameter value and variance/covariance could then have a catastrophic effect on high level laboratory inter-comparison work or in safety critical experimental measurements as parameter uncertainties of equipment and instruments are typically utilized in many quality engineering systems to specify tolerances and safety limits for equipment and personnel.

Due to these specific limitations particularly for scientific measurement work at national laboratories the total least squares problem was subsequently revisited by Krystek and Anton [9] who simplified the original two-dimensional problem into an equivalent one-dimensional optimization which is amenable with a set of simplified mathematical tools. This initial approach was then later generalized to solve the full problem of errors in both coordinates with correlation effects in a later paper by Krystek and Anton [10]. This approach is now termed as a weighted total least squares with correlation (WTLSC) analysis approach and it makes an assumption that the covariance matrix U_k for individual experimental data points for k = 1, … , n are available prior to optimizing the merit function which is explicitly specified as (6)

Part of the reasoning as to why the WTLSC approach does not appear to have been more widely adopted is that in the specific area of scientific metrology that the calculation of the covariance matrix is extremely mathematically complicated when using standard calculus based approaches, and is therefore only tractable from a full Monte Carlo pure numerical simulation approach as discussed by Harris [11]. As a result, if the covariance matrix is unavailable then the WTLSC approach becomes problematic in which case the covariance must either be calculated or estimated, or alternately modified by setting the correlations to zero such that the algorithm implementation then simplifies to the earlier weighted total least squares WLS algorithm.

Whilst the absence of a covariance matrix may be viewed as technical limitation when implementing a GLS regression where the standard statistical formulation comprises of determining the value of the regression parameter β which optimizes the residuals ε in a linear regression model y = X β + ε, it has subsequently been determined by researchers that this is not an insurmountable challenge as mechanisms exist to infer a likely covariance matrix. For this problem the covariance matrix Cov[ε |X] = Ω, is known to take the general form(7)

In order to account for the unavailability of the actual covariance matrix Ω which may be unknown, or impractical to estimate, one may nevertheless get a consistent estimate of Ω usually written as that is termed the feasible generalized least squares (FGLS) estimator using either iteration based commercial routines with software such as Matlab [12], or with newer support vector machines (SVMs) approaches that are currently under development such as recently reported results from Miller and Startz [13].

2.3 Recent statistical regression developments

Recently work reported in the literature by Wuethrich and Souiyam [14] building on earlier work by Otal and Yardin [15] has investigated the use of Monte Carlo based simulations for Ordinary Least Squares (OLS), Weighted Least Squares (WLS) and Generalized Least Squares (GLS) regression schemes in order to determine the parameter values of a pressure balance model of the form S = bP + a, where P is an input of applied pressure, S is an output of cross-floated area, and a and b are model parameters whose values and uncertainties must be determined. Conventionally the pressure balance model has the form S = A₀(1 + λP) where A₀ is a zero-pressure area and λ is a distortion coefficient and the regression analysis fits the values for a = A₀ and the product b = A₀ λ for i = 1, … , n data-points such that(8) (9)

For the above regression scheme the distortion coefficient is recovered as λ = b/a after the residuals ε_i are minimized in the matrix equation X C = Y + ε. This minimization is implemented through the use of a coefficient matrix V that is used to incorporate the known information such that(10) (11)

In an OLS scheme no weighting is applied and a coefficient matrix V consists of just a diagonal matrix with unity on each of the diagonal elements such that V_i,i = 1, i = 1,…,n in a WLS scheme the coefficient matrix consists of a diagonal matrix where the elements on each of the diagonal elements are proportional to the uncertainties of the areas i.e. V_i,i ∝ u²(S_i), and finally in a GLS scheme the diagonal elements of V equals the variances of the areas such that V_i,i = u²(S_i) whilst the off-diagonal elements equal the covariances such that V_i,j = cov(S_i, S_j),  i, j = 1, … , n,  i ≠ j.

The general approach of Wuethrich and Souiyam [14] is to solve equation (9) with nominal values of the input quantities at specified pressures P_i for the equivalent homogeneous equation in order to determine corresponding values of S_i for a sequence of pressures P₁, … , P_n with i = 1, … , n, and to then use this information to determine the regression parameter C_j = [(A₀)_j, b_j]^T from the matrix equation with an appropriate choice of weighting matrix V. This process is repeated M = 10000 times with slightly perturbed values of the inputs consistent with the underlying uncertainties to sequentially generate samples C_j for j = 1, … , M that then results in a sequence C₁, … , C_M of regression parameters, which are then post-processed to calculate expected values, variances, and the covariance of the regression parameters A₀ and b = A₀ λ where the distortion coefficient is extracted as λ_j = b_j/(A₀)_j. The pair of values ((A₀)_j, λ_j) is then further post-processed according to the GUM [1] for two sampled random variables q_j and r_j to calculate the covariance that is estimated as(12) (13) (14)

Investigations performed by Wuethrich and Souiyam [14] has concluded that a GLS regression produces high accuracy results, followed by the WLS regression which produces medium accuracy results, and finally by the OLS regression which produces low accuracy results. As a result the GLS approach may be considered the best overall approach however the key difficulty with implementing this approach is in the estimates of the covariances cov(S_i, S_j) that are required, which are either mathematically complex to determine or which must be numerically estimated with a large number of Monte Carlo simulations.

The key differences in the newer approach of Wuethrich and Souiyam [14] and an earlier approach of Ramnath [16] is that in the latter work a multivariate GUM Supplement 2 approach [17] was utilized that modelled the various inputs x by treating the two parameters A₀ and λ as a single vector y = [A₀, λ]^T with an equation h(y, x) = 0 where a multivariate PDF g_y(η_y) of the output was to be determined from a multivariate input g_x(ξ_x) by numerically solving the equation h(η_y, ξ_x) = 0 for a large number M of Monte Carlo simulation events.

In this approach sampled multivariate values of the inputs consistent with the joint PDFs of the inputs x were used to solve the homogeneous vector equation corresponding to and a sequence of pressures P₁, … , P_n and areas S₁, … , S_n were generated that was then used with an OLS approach to calculate sampled values of ((A₀)_j, λ_j) from the joint PDF g_y(η_y). The sampled values were then used to construct a bivariate copula that modelled u²(A₀), u²(λ) and cov(A₀, λ), however a limitation of this approach was that it required a large number of Monte Carlo simulations, was mathematically complex, and thus more suited to scientific metrology applications in national metrology institutes.

Differences in the two approaches is that firstly the method of Wuethrich and Souiyam [14] requires estimates for the uncertainties of the areas u²(S₁), … , u²(S_n) in the case of a WLS approach and additionally that of the covariances cov(S_i, S_j),  1i, jn for each Monte Carlo event j, and secondly makes an assumption of a Gaussian correlation between A₀ and λ. On the other hand the method of Ramnath [16] firstly avoids the explicit need for calculating u²(S_i) and cov(S_i, S_j) at each Monte Carlo event j as the data is considered as statistical samples, and secondly does not assume a Gaussian correlation between A₀ and λ since bivariate copulas can model non-Gaussian correlation effects.

It should be noted that the calculation of the covariances cov(S_i, S_j) that are required in a GLS regression whilst theoretically possible requires an additional matrix formulation solution as outlined in the GUM Supplement 2 [17], pp. 15–17. For an explicit matrix equation y = f(x) the covariance matrix of the output is and for an implicit matrix equation h(y, x) = 0 the covariance matrix of the output is solved from the matrix equation , where U_x is a covariance matrix of the input, and and are corresponding sensitivity or Jacobian matrices.

Due to the combination of algebraic and numerical ill-conditioning complexity that is present for determining a multivariate covariance U_y, this calculation is usually only tractable from a pure numerical based Monte Carlo simulation using an appropriate homogeneous equation for i = 1, … , n as previously discussed.

3 Mathematical modelling

We consider the standard case of a straight line equation y = ax + b with inputs x_i, outputs y_i, variances u²(x_i) and u²(y_i), and covariances cov(x_i, y_i) for measurements i = 1, … , n. This system was originally analysed by Krystek and Anton [10] who termed this system as a weighted total least squares with correlation system and developed a WTLSC algorithm that is able to calculate the variances u²(a) and u²(b) and the covariance cov(a, b). An earlier paper by Ramnath [18] utilized this scheme to quantify the variances u²(A₀) and u²(λ) but did not explicitly report on the covariance cov(A₀, λ) term.

Considering the straight line equation S = aP + b = A₀(1 + λ) it immediately follows that(15) (16)

From the above set of equations it follows that the variances are(17) (18)

Using a = A₀ λ it follows that

Expanding out the partial derivatives in the above equation and solving for cov(A₀, λ) then yields (19)

4 Numerical simulations

Numerical simulations are performed using the summary of statistical data reported in Table 1 which is based on an earlier full M = 1000 Monte Carlo simulation reported in detail in Ramnath [16]. In this table the columns are the values of P_k, u(P_k), S_k, u(S_k) and ρ(P_k, S_k) for k = 1, … , 10, where ρ(P_k, S_k) is the dimensionless correlation that is defined in terms of the covariance cov(P_k, S_k) by the formula(20)

Results are obtained by running the Octave/Matlab implementation of the WTLSC algorithm reported in Figure 1 with the script in Figure 2 as summarized in Table 2. All of the results reported in this table were directly obtained from a Octave/Matlab program with the exception of the dimensionless correlation coefficient ρ(A₀, λ) as the available number of Monte Carlo simulation events at M = 1000 was too small to obtain meaningful statistical results with the use of equation (11).

Instead the correlation was estimated by saving the processed data from the Octave/Matlab program as a txt file and fitting a Gaussian copula which is defined in terms of the covariance using the following RStudio computer code [19]:

library(copula)
library(VineCopula)
xydata <- read.table(''C:/mydata.txt'')
x1 <- xydata[, 1]
x2 <- xydata[, 2]
fncF1 <- ecdf(x1)
fncF2 <- ecdf(x2)
u1 <- fncF1(x1)
u2 <- fncF2(x2)
C21 <- BiCopSelect(u1, u2, familyset = c(1))
rho <- C21$par

This trick was utilized since a calculation of a covariance/correlation is mathematically equivalent to the assumption of a Gaussian bivariate copula to model the coupling effect between two random variables, and the VineCopula statistical package [20] uses more advanced optimization algorithms to estimate the parameters for a covariance/correlation. In general if a large number of Monte Carlo simulation events for example M = 10000 is available then equation (11) can be used directly, and this observation illustrates the essential need for a large number of simulations to adequately apply a GLS regression in order to obtain a statistically meaningful result particularly when a relatively small number of scattered data-points from a small number of Monte Carlo simulations is available for data analysis as shown in Figure 3.

It may be observed that in the case of a bivariate PDF that the joint PDF is then a three-dimensional function and as a result it becomes difficult to visualize the three PDFs g_OLS(x), g_WLS(x) and g_GLS(x) corresponding to the OLS, WLS and GLS regressions as these function overlap and intersect with one another. As a result in order to simplify the interpretation of the numerical results level curves of the PDFs are extracted as shown in Figure 4 so that the relative shape and sharpness of the PDFs may be meaningfully compared later in the paper.

Due to the complexity of the underlying pressure balance equations a large number of Monte Carlo simulations such as say M = 10000 is only feasible if the equations are implemented in an appropriate computer programming environment such as a Visual Basic for Applications (VBA) script within a Excel spreadsheet or a Matlab script, and as a result may present onerous technical computing demands on a pressure metrologist in an industrial calibration laboratory.

The results from Table 2 are visualized in Figure 5 from which it is confirmed that a OLS regression gives the lowest accuracy and a WLS gives a medium accuracy. The relative errors for the zero-pressure area , distortion coefficient , correlation , distortion coefficient uncertainty , and zero-pressure uncertainty , may be defined as (21)where Q represents a selection of a quantity and T denotes a selection of a technique.

Considering the relative errors for the WTLSC regression when compared to predictions from the GLS which is considered as the most accurate technique, it is observed that the relative errors are , , , , and respectively for the representative pressure balance data that was analysed. As a result it may be concluded that a WTLSC algorithm gives essentially the same quality of a results as full GLS Monte Carlo simulation if there are available estimates for the correlation terms cov(P_i, S_i) for the data points i = 1, … , n.

The key benefits of utilizing a WTLSC algorithm is that it is a relatively straightforward implementation of using a Octave/Matlab function file FitWTLSC with about 20 lines of computer code as outlined in Figure 2, and does not require extensive or advanced programming skills for a pressure metrologist to implement. Limitations of the method are nevertheless present in that estimates for the correlation terms ρ(P_i, S_i),  i = 1, … , n are technically necessary to utilize the WTLSC algorithm.

This issue may be addressed in one of two practical ways, namely (i) setting zero correlations for all of the cross-float data-points such that ρ(P_i, S_i) = 0 for i = 1, … , n for which the WTLSC algorithm will reproduce results for a WLS algorithm as a conservative compromise, or (ii) performing a small number of say 50 or 100 repeat measurements in an Excel spreadsheet which is feasible without the need for any specialist VBA scripting programming within a spreadsheet with appropriate expert physical judgement to generate (P_i, S_i) data to fit a Gaussian copula with the 11 lines of RStudio code previously discussed to obtain an approximate estimate for the correlation.

Referring to the results in Figure 3 it is observed that the PDFs of the OLS, WLS and GLS regressions may be compared by examining the behaviour of the PDFs along the semi-major and semi-minor axes of the hyper-ellipsoid for the confidence region. When this analysis was performed as visualized in Figure 4 it was observed that the WLS regression performed relatively well when compared to the GLS regression. As a result it may be concluded that the limitations of the potential absence of a covariance matrix is not excessively detrimental.

Fig. 1 Gnu Octave/Matlab function file to implement weighted total least squares with correlation algorithm. The x input and its standard uncertainty u(x) have the same physical units, and similarly the y output and its standard uncertainty u(y) also have the same physical units, whilst the correlation ρ(x, y) = cov(x, y)/(u(x) × u(y)) is dimensionless.

Fig. 2 Gnu Octave/Matlab script file using the function file FitWTLSC.m to calculate the pressure balance parameter values and associated variance and covariance information in a typical pressure laboratory calibration certificate for a client pressure balance instrument. The data in the loaded file is from Table 1 and may conveniently be extracted from an Excel spreadsheet of measurement results and as saved a txt file.

Fig. 3 Comparison of ordinary least squares (OLS), weighted least squares (WLS) and generalized least squares regression analysis from Monte Carlo data. Hyper-ellipsoids are generated to indicate the corresponding 95% confidence regions (CRs) and reveals that the highest accuracy scheme is that of a GLS which has the smallest CR which corresponds to a high degree of certainty of the range of possible values, followed by a WLS which has a slightly lower accuracy compared to a GLS, and finally a OLS which has the lowest accuracy as the CR is unduly large.

Fig. 4 Comparison of joint PDFs from OLS, WLS and GLS regressions. The confidence region may be visualized as a three dimensional volume underneath the envelope covered by the joint PDF and bounded by the hyper-ellipsoids, and indicates that the joint PDFs may fit inside one another. Due to the difficulties with visualizing three different joint PDFs intersecting within each other simultaneously, cuts along the major axis (orange line in previous figure) and minor axis (magenta line in previous figure) of the ellipsoids confidence regions are taken which produce level curves to more easily visualize the shape, orientation and size of the different OLS, WLS and GLS regression schemes. From these level curves it is observed that for level curves along the ellipsoid major axis that the GLS regression (green curve) has the narrowest/sharpest curve which corresponds to the highest confidence or accuracy in the estimates of the parameters, then the WLS regression (blue curve) which has a medium spread/sharpness, and finally the OLS regression (red curve) which is relative wide and flat which corresponds to less certainty/accuracy of the parameters in the regression analysis. A broadly similar behaviour is observed for level curves cut along the minor ellipsoid axis, and as a result it is observed and concluded a OLS regression has the least accuracy, a WLS regression has a medium accuracy, and a GLS regression has the highest accuracy. From the above graphs it may be concluded that a straight line regression for pressure metrology work should ideally utilize a GLS regression (green line) however if no covariance information is available the a WLS regression (blue line) can give a reasonable estimate particularly for industrial pressure metrology laboratories.

Fig. 5 Comparison of relative errors of the quantities predicted by OLS, WLS and WTLSC regression approaches relative to predictions from a GLS regression approach. These graphs indicate that the relative error between a WTLSC regression and a GLS regression is relatively insignificant, and as a result a pressure metrologist instead of applying a more complicated GLS regression can simply apply a WTLSC regression with the aid of the supplied methodology developed in this paper.

5 Conclusions

Based on the research reported in this paper the following conclusions were determined:

• Numerical simulations have been performed that demonstrate that a weighted total least squares with correlation (WTLSC) algorithm gives essentially equivalent results for straight line parameters as a more advanced generalized least squares (GLS) algorithm with additional Monte Carlo simulations for parameter uncertainty estimation

• A method has been developed for the WTLSC algorithm that is able to analytically extract the variance/covariance information of the straight line parameters without the need for additional and more advanced Monte Carlo simulations

6 Influences and implications

Based on the research reported in this paper the following influences and implications are:

• Pressure metrologists working in industrial calibration laboratories may now utilize the WTLSC algorithm to conveniently and easily quantify the expected values and variances/covariance of the zero-pressure area A₀ and distortion coefficient λ

• A computer program FitWTLSC.m is now freely and publicly available to conveniently implement the WTLSC algorithm in Gnu Octave or Matlab for calculating the straight line parameters and their associated variance/covariance information and may now also be used by other metrologists from different fields

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 Hyper-ellipsoidal bivariate confidence region

Let random variables x₁, … x_M ∈ ℝ^N be grouped as Ω = [x₁| … |x_M] where x_j,  j = 1, … , M are N × 1 vectors, Ω is a N × M matrix, and M is a suitably large positive integer number of Monte Carlo simulation events. Then the expected value μ a N × 1 vector and the covariance matrix U a N × N matrix may be estimated with the GUM Supplement 2 [17] as μ and U. Under these circumstances the joint PDF, unless higher order statistics information over and above the expected and covariance values are provided, is then assumed to follow a multivariate Gaussian distribution g(x) ~ (μ, U).

The corresponding joint PDF is then g(x) and consists of hyper-ellipsoids within a ℝ ^N dimensional space that are centred at μ and have an orientation determined from characteristics of U such that(A1) (A2) (A3)

In the special case of bivariate data where , and , the hyper-ellipsoid will lie in the x₁ x₂-plane within a 2 dimensional space, and the joint PDF will be (A4)

For the above equation the terms are , , and since the covariance matrix may equivalently be written as . According to the GUM Supplement 2 [17] the 100p% coverage region is a hyper-ellipsoid of the form(A5)

An earlier analysis by Ramnath [21] determined that the ellipsoid semi-axis values a₁ and a₂ could be approximately determined by forming the matrix(A6)

Once this matrix was formed the two eigenvectors v₁ and v₂ along with the corresponding eigenvalues λ₁ and λ₂ were calculated from the equation(A7)

The ellipse semi-axis values were then calculated as(A8) (A9)and the eigenvectors v₁ and v₂ then define the orientations of the ellipsoid principal axes following the analysis by Lee [22]. The disadvantage of this approach is that an appropriate software routine for determining eigenvectors and eigenvalues is necessary, however an advantage of this approach is that it is generalisable to higher dimensions.

An alternative direct algebraic approach is to directly substitute the statistical information μ and U into the formula equation (A5) and expand out so that

(A10)

The above equation is in the form of the generalized ellipse equation as discussed by Weisstein [23] that is centred at x_o and y_o, with semi-major and semi-minor axes a and b and rotated by an angle θ of the form(A11)

The algebraic coefficients may now be expressed in terms of the statistical terms as(A12) (A13) (A14) (A15) (A16) (A17)

The corresponding ellipse geometry parameters may now be algebraically calculated as (A18) (A19) (A20)

and the ellipse rotation is given as(A21)

The hyper-ellipsoid may then be conveniently plotted using the parametrization 0 ≤ t ≤ 2π so that the initial boundary points x_i and y_i are(A22) (A23)

Next the boundary points are rotated by the angle θ by using a rotation matrix to give x_r and y_r such that(A24)

Finally the points are translated as(A25) (A26)

This algorithm is used to produce the hyper-ellipsoids in Figure 3 and may be used to determine the points along the ellipse major and minor axes. As a result if these lines are used a “cut” through the joint PDF may be used to study how the level curves of the different OLS, WLS, and GLS regressions behave as it is too difficult to easily visualize the three dimensional joint PDFs in Figure 4 even with transparency effects. Once the lines are determine the x₁ and x₂ points are known from a simple straight line fit so the x₁ and x₂ points may then be substituted in the bivariate joint PDF in equation (A4) which will produce a density value. The set of points (x₁, x₂, g(x₁, x₂)) may then be plotted as a three-dimensional curve in order to visualize the shape and sharpness of the corresponding level curve of the PDF. This algebraic approach is only suitable for bivariate simulations and becomes impractical in higher dimensions, for which the eigenvector/eigenvalue approach is then necessary.

Appendix B Simplified worked example

In order to illustrate the method suppose that the experimental data has n = 5 points of generated pressures from 100  MPa through to 500  MPa in 100  MPa steps and cross-floated areas, which will give an experimental control variable of x = [x₁, x₂, x₃, x₄, x₅] and an experimental observed variable of y = [y₁, y₂, y₃, y₄, y₅]. Writing down the basis equation with the experimental data (for one particular Monte Carlo event) as follows.

Let

Considering an ordinary least squares (OLS) regression the estimate for the matrix V is justwhere W = V⁻¹.

Implementing the formula in Matlab or GNU Octave with the above information gives

Next considering a WLS regression it is known that the V_i,i ∝ u²(S_i) so using the data from Table 1 let u(S₁) =0.0000990361, u(S₂) = 0.0001110578, u(S₃) = 0.0001766929, u(S₄) = 0.0002740451 and u(S₅) = 0.0004087028, and then simply let

.

Using Matlab/Octave to implement the matrix operations then produces

The above process is repeated for all of the experimental data points and then the results can be post-processed in order to obtain the statistical characteristics for the OLS, WLS and GLS regressions as shown in Figure 3.

Actual physical data is first obtained from the experimental measurements in the laboratory and this data is used as inputs into the instrument models to produce the data X and Y in this example where the computer program outlines the entire logical calculation process. As an example the piston temperature may have been measured as say 20.4^∘C with a standard uncertainty of ±0.1^∘C with an assumed Gaussian distribution, then random temperatures consistent with the physical experimental uncertainty are sampled and used as inputs into the physical model. This is due to the fact that modern measurement science has generally adopted a Bayesian statistics framework for measurement uncertainty and as a result no physical measurement is deterministic since all physical measurements, with the partial exception of certain physical constants, have an underlying probability distribution i.e. any physical measurement is considered as a statistical sample from that random variables probability density function. Once sufficient samples have been generated and used as inputs/outputs into the measurement model then the OLS, WLS or GLS regressions may be implemented as outlined in Figure B1.

Fig. B1 Gnu Octave/Matlab function file to implement OLS, WLS and GLS regression analysis.

The general steps to perform the regression analysis are therefore as follows:

• Obtain laboratory nominal experimental measurement data x that has an associated joint PDF g_X(x)

• Sample a random variable ξ from the joint PDF g_X(x)

• Solve the associated measurement equation h(η, ξ) = 0 to obtain the random variable η associated with ξ, and repeat this several times to produce the Monte Carlo data Ω

• Step through the data input ξ and data output η i.e. the inputs x and outputs y and substitute this in the regression equation C = (X^T W X)⁻¹(X^T W Y)

• Solve the regression equation for C and then post-process all the results to work out the statistical characteristics of A₀ and b namely u²(A₀), u²(b) and cov(A₀, b)

• Simply apply the new algebraic formula to work out u²(A₀) = u²(b), , and

The above computer code routines for performing the statistical regressions may be observed to be relatively lengthy and complicated, and may be compared to the WTLSC implementation which is simply the 20 lines of Matlab/Octave code in Figure 2. In this approach there is no need to perform any Monte Carlo simulations and the actual physical experimental data, in this simplified problem the n = 5 data points, may be used directly as inputs into the Octave/Matlab program FitWTLSC.m to directly extract the uncertainty information of the parameters A₀ and λ and obtain u²(A₀), u⁽λ) and cov(A₀, λ).

References

All Tables

All Figures

	Fig. 1 Gnu Octave/Matlab function file to implement weighted total least squares with correlation algorithm. The x input and its standard uncertainty u(x) have the same physical units, and similarly the y output and its standard uncertainty u(y) also have the same physical units, whilst the correlation ρ(x, y) = cov(x, y)/(u(x) × u(y)) is dimensionless.
In the text

Fig. 2 Gnu Octave/Matlab script file using the function file FitWTLSC.m to calculate the pressure balance parameter values and associated variance and covariance information in a typical pressure laboratory calibration certificate for a client pressure balance instrument. The data in the loaded file is from Table 1 and may conveniently be extracted from an Excel spreadsheet of measurement results and as saved a txt file.

In the text

Fig. 3 Comparison of ordinary least squares (OLS), weighted least squares (WLS) and generalized least squares regression analysis from Monte Carlo data. Hyper-ellipsoids are generated to indicate the corresponding 95% confidence regions (CRs) and reveals that the highest accuracy scheme is that of a GLS which has the smallest CR which corresponds to a high degree of certainty of the range of possible values, followed by a WLS which has a slightly lower accuracy compared to a GLS, and finally a OLS which has the lowest accuracy as the CR is unduly large.

In the text

Fig. 4 Comparison of joint PDFs from OLS, WLS and GLS regressions. The confidence region may be visualized as a three dimensional volume underneath the envelope covered by the joint PDF and bounded by the hyper-ellipsoids, and indicates that the joint PDFs may fit inside one another. Due to the difficulties with visualizing three different joint PDFs intersecting within each other simultaneously, cuts along the major axis (orange line in previous figure) and minor axis (magenta line in previous figure) of the ellipsoids confidence regions are taken which produce level curves to more easily visualize the shape, orientation and size of the different OLS, WLS and GLS regression schemes. From these level curves it is observed that for level curves along the ellipsoid major axis that the GLS regression (green curve) has the narrowest/sharpest curve which corresponds to the highest confidence or accuracy in the estimates of the parameters, then the WLS regression (blue curve) which has a medium spread/sharpness, and finally the OLS regression (red curve) which is relative wide and flat which corresponds to less certainty/accuracy of the parameters in the regression analysis. A broadly similar behaviour is observed for level curves cut along the minor ellipsoid axis, and as a result it is observed and concluded a OLS regression has the least accuracy, a WLS regression has a medium accuracy, and a GLS regression has the highest accuracy. From the above graphs it may be concluded that a straight line regression for pressure metrology work should ideally utilize a GLS regression (green line) however if no covariance information is available the a WLS regression (blue line) can give a reasonable estimate particularly for industrial pressure metrology laboratories.

In the text

Fig. 5 Comparison of relative errors of the quantities predicted by OLS, WLS and WTLSC regression approaches relative to predictions from a GLS regression approach. These graphs indicate that the relative error between a WTLSC regression and a GLS regression is relatively insignificant, and as a result a pressure metrologist instead of applying a more complicated GLS regression can simply apply a WTLSC regression with the aid of the supplied methodology developed in this paper.

In the text

	Fig. B1 Gnu Octave/Matlab function file to implement OLS, WLS and GLS regression analysis.
In the text