Parameters to be found

Criterion function. Robust estimates

Algorithms

Stop criteria

Residual variance and criterion c²

Skewness, kurtosis, mean residual and residual mean

Estimation of errors and correlation coefficients for parameters

Redundancy of a model as decided from the analysis of Jacoby matrix

To the Help main page

1. logarithms p of unknown equilibrium constants;
2. intensity factors W of reagents (e.g. molar absorptivities in spectrophotometry) for L analytical positions (e.g. wavelength in spectrophotometry).

Total number of parameters to be found z = p + L ґ W.

Here lg b_i⁰ is the invariable contribution into lg b_i;

t_iu are elements of a pґ p matrix T transforming parameters LK_u into lg b_i.

See description of input data

For simplicity at this section we shall assume lg b_i⁰ = 0 and matrix T to be the unit one. Hence, lg b_i are considered as parameters to be calculated.

A vector q (composed of unknown b _i, a _li) that gives minimum to the chosen criterion function U:

|q > = arg min U(q)

is adopted as a proper estimation of sought-for parameters.

If all measured A_lk are independent and the density of their errors e is known, then the criterion function is chosen according to maximum likelihood principle.

If the density is normal, parameters may be estimated using the nonlinear least-square method as vector q that gives minimum to functional

where w_lk is the statistical weight corresponding to variance estimate s²(D_lk),

D_lk is the discrepancy

A_lk is a property of the equilibrium system.

If the distribution of e deviates from the normal one, estimates obtained using the least squares method (LSM-estimates) are biased, insignificant and inefficient. Alternative ways of estimating parameters that are less sensitive to the violation of the LSM basic assumptions are termed robust.

In developing a robust criterion for parameter identification a distribution of e is assumed with "longer" tails in comparison with the normal distribution (i.e. a fractional mass is moved from the central region to the tails). These tails increase the probability for e considerably deviated from zero. A quantitative measure of tail length is the kurtosis g₂: g₂=0 for the normal distribution and g ₂>0 if tails are longer.

Let measurement errors e have density

р(e) = [(100 - d)Ч j(e) + d Ч h(e)] / 100,

where j (e) is the standard density with zero mean and unity variance (density of main observations), h(e) is the density of outliers, d is the intensity of outliers (%). This hypothesis about experimental errors is termed "the model of outliers (gross errors)". It was shown by Huber [Huber P.J. Robust Statistical Procedures, SIAM, 1996] that the maximum likelihood estimate corresponds in this case to the minimum of criterion

where x _k is the weighted discrepancy (w_k^1/2Ч D_k), r(x_k, q) is a loss function. This formula is given for the case of one analytical position, L =1.

Estimates corresponding to the minimum of the above functional are termed M-estimates. They are non-biased, significant and highly efficient provided that function r(x_k, q) increases slowly than x². M-estimates are more robust with respect to outliers than LSM-estimates. The loss function introduced by Huber has the form

where c is a constant. It depends on the intensity of outliers d and is found as a solution to equation

M-estimates keep their optimal properties in the case of small numbers of experiments and also robust to asymmetrical contaminations of the normal distribution of errors.

Criterion M is a hybrid of criteria corresponding to LSM and least modules method (LMM). When c=µ (intensity of outliers d =0) it becomes

whereas when d ® 100% (c® 0) it is equivalent to

M-estimates are not invariant with respect to scaling. To treat data with an arbitrary variance of density (not only unit density) a minimization problem

need to be solved in which unknown scale parameter s is estimated along with parameters q . The value of a is found as

where z is the number of sought-for parameters, j(x) is the standard normal density.

The variance of density of main observations s² is an estimate of residual variance provided that distribution of errors is normal (gross errors are absent):

The problem of parametric identification can be easily generalized for the multidimensional case when L properties of an equilibrium system are registered at each experimental point. In this case the criterion to be minimized is

Note that a high flexibility is inherent in the method of data analysis based on M-estimates due to additional variable, percentage of outliers d , in comparison with the method of least squares. A priori information about d is not available, and variation of this parameter may be considered as an adaptation of procedures of parametric identification to treated data.

To the top of the page

Algorithms

Stop criteria

Point estimates of unknown parameters |lg b_i>, <a_li> are found as values that give the minimum to the criterion function . Along with them the scale parameter s is estimated.

Newton and Gauss-Newton methods are realized in CLINP 2.1 for minimizing . The algorithm is as follows. After m iterations of s vector q is refined proceeding from s ^(m), and then new value of s is calculated. Iterations are repeated until the difference between two subsequent values of s will be not more than 1%.

Iterative scheme that makes the current values of sought-for parameters lgb_c (c means current) more precise at fixed value of s is represented by the formula

lg b⁺ = lg b_c + l Ч s,

in which vector

s = -H^-1 Ч g

specifies the direction of movement from lg b_c to lg b⁺, new values of parameters at next iteration, Н is Hess matrix, g is the gradient vector for , 0<lЈ1 is the step. Hess matrix is calculated exactly if Newton method is used and is approximated if Gauss-Newton method is used. M-estimates of unknown intensity factors a_li are calculated at each iteration of lg b using linear regression method.

Global convergence of the algorithm and fast local convergence are ensured by

1) correction of Hess matrix of indefinite signs based on its spectral decomposition;

2) bypass of saddle points in the direction of negative curvature;

3) linear strategy of searching l beyond vicinities of saddle points (Armijo-Goldstein rule [Armijo L. Pacific J. Math. 16 (1966) 1; Goldstein A. A. Constructive Real Analysis, Harper&Row, New York, 1967]).

Simplified calculations of Hess matrix in Newton-Gauss method may result in the loss of local convergence a) in the case of big discrepancies because of inappropriate choice of zero approximation for sought-for parameters or inappropriate choice of the set of chemical reactions included into the model; b) if minimized criterion is essentially nonlinear.

To the top of the page

Two main stop criteria are used in iterating lgb

The first one is constructed using the gradient of minimized function g(lg b). Vector g_r(lg b⁺) of the relative gradient is calculated during iterations, its components being

where is the value of minimized criterion. When norm ||g_r(lgb⁺)|| is small enough (smaller than 10^-7 –10^-4) iteration process is stopped.

Iteration process may also be stopped for another reason: when all increments lЧs_i are small. Fulfillment of the following condition

is the criterion in this case, where

typb is the "typical value" of sought-for parameters (recommended: 1-10),

e_step is the minimal admissible relative iteration step (recommended 10^-6 – 10^-4 ).

In addition, users may restrict the number of iterations by specifying the maximum number Itermax.

If Itermax=0, the program will calculate the equilibrium composition of a system, weighted discrepancies and characteristics of model adequacy for parameter values specified by a user.

To the top of the page

Residual variance and criterion c²

Skewness, kurtosis, mean residual and residual mean

If weighted discrepancies x_lk = w_lk^1/2ЧD_lk are independent random variables normally distributed with zero mean and unit density, then the residual variance is

where f = NЧL- z, z is the total number of calculated parameters. This random variable is distributed as c²/f with f degrees of freedom and unit mean. A model may be considered as adequate one if inequality

,

is fulfilled, where c_f²(a) is the 100a -percentage point for distribution c ² with f degrees of freedom. This condition is exact only for models linear in parameters. It is also used for checking the adequacy of complex formation models nonlinear in parameters, the significance level a being considered as approximate.

When Huber's M-estimates are used in parametric identification of models, resulted values of parameters a_li and lg b_i may not correspond to the minimum of residual variance. Accordingly, the conventional procedure of checking the model adequacy needs to be relaxed. This may be done by reducing f.

If a distribution of weighted discrepancies x_lk is differ from the normal one, its kurtosis being g₂, corrected value of f may be found using formula

f = (NЧ L - z) Ч {1 + 0.5Ч g_{2 Ч} (NЧ L - z)/NЧ L }^-1,

where int means the integer part of a number. Instead of g₂ its sampling estimate is used in CLINP 2.1.

To the top of the page

Proceeding from values of weighted discrepancies x_lk sampling estimates may be found for parameters of their distribution. CLINP 2.1 calculates:

where m₂ and m₃ are second (variance) and third central moments (respectively) of the distribution;

g₂ = m₄ / (m₂)² - 3,

where m₄ is the fourth central moment;

If x_lk are distributed normally with zero mean and unit density, expectations of the above quantities are as follows:

Sampling values may be compared with percentage points of distributions of corresponding statistics (Tables 1-3).

If for a given significance level a sampling estimates do not exceed their 100a -percentage points, the model may be considered as adequate one.

To the top of the page

Table 1. Percentage points of distribution

Table 2. Percentage points of distribution

Table 3. Percentage points of g₂ distribution

To the top of the page

After the iteration process is finished, the program calculates symmetrical covariance matrix D(lg b)of calculated parameters. Its diagonal elements are variance estimates for components of vector lg b – s²(lg b_j), other elements are their covariances cov(lg b_i, lg b_j):

This matrix is used for estimating confidence limits of parameters and their correlation coefficients.

Partial correlation coefficients r_ij characterize the mutual influence of parameters lg b_i and lg b_j provided that all other parameters are fixed:

General correlation coefficients s_ij serve as a measure of correlation of parameters lg b_i and lg b_j, when other parameters are considered as adjustable:

Multiple correlation coefficients R_i measures the degree of correlation of a given parameter with all other parameters lg b_i:

Correlation coefficient may vary in the interval [-1;1]. Zero values mean total independence of parameters, values ± 1 correspond to linear dependence of parameters (a model is overdetermined). The reconsideration of a chemical model (exclusion of one of complexes) is believed to be of need when the absolute value of corresponding correlation coefficient exceeds a threshold value chosen in the interval 0.95-0.98.

The program shows results of singular decomposition of Jacoby matrix

J = ¶ A / ¶ lg b = U S V^T,

where U and V  are orthogonal matrices, S is a rectangular matrix consisted of diagonal and zero blocks. Diagonal elements of S are singular numbers k_i.

Small singular numbers (k_i/k_max< 10^-4 – 10^-6) correspond to linear combinations of sought-for parameters lg b_j. The criterion function is insensitive to them.

Coefficients of these linear combinations are situated in columns of matrix V with numbers corresponding to that of small singular values. Not infrequently these coefficients provide sufficient information for detecting an odd complex. In more involved cases this may be easily done by analyzing equilibrium concentrations of reagents (an unresolvable linear combination of logarithms of constants appears when expressions of sought-for constants include concentrations of reagents not represented in the system).

To the top of the page

To the Help main page