Confidence Intervals

The estimation of confidence intervals from least squares fits of nonlinear dynamic models against a collection of data sets over time is not standardised. The method used here follows closely the arguments given by Rawlings and Ekerdt [81] and Bard [5]. It is illustrated and extended here for the sake of clarity and the application to the experimental data obtained in this work.

The region of confidence around a fitted pair of parameters (e. g. X₁^∗ and X₂^∗) can be vi-sualised as the area inside an ellipse centred at X^∗ as sketched in Fig. 3.8. The ellipse is described by a covariance matrix derived from the data and the fit. Usually individual con-fidence intervals are given to each parameter of a parameter set. This corresponds to a box, which circumscribes the ellipse and thus overestimates the confidence region. An ellipse can be described by Eq. 3.11 and the corresponding confidence box by Eq. 3.12.

X^TAX =b (3.11)



 X₁^∗±

q b A⁻₁₁¹ X₂^∗±

q b A⁻₂₂¹



 (3.12)

X^∗ X1

X₂

q b A⁻₁₁¹ q

b A⁻₂₂¹

Figure 3.8.: The geometry of an ellipseX^TA X =bas confidence region around a parameter setX^∗.

When parameters are computed via a least squares optimisation of a nonlinear model, the description of all parameter sets inside the confidence interval can be formulated as Eq. 3.13 [81]. Its form is analogous to the ellipse described in Eq. 3.11.

(X−X^∗)^T ·H_X=X∗·(X−X^∗)≤2·s²·n_p·F(np,n_d−n_p,α) (3.13) with the (true/fitted) parameter set (X^∗)X, the variances², the Fisher probability functionF, the number of model parametersn_p, the number of data pointsn_d, the confidence levelα, and the HessianH of the objective function atX^∗.

H can be approximated asH ≈2J^TJ with the JacobianJ, as is done in the Gauß-Newton optimisation algorithm for solving least squares problems. In this application, the Jacobian can also be interpreted as the sensitivity of the model fit against the variation of parameters.

The confidence interval, or box for each parameter of the set can therefore be computed with Eq. 3.14, derived from Eqs. 3.11 through 3.13.

box_i= q

(2·s²·n_p·F(np,n_d−n_p,α))·(2J^TJ)⁻¹ (3.14)

On Estimating The Variance

If the variance does not change from data point to data point, the variance can be estimated from

s²=RSS(X^∗)

n_d−n_p (3.15)

3.3. Model Discrimination

0 100 200 300 400

0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1

t / [s]

normalised mass / [−]

synthetic data, s² of all datapoints actual data

0 100 200 300 400

0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1

t / [s]

normalised mass / [−]

synthetic data, s² of 1/4 of all datapoints actual data

Figure 3.9.: Visualisation of the estimated data variance by generating synthetic data from the model and noise withs². Left: Using all data points fors² estimation. Right:

Using only the “worst quarter” by deviation, thus better covering the data’s spread.

The actual data plotted are of Fe₃O₄ reductions at y_CO2 =0.52,Y_C =0.17,T = 750^◦C.

withRSS(X^∗)as the residual sum of squares at optimum.

Problems arise with a variance changing over time or data points, as is the case e. g. with the TGA data: All data points at t=0s have the same value of 1 as all data sets have been normalised w. r. t. to their initial values. Additionally, they have very similar values at the end of the data sets at steady state, since the reactions are limited by thermodynamic constraints.

It is between those two boundaries where the variance is almost zero, that the data sets differ significantly because of differing reaction rates.

Another problem is the prerequisite for Eq. (3.15) that alln_d data points are independent from each other. For the data points taken in a time series of a TGA experiment this is not the case. They do depend on the previous states of the system, which are represented by the previous measurements.

Therefore, if the estimation in Eq. 3.15 is applied to the data, the variance is underestimated as is visualised in Fig. 3.9. To estimate a more representative variance, only the quarter of all data points was used, with the highest deviation between fitted model and measurement (Eq. 3.15). This typically selected the second quarter of measurements in the time series.

s²= RSS_top25%(X^∗)

1/⁴n_d−n_p (3.16)

Figure 3.10.: The confidence intervals calculated from Eq. 3.14 might not be correct, because the applied kinetic model is strongly non-linear. The computed confidence inter-vals were therefore tested via the Monte Carlo method.

Testing The Estimated Confidence Intervals

During the process of estimating the confidence intervals, a number of approximations were applied. A statistical method for calculating confidence intervals for parameters of linear mod-els fitted to independent data points was expanded for nonlinear modmod-els and non-independent data points. The models were effectively linearised atX^∗for the determination of the Jacobian.

Furthermore, the Gauß-Newton approximation of the Hessian was used.

To test the computed confidence intervals, Rawlings [81] proposes the application of a Monte Carlo study for each parameter set. This method generates a large number of synthetic data sets from fitted kinetic parameters and tests, whether new kinetic parameters obtained from these synthetic data sets are within the original confidence intervals. The procedure includes the following steps:

1. Generate a synthetic data set from the model with original, fitted parameters and Gauss-ian noise with varGauss-iances².

2. Fit a new parameter set to that data set by least squares optimisation.

3. Repeat steps 1 & 2 a large number of times, e. g. 500.

4. Test, whether more than 95 % of the new fitted parameter sets from step 2 are within the α =0.95 confidence interval of the original fit (see Eq. 3.13). In a linear model with independent measurements this is always the case.

This Monte Carlo study is also illustrated in in Fig. 3.10 and 3.11.

3.3. Model Discrimination

Figure 3.11.: Example for Monte Carlo study with 20 synthetic data sets generated from pa-rameter set X^∗ with 95 % confidence interval. Parameter sets generated from synthetic data sets are marked by dots. In this illustration the study indicates that the stated confidence interval is not underestimated because no less than 95 % of generated parameter sets are in the original confidence ellipse.

All given parameter sets passed the Monte Carlos study. Approximately 95 % of the cor-responding synthetic parameter sets were within the 95 % confidence ellipse and much more inside the (larger) rectangular confidence box. An exemplary visualisation is done for the Avrami model fit of the reduction of Fe₃O₄ to FeO in Fig. 3.12. The estimated confidence intervals are therefore not impaired by the simplifications used in calculating them.