M-estimators based on Huber and Andrews sine loss functions were successfully used for approximation of heat capacities and heat contents of K-substituted natrolite and petalite by means of the weighted sum of Einstein functions. It automatically excluded outliers for petalite and narrow peak of lambda-transition for K-natrolite. k=1 k=1 where n is the number of points, rk are residuals, ycalc and yexp are calculated k k and experimental values respectively, β is the model parameters column vector, ρ(t) is the loss function, σ is the scaling factor, ωk are statistical weights. However, minimization of eq. 1 requires special algorithms and much more computational power than the least squares method. The aim of this work is to demonstrate the applicability of M-estimators for approximation of heat capacities and heat contents of individual substances. Experimental data for petalite LiAlSi4O10 and K-substituted natrolite (K-natrolite) ⋆ Supported by the Russian Foundation for Basic Research for providing financial support (Grant No. 20-03-00575) and by the “Chemical Thermodynamics and Theoretical Materials Science” program (No. 121031300039-1).
Evaluation of parameters of thermodynamic models from experimental data is
a very common problem of nonlinear optimization. It is usually based on the
weighted non-linear least squares method. Selection of the statistical weights is
a complex problem due to diferent accuracy of experimental data and possible
presence of systematic errors and outliers. Diferent schemes of their automatic
selection were suggested [
F (β ) = n X ρ
ycalc(β ) − yexp ωk k k Na0.01K1.85Mg0.01Ca0.04Al1.96Si3.04O10 · 2.72H2O were be used as examples. (1)
In this work the the iterative reweighted least squares (IRLS) algorithm
combined with Levenberg-Marquardt type regularization technique was used for
finding eq. 1 minimum. IRLS for robust regression was suggested by Mudrov, Kushko
et al. [
n = X ψ k=1 rkωk σ
ωk ∂rk = 0; ψ (t) ≡ ρ(˙t) ≡
· σ ∂β i
∂ρ(t)
∂t
IRLS uses two simplifications to get rid of the second derivatives. The first one
is linearisation of the deviations rk near the initial approximation β ◦ :
rk(β ) = rk(β ◦ ) +
m
X ∂rk (β j − β j◦ ) ⇒ r = r◦ + J p; p = β − β ◦
j=1 ∂β j
where J is Jacobian (n × m matrix), m is the number of parameters, r is the
residuals column vector, The second step is exclusion of the ρ(¨t) function by
introduction of so called weight function w(t) = ψ (t)/t [
J ⊤W r◦ ; C = (r◦ )⊤W r◦ n − m (J W J )−1 where W is the n × n diagnonal matrix with the Wkk = ωk2 · w(ωkrk/σ ) elements.
The IRLS algorithm was embedded into the CpFit program [
σ = Φ −1 (0.75) · median |r| = 1.483 · median |r|
where Φ −1 (x) is the inverse cumulative distribution function for standard
normal distribution.
3. Run the IRLS iterations combined with Levenberg-Marquardt type
regularization technique using the given data and the loss function ρ(t).
(2)
(3)
(4)
(5)
(6)
The ρ(t) = 0 .5t2 (i.e. the least squares method), Huber and Andrews sine loss
function were used in this work, their ρ(t), ψ (t) and w(t) are given in Table 1
and at Figure 1. The tuning constants a for 95% asymptotic eficiency in the
case of normal distribution of errors were taken from Holland and Welsch [
Andrews sine and Huber functions are piecewise. They turn into AT 2 at smaller t and to constant and linear functions respectively at larger t. This reduces values of their weight functions w(t) at larger t and influence of outliers. For Andrews sine function w(t) reaches 0 for finite values of t. It causes exclusion of potential outliers from the optimization. In the case of Huber function w(t) is always positive.
Huber function is convex and Andrews sine function is not (see Figure 1).
The latter one belongs to redescending M-estimators that have non-convex ρ(t),
ψ (t) with local extrema and limt→∞ ψ (t) = 0. They allow to totally exclude
outliers from the optimization but may lead to non-convex objective function
(see eq. 1) even in the case of linear regression. This increases the possibility of
reaching local minimum instead of global and requires more careful selection of
the initial approximation [
Experimental data for K-substituted natrolite and petalite heat capacity and
heat content were considered in this work. They are summarized in Table 2.
These data were already approximated earlier by Voskov et al. [
0 HT −
H0 =
Cp(T ) dT = θ i T ;
CE(x)
R m X α iHE i=1 θ i T = ;
3x2ex (ex − 1)2
HE(x)
RT =
3x
ex − 1
(7)
(8)
where m is the number of terms, R is the universal gas constant, CE(x) is Einstein
function, α i and θ i are model parameters that are found by the minimization of
eq. 1. They may be considered as a crude approximation of phonon spectrum
but due to anharmonism and possible Schottky anomalies they are closer to
adhoc parameters. However the approximation based on the least squares method
required manual exclusion of low-temperature outliers for petalite [
The results of approximation for both substances are shown at Figure 2. Higher uncertainties at T < 25 K are due to less accurate experimental data, see Table 2. The corresponding α i and θ i values are given in Tables 3 and 4. Extra digits are left intentionally: parameters confidence intervals because parameters are correlated to each other, i.e. C from eq. 5 is not diagonal. This is typical for linear and nonlinear regression. Number of terms for petalite is not equivalent for diferent models because attempts to increase number of terms up to 5 in all models caused ill conditioned optimization problems.
80 60 40 20 0 -20 0
The least squares method is sensitive to the Cp anomaly and outliers: the
obtained models are not accurate and undergo oscillations. M-estimators based
on Huber and Andrews sine function are much less sensitive to them. Standard
“baseline” (i.e. not taking into account Cp anomalies) entropies S2◦,B98L.15 were
calculated for all models, see Table 5. For K-natrolite uncertainties are 1.1%, 0.5%
and 0.3% for quadratic, Huber and Andrews loss functions; the reference value
S2◦,B98L.15 = 437.7 J · (mol · K)−1 was taken from [
Although parameters confidence intervals were controlled to avoid overfitting, k-fold cross-validation with k = 5 was made for all models. The results are present in Table 5, all standard errors were estimated by means of eq. 6. s were calculated for parameters from Tables 3 and 4 before cross-validation. stest and strain were evaluated as mean standard errors for test and training sets respectively; stest and strain were estimated by means of eq. 6.
For Andrews sine functions obtained sCp and s∆H are close to the standard
errors of existing models: for K-natrolite sCp = 0.68% (restored from model
parameters from [
Robust regression based on M-estimators and the IRLS algorithm were successfully applied for approximation of isothermal heat capacity and heat content of K-substituted natrolite and petalite. This approach allowed to automatically exclude outliers. However, it is not designed for estimation of random and systematic errors of diferent data series, and may be combined with other schemes of statistical weights assignment if required. It also can’t replace critical data evaluation of available experimental data but may help to find anomalies and outliers. 6
CpFit program is available at the site of Laboratory of Chemical Thermodynamics (http://td.chem.msu.ru). Data files for K-natrolite and petalite are published as Mendeley data set (http://dx.doi.org/10.17632/gbgnkr3f2x.1).