1+ .. _covariancesection :
2+
13Covariance Matrix Estimation
24============================
35
4- The uncertainty in the estimated parameters is quantified using the covariance matrix.
5- The diagonal of the covariance matrix contains the variance of the estimated parameters.
6- Assuming Gaussian independent and identically distributed measurement errors, the
7- covariance matrix of the estimated parameters can be computed using the following
8- methods which have been implemented in parmest.
6+ The goal is parameter estimation (see :ref: `driversection ` Section) is to estimate unknown model parameters
7+ from experimental data. When the model parameters are estimated from the data, their accuracy is measured by
8+ computing the covariance matrix. The diagonal of this covariance matrix contains the variance of the
9+ estimated parameters which is used to calculate their uncertainty. Assuming Gaussian independent and identically
10+ distributed measurement errors, the covariance matrix of the estimated parameters can be computed using the
11+ following methods which have been implemented in parmest.
912
10131. Reduced Hessian Method
1114
1215 When the objective function is the sum of squared errors (SSE) defined as
13- :math: `\text {SSE} = \sum _{i = 1 }^n
14- \left (\boldsymbol {y}_{i} - \hat { \boldsymbol {y}} _{i}\right )^2 `,
16+ :math: `\text {SSE} = \sum _{i = 1 }^{n}
17+ \left (\boldsymbol {y}_{i} - \boldsymbol {f}( \boldsymbol {x} _{i}; \boldsymbol { \theta }) \right )^2 `,
1518 the covariance matrix is:
1619
1720 .. math ::
1821 \boldsymbol {V}_{\boldsymbol {\theta }} = 2 \sigma ^2 \left (\frac {\partial ^2 \text {SSE}}
19- {\partial \boldsymbol {\theta } \partial \boldsymbol { \theta } }\right )^{-1 }_{\boldsymbol {\theta }
22+ {\partial \boldsymbol {\theta }^ 2 }\right )^{-1 }_{\boldsymbol {\theta }
2023 = \hat {\boldsymbol {\theta }}}
2124
2225 Similarly, when the objective function is the weighted SSE (WSSE) defined as
23- :math: `\text {WSSE} = \frac {1 }{2 } \sum _{i=1 }^{n}
24- \left (\boldsymbol {y}_{i} - \hat {\boldsymbol {y}}_{i}\right )^\text {T}
25- \boldsymbol {\Sigma }_{\boldsymbol {y}}^{-1 }
26- \left (\boldsymbol {y}_{i} - \hat {\boldsymbol {y}}_{i}\right )`,
26+ :math: `\text {WSSE} = \frac {1 }{2 } \sum _{i = 1 }^{n} \left (\boldsymbol {y}_{i} -
27+ \boldsymbol {f}(\boldsymbol {x}_{i};\boldsymbol {\theta })\right )^\text {T} \boldsymbol {\Sigma }_{\boldsymbol {y}}^{-1 }
28+ \left (\boldsymbol {y}_{i} - \boldsymbol {f}(\boldsymbol {x}_{i};\boldsymbol {\theta })\right )`,
2729 the covariance matrix is:
2830
2931 .. math ::
3032 \boldsymbol {V}_{\boldsymbol {\theta }} = \left (\frac {\partial ^2 \text {WSSE}}
31- {\partial \boldsymbol {\theta } \partial \boldsymbol { \theta } }\right )^{-1 }_{\boldsymbol {\theta }
33+ {\partial \boldsymbol {\theta }^ 2 }\right )^{-1 }_{\boldsymbol {\theta }
3234 = \hat {\boldsymbol {\theta }}}
3335
3436 Where :math: `\boldsymbol {V}_{\boldsymbol {\theta }}` is the covariance matrix of the estimated
3537 parameters :math: `\hat {\boldsymbol {\theta }}`, :math: `\boldsymbol {y}` are observations of the measured variables,
36- :math: `\hat {\boldsymbol {y}}` are model predictions of the measured variables,
37- :math: `n` is the number of experiments, and :math: `\boldsymbol {\Sigma }_{\boldsymbol {y}}` is the measurement error
38- covariance matrix, whose leading diagonal contains the inverse of the variance of the measurement errors,
39- :math: `\sigma ^2 `. When the standard deviation of the measurement error is not supplied by the user, parmest
40- approximates :math: `\sigma ^2 ` as: :math: `\hat {\sigma }^2 = \frac {1 }{n-l} \sum _{i=1 }^{n} e_i^2 `, where :math: `l` is
41- the number of fitted parameters, and :math: `e_i` is the residual between the data and model prediction for
42- experiment :math: `i`.
38+ :math: `n` is the number of experiments, :math: `\boldsymbol {\Sigma }_{\boldsymbol {y}}` is the measurement error
39+ covariance matrix, and :math: `\sigma ^2 ` is the variance of the measurement error. When the standard deviation of
40+ the measurement error is not supplied by the user, parmest approximates :math: `\sigma ^2 ` as:
41+ :math: `\hat {\sigma }^2 = \frac {1 }{n-l} \sum _{i=1 }^{n} e_i^2 `, where :math: `l` is the number of fitted parameters,
42+ and :math: `e_i` is the residual between the data and model for experiment :math: `i`.
4343
4444 In parmest, this method computes the inverse of the Hessian by scaling the
4545 objective function (SSE or WSSE) with a constant probability factor, :math: `\frac {1 }{n}`.
4646
47472. Finite Difference Method
4848
4949 In this method, the covariance matrix, :math: `\boldsymbol {V}_{\boldsymbol {\theta }}`, is
50- calculated by applying the Gauss-Newton approximation to the Hessian,
50+ computed by differentiating the Hessian,
5151 :math: `\frac {\partial ^2 \text {SSE}}{\partial \boldsymbol {\theta } \partial \boldsymbol {\theta }}`
5252 or
53- :math: `\frac {\partial ^2 \text {WSSE}}{\partial \boldsymbol {\theta } \partial \boldsymbol {\theta }}`,
54- leading to :
53+ :math: `\frac {\partial ^2 \text {WSSE}}{\partial \boldsymbol {\theta } \partial \boldsymbol {\theta }}`, and
54+ applying Gauss-Newton approximation which results in :
5555
5656 .. math ::
5757 \boldsymbol {V}_{\boldsymbol {\theta }} = \left (\sum _{i = 1 }^n \boldsymbol {G}_{i}^{\text {T}}
5858 \boldsymbol {\Sigma }_{\boldsymbol {y}}^{-1 } \boldsymbol {G}_{i} \right )^{-1 }
5959
60- This method uses central finite difference to compute the Jacobian matrix,
61- :math: `\boldsymbol {G}_{i}`, for experiment :math: `i`, which is the sensitivity of
62- the measured variables :math: `\boldsymbol {y}` with respect to the parameters, :math: `\boldsymbol {\theta }`.
63- :math: `\boldsymbol {\Sigma }_{\boldsymbol {y}}` is the measurement error covariance matrix, whose leading diagonal
64- contains the inverse of the variance of the measurement errors, :math: `\sigma ^2 `.
60+ where
61+
62+ .. math ::
63+ \boldsymbol {G}_{i} = \frac {\partial \boldsymbol {f}(\boldsymbol {x}_i;\boldsymbol {\theta })}
64+ {\partial \boldsymbol {\theta }}
65+
66+ This method uses central finite difference to compute the Jacobian matrix, :math: `\boldsymbol {G}_{i}`,
67+ for experiment :math: `i`.
68+
69+ .. math ::
70+ \boldsymbol {G}_{i}[:,\, k] \approx \frac {\boldsymbol {f}(\boldsymbol {x}_i;\theta _k + \Delta \theta _k)
71+ \vert _{\hat {\boldsymbol {\theta }}} - \boldsymbol {f}(\boldsymbol {x}_i;\theta _k - \Delta \theta _k)
72+ \vert _{\hat {\boldsymbol {\theta }}}}{2 \Delta \theta _k} \quad \forall \quad \theta _k \, \in \,
73+ [\theta _1 ,\cdots , \theta _p]
6574
6675 3. Automatic Differentiation Method
6776
@@ -71,23 +80,24 @@ methods which have been implemented in parmest.
7180 \boldsymbol {V}_{\boldsymbol {\theta }} = \left ( \sum _{i = 1 }^n \boldsymbol {G}_{\text {kaug},\, i}^{\text {T}}
7281 \boldsymbol {\Sigma }_{\boldsymbol {y}}^{-1 } \boldsymbol {G}_{\text {kaug},\, i} \right )^{-1 }
7382
74- However, this method uses the model optimality (KKT) condition to compute the Jacobian matrix,
75- :math: `\boldsymbol {G}_{\text {kaug},\, i}`, for experiment :math: `i`.
83+ However, this method uses implicit differentiation and the model-optimality or Karush–Kuhn–Tucker (KKT) conditions
84+ to compute the Jacobian matrix, :math: `\boldsymbol {G}_{\text {kaug},\, i}`, for experiment :math: `i`.
85+
86+ .. math ::
87+ \boldsymbol {G}_{\text {kaug},\, i} = \frac {\partial \boldsymbol {f}(\boldsymbol {x}_i,\boldsymbol {\theta })}
88+ {\partial \boldsymbol {\theta }} + \frac {\partial \boldsymbol {f}(\boldsymbol {x}_i,\boldsymbol {\theta })}
89+ {\partial \boldsymbol {x}_i}\frac {\partial \boldsymbol {x}_i}{\partial \boldsymbol {\theta }}
7690
77- The covariance matrix calculation is only supported with the built-in objective
78- functions "SSE" or "SSE_weighted".
91+ The covariance matrix calculation is only supported with the built-in objective functions "SSE" or "SSE_weighted".
7992
8093In parmest, the covariance matrix can be computed after creating the
8194:class: `~pyomo.contrib.parmest.experiment.Experiment ` class,
82- defining the :class: `~pyomo.contrib.parmest.parmest.Estimator ` object, and estimating the unknown parameters using
83- :class: `~pyomo.contrib.parmest.parmest.Estimator.theta_est ` (these steps were addressed in the
84- :ref: `driversection ` Section).
95+ defining the :class: `~pyomo.contrib.parmest.parmest.Estimator ` object,
96+ and estimating the model parameters using :class: `~pyomo.contrib.parmest.parmest.Estimator.theta_est `
97+ (all these steps were addressed in the :ref: `driversection ` Section).
8598
8699To estimate the covariance matrix, with the default method being "finite_difference", call
87- the :class: `~pyomo.contrib.parmest.parmest.Estimator.cov_est ` function as follows (note that this call should be made
88- after creating the :class: `~pyomo.contrib.parmest.experiment.Experiment ` class and the list of
89- :class: `~pyomo.contrib.parmest.experiment.Experiment ` objects as shown in the :ref: `driversection ` Section):
90- :
100+ the :class: `~pyomo.contrib.parmest.parmest.Estimator.cov_est ` function as follows:
91101
92102.. testsetup :: *
93103 :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
0 commit comments