Interpolating between discrete calibration points: least squares curve fitting

ABSTRACT
In preceding publications, the accuracy of interpolation using reference functions for temperature sensors (thermocouples and platinum resistance thermometers (PRTs)), as well as a crude approach to interpolation uncertainty in the absence of any knowledge of the interpolating function, have been discussed. This article discusses how to fit curves to PRT calibration results using the linear least squares technique (implemented using matrix functions in Microsoft Excel or Libreoffice Calc), both directly to the measured data and to deviations from a reference function. (An overview of the fitting method may be found in [Numerical Recipes in Fortran 77, Chapter 15. Modeling of data].)

INTRODUCTION

In the preceding articles, we saw that
(i) industrial PRT calibration results could be interpolated “piece-wise” (between any two calibration points a few hundred degrees Celsius apart) to an accuracy of around 0.05 °C, when expressed as deviations from the Callendar-van Dusen reference functions, and,
(ii) without any knowledge of the accuracy of the interpolating function, the accuracy of interpolated values could be grossly estimated as 0.15 to 0.5 °C (for the above-mentioned spacing between calibration points).
While the latter approach can be applied generically to almost any calibration data, it has the drawback that the uncertainty of interpolated values is unrealistically large when working close to a calibration point. In the present article, the most accurate approach to the problem will be implemented, namely, to fit a curve to the complete set of calibration results. Firstly, a Callendar equation will be fitted directly to the measured data, and, secondly, a quadratic polynomial will be fitted to the deviations of the measured data from the ITS-90 PRT reference function. As the ITS-90 reference function deals with most PRT non-linearity, it is expected that the latter approach will produce the more accurate interpolation.

THE MODEL EQUATION

The Callendar equation, applicable above 0 °C, is: Rt = R0∙(1 + A∙t + B∙t^2), where Rt is the resistance (in Ω) at temperature t (in °C) and R0 is the resistance (in Ω) at the ice point (0.00 °C). Rearranging the equation, A∙t + B∙t^2 = Rt/R0 – 1, indicating that the independent variable x = t and the dependent variable y = Rt/R0 – 1.
(As the sensitivity. d(Rt/R0)/dt, will be required to convert uncertainties from temperature units, we also note that d(Rt/R0)/dt = A + 2∙B∙t. We will use the coefficients from IEC 60751, namely, A = 3.9083e-3 and B = -5.775e-7, to calculate sensitivities below.)

The measured data are as follows:

IPRT calibration data.
IPRT calibration data.

To find R0 from R(0.01 °C): R0 = R(0.01 °C) + (0.00 °C – 0.01 °C) ∙ 0.391 Ω/°C. Uncertainties are converted to the same units as y. (The 1st data point is used to determine R0, so only the four subsequent data points are used in the fit.)
IPRT data converted to appropriate units.
IPRT data converted to appropriate units.

The general model equation is c1∙f1(x) + c2∙f2(x) + … = y. For the Callendar equation, c1 = A, f1(x) = x, c2 = B and f2(x) = x^2.
The four data points lead to four simultaneous equations, namely:
c1∙f1(x1) + c2∙f2(x1) = y1
c1∙f1(x2) + c2∙f2(x2) = y2
c1∙f1(x3) + c2∙f2(x3) = y3
c1∙f1(x4) + c2∙f2(x4) = y4
Expressed in matrix notation, they are:
The unweighted simultaneous equations, expressed in matrix notation, are: A x c = y.
The unweighted simultaneous equations, expressed in matrix notation, are: A x c = y.

WEIGHTING FACTORS

Now, apply the weighting factors 1/ui to each equation (smaller std uncertainty => larger weight). The weighted simultaneous equations, expressed in matrix notation, are: W x A x c = W x y.
Weighting_matrix_W

THE SOLUTION – NORMAL EQUATIONS

The optimal (least squares) solution to this over-determined system is obtained by left-multiplying by the matrix transpose (WA)T = ATxWT, to obtain the normal equations:
Normal_equations
Matrix WTxW simply contains 1/variance at each calibration point:
Matrix_WtW
Coeff_c
Then, to find the fitted y-value for any value of x:
y_for_x
Variance_of_y
Finally, expanded uncertainty of y = 2*√variance(y): this is called the “propagated uncertainty” of y at x.

The numerical implementation of this “linear least squares” curve-fitting technique follows, using the Excel or Calc matrix functions TRANSPOSE(), MMULT() and MINVERSE():
Numerical_implementation
Now to check the fitted curve against the measured data (residual = measured – fitted):
Residuals
Are the residuals small enough, or, does the model fit the data well enough, relative to the uncertainties?: The chi-squared statistic, chi^2 = [(residual_1/u1)^2 + (residual_2/u2)^2 + …], will tell us:
Chi-squared
Chi-squared is larger than the degrees of freedom (= number of data points – number of fitted parameters = 4-2), so either the uncertainties are underestimated, or the model does not represent the data well (relative to the uncertainties). If chi^2 ≲ d.o.f., then the fit is “good enough”. (This is only strictly true if the uncertainties at different points are uncorrelated, but we will assume that this is the case.)

DEVIATION FROM ITS-90 REFERENCE FUNCTION

Perhaps direct fitting to the Callendar equation is not good enough, at this level of uncertainty. Let’s try deviations from the ITS-90 reference function – the reference function should take care of much of the PRT’s non-linearity, leaving us with smoother data to fit.

For the range 0 to 660 °C, the ITS-90 document gives the deviation function as: W-Wr = a∙(W-1) + b∙(W-1)^2 + c∙(W-1)^3, where W = Rt/Rtp and Wr is the value of the reference function.
It is good practice in curve fitting to use as few coefficients as will adequately fit the data, so we will try two, a and b.
ITS-90_dev_data
Dev_implementation
Now to check the fitted curve against the measured data (residual = measured – fitted):
Resid_dev
It can be seen that this fit is better than the previous one. (Residuals are around half the size.)
Chi2_dev
Chi-squared is smaller than the degrees of freedom, indicating that the curve fits the data well, within the uncertainties.

Below is a table of fitted values. Because of the structure of the ITS-90 functions, the uncertainty of the curve is zero at the water triple point (WTP, 0.01 °C). The uncertainty at WTP is added to the uncertainty of the curve, in the rightmost column.
Fitted_table

CONCLUSIONS

∙Both curves, fitted directly to measured PRT data, and to deviations from a reference function, represent the behaviour of the UUT better than do the individual data points. (The curves take into account all the data points, thereby potentially “smoothing out” random errors in measurement.)
∙The ITS-90 deviation function fits the data better, as is expected when using a good reference function.
∙When weighted least squares fitting is performed, the covariance matrix provides a statistically justified manner of propagating uncertainty to intermediate points, which results in small uncertainties. (To be strictly correct, we should have considered possible correlation between data points, but, as long as the dominant uncertainty component(s) are uncorrelated, our approach is reasonable.)

(Contact the author at LMC-Solutions.co.za.)

Interpolating between discrete calibration points: the effect on uncertainty – Addendum

ABSTRACT
This publication follows “Interpolating between discrete calibration points: the effect on uncertainty” of September 2025, describing a more general approach to interpolation of calibration uncertainty (not specific to a particular type of sensor), drawing a tentative conclusion regarding the spacing of calibration points versus resultant uncertainty, and applying this approach to digital thermometer, PRT and thermocouple calibration data.

INTERPOLATING UNCERTAINTY “BY EYE”

Here are two sets of calibration results, reporting correction or error of the Unit Under Test (UUT), and expanded uncertainty, at various temperatures:

Calibration results of a digital thermometer with type K thermocouple sensor, and an industrial PRT.
Calibration results of a digital thermometer with type K thermocouple sensor, and an industrial PRT.

Looking at the left-hand (digital thermometer) results, what do we observe?:
1. The correction at all three temperatures is effectively constant, relative to the calibration uncertainty.
2. The UUT is a thermocouple thermometer with a range of -50 to 300 °C. So, none of the calibration temperatures appears to be “special”. (In this context, “special” temperatures are those where the thermometer might be adjusted to have small corrections, for example, the ends of the operating range, room temperature (where all signal comes from Cold Junction Compensation and none from the sensor), and temperatures where the measuring electronics change range. For a liquid-in-glass thermometer, “special” temperatures would be scale pointing marks.)
These observations lead us to a (tentative) conclusion: The thermometer correction is stable at three “random” temperatures, so the correction at intermediate temperatures can probably be estimated with confidence, without any increase in uncertainty.

Now, looking at the right-hand (IPRT) results, we see:
1. The error varies significantly between calibration points, relative to the uncertainty.
2. The error varies fairly linearly with temperature, though deviation from a straight line is sometimes larger than the calibration uncertainty (at 232 °C, in this case):

Calibration results of the above industrial PRT, graphed.
Calibration results of the above industrial PRT, graphed.

What may we conclude from the IPRT data?: If we wish to estimate the error at an intermediate temperature by linear interpolation, the uncertainty at this intermediate temperature should probably be larger than that at the neighbouring calibration temperatures.

INTERPOLATING UNCERTAINTY – NUMERICAL ESTIMATE

How may we numerically estimate the additional uncertainty caused by interpolation?: If we assume that the value at the interpolated point lies between the two neighbouring calibration values, with an equal probability anywhere in that range, the uncertainty caused by interpolation may be estimated as |corr_1 – corr_2|, as the full-width of a rectangular distribution. (Note: The assumption that the value lies between the two neighbouring calibration values is not necessarily correct for “badly behaved” instruments.) To obtain standard uncertainty, u_interp, divide this by 2√3. Or, for expanded uncertainty, U_interp = u_interp*2 = |corr_1 – corr_2| / √3. Total uncertainty U_tot = √[U_cal^2 + U_interp^2].
This technique is applied for the digital thermometer (using a more recent, larger, data set) and the IPRT mentioned above, with values interpolated between low, middle and high temperatures (in bold) being compared to measured values at intermediate temperatures:

Linear interpolation between bold calibration results, for digital thermometer and IPRT.
Linear interpolation between bold calibration results, for digital thermometer and IPRT.

Linear interpolation between the bold calibration results differs from the measured corrections or errors at intermediate temperatures (residual = measured – interpolated) by less than the estimated total uncertainty, suggesting that this approach to interpolating uncertainty is “safe”, at least for these two instruments.
For the digital thermometer, the uncertainty of interpolated values is essentially the same as that of neighbouring calibration points (differs less than 5%), while interpolated values have significantly (~ten times) larger uncertainty for the IPRT. This seems like a reasonable approach, without having any deeper knowledge of the interpolating function (reference function) being used. (The only assumption is that the correction or error varies “more-or-less” monotonically between calibration points.) Note that interpolation may be much more accurate if one does have deeper knowledge of the interpolating function, as reported in the preceding paper, or if one performs a least squares fit to the full set of data.

The digital thermometer results above were stable, and those of the IPRT close to linear. What about non-linear results? The tables and graphs below show results for a type R thermocouple (linear) and a type K thermocouple (non-linear):

Linear interpolation between bold calibration results, for type R and type K thermocouples.
Linear interpolation between bold calibration results, for type R and type K thermocouples.

Calibration results for type R and type K thermocouples, graphed.
Calibration results for type R and type K thermocouples, graphed.

It can be seen that, though the type K’s results are non-linear (and even slightly non-monotonic) at low temperatures, the residual (= measured – interpolated) is always smaller than the total uncertainty.

SPACING BETWEEN CALIBRATION POINTS

If one wishes to add negligible uncertainty from interpolation, how far apart may the calibration points be? Consider the following “generic” calibration results:

Generic calibration results; correction or error y vs independent variable x.
Generic calibration results; correction or error y vs independent variable x.

If the correction or error, y, changes by at most half the expanded uncertainty, between one calibration point and the next, then the total uncertainty at interpolated points is essentially the same as at the calibration points (differs less than 5%).

CONCLUSION

∙The correction of a UUT at a point intermediate between calibration points may be estimated in a simple manner, by interpolating linearly between the two adjacent calibration values. The uncertainty at this interpolated point may be grossly estimated by adding the difference between the two calibration values in quadrature to the calibration uncertainty. For UUTs with stable corrections, this adds negligibly to the uncertainty, but if the correction varies significantly between calibration points, the resultant uncertainty at the intermediate point is much larger.
∙If the user specifies the calibration points, he takes responsibility for the uncertainty between points [Petersen, “Principles for Calibration Point Selection”, NCSLI Measure, Volume 8, No 3, 2013].
∙If the user only specifies the range, the calibration laboratory should suggest calibration points based on
– instrument manufacturer’s recommendation
– understanding of the operating principles of the instrument
– historical experience with this type or model of instrument (“type testing”)
– in the absence of further information, sufficient points that the UUT correction changes by at most half the expanded uncertainty, between one calibration point and the next (unless some uncertainty is added for interpolation)
– ideally ,the calibration laboratory should be confident enough to include a statement in the certificate such as “Effects due to interpolation are considered negligible over the calibrated range.”
∙The user should agree to the suggested calibration points during contract review.

(Contact the author at LMC-Solutions.co.za.)