Updated covariance.rst and Xinhong comment on parmest.py

Author: slilonfe5

doc/OnlineDocs/explanation/analysis/parmest/covariance.rst

@@ -1,16 +1,88 @@
 Covariance Matrix Estimation
 =================================
 
-If the optional argument ``calc_cov=True`` is specified for :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`, 
-parmest will calculate the covariance matrix :math:`V_{\theta}` as follows:
-
-.. math::
-   V_{\theta} = 2 \sigma^2 H^{-1} 
-
-This formula assumes all measurement errors are independent and identically distributed with 
-variance :math:`\sigma^2`. :math:`H^{-1}` is the inverse of the Hessian matrix for an unweighted 
-sum of least squares problem. Currently, the covariance approximation is only valid if the 
-objective given to parmest is the sum of squared error. Moreover, parmest approximates the 
-variance of the measurement errors as :math:`\sigma^2 = \frac{1}{n-l} \sum e_i^2` where :math:`n` is 
-the number of data points, :math:`l` is the number of fitted parameters, and :math:`e_i` is the 
-residual for experiment :math:`i`.
+The uncertainty in the estimated parameters is quantified using the covariance matrix.
+The diagonal of the covariance matrix contains the variance of the estimated parameters.
+Assuming Gaussian independent and identically distributed measurement errors, the
+covariance matrix of the estimated parameters can be computed using the following
+methods which have been implemented in parmest.
+
+1. Reduced Hessian Method
+
+    .. math::
+       V_{\boldsymbol{\theta}} = 2 \sigma^2 \left(\frac{\partial^2 \text{SSE}}
+        {\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}\right)^{-1}_{\boldsymbol{\theta}
+        = \boldsymbol{\theta}^*}
+
+    Where SSE is the sum of squared errors, WSSE is the weighted SSE,
+    :math:`\boldsymbol{\theta}` are the unknown parameters, :math:`\boldsymbol{\theta^*}`
+    are the estimate of the unknown parameters, and :math:`\sigma^2` is the variance of
+    the measurement error. When the standard deviation of the measurement error is not
+    supplied by the user, parmest approximates the variance of the measurement error as
+    :math:`\sigma^2 = \frac{1}{n-l} \sum e_i^2` where :math:`n` is the number of data
+    points, :math:`l` is the number of fitted parameters, and :math:`e_i` is the residual
+    for experiment :math:`i`.
+
+2. Finite Difference Method
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left( \sum_{r = 1}^n \mathbf{G}_{r}^{\mathrm{T}} \mathbf{W}
+        \mathbf{G}_{r} \right)^{-1}
+
+    This method uses central finite difference to compute the Jacobian matrix,
+    :math:`\mathbf{G}_{r}`, which is the sensitivity of the measured variables with
+    respect to the parameters, `\boldsymbol{\theta}`. :math:`\mathbf{W}` is a diagonal
+    matrix containing the inverse of the variance of the measurement errors,
+    :math:`\sigma^2`.
+
+3. Automatic Differentiation Method
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left( \sum_{r = 1}^n \mathbf{G}_{\text{kaug},\, r}^{\mathrm{T}}
+        \mathbf{W} \mathbf{G}_{\text{kaug},\, r} \right)^{-1}
+
+    This method uses the model optimality (KKT) condition to compute the Jacobian matrix,
+    :math:`\mathbf{G}_{\text{kaug},\, r}`.
+
+In parmest, the covariance matrix can be calculated after defining the
+:class:`~pyomo.contrib.parmest.parmest.Estimator` object and estimating the unknown
+parameters using :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`. To
+estimate the covariance matrix, call
+:class:`~pyomo.contrib.parmest.parmest.Estimator.cov_est` and pass it the number
+of data points, e.g.,
+
+.. testsetup:: *
+    :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
+
+    # Data
+    import pandas as pd
+    data = pd.DataFrame(
+        data=[[1, 8.3], [2, 10.3], [3, 19.0],
+              [4, 16.0], [5, 15.6], [7, 19.8]],
+        columns=['hour', 'y'],
+    )
+    num_data = len(data)
+
+    # Create an experiment list
+    from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import RooneyBieglerExperiment
+    exp_list = []
+    for i in range(data.shape[0]):
+        exp_list.append(RooneyBieglerExperiment(data.loc[i, :]))
+
+.. doctest::
+    :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
+
+    >>> import pyomo.contrib.parmest.parmest as parmest
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
+    >>> obj_val, theta_val = pest.theta_est()
+    >>> cov = pest.cov_est(cov_n=num_data)
+
+Optionally, one of the three methods; "reduced_hessian", "finite_difference",
+and "automatic_differenciation_kaug" can be supplied for the covariance calculation,
+e.g.,
+
+.. doctest::
+    :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
+
+    >>> cov_method = "reduced_hessian"
+    >>> cov = pest.cov_est(cov_n=num_data, method=cov_method)

pyomo/contrib/parmest/parmest.py

@@ -41,6 +41,7 @@
 import re
 import importlib as im
 import logging
+import warnings
 import types
 import json
 from collections.abc import Callable
@@ -259,9 +260,7 @@ def SSE_weighted(model):
     _check_model_labels_helper(model, logging_level=logging.ERROR)
 
     # Check that measurement errors exist
-    if hasattr(model, "measurement_error"):
-        pass
-    else:
+    if not hasattr(model, "measurement_error"):
         raise AttributeError(
             'Experiment model does not have suffix "measurement_error". '
             '"measurement_error" is a required suffix for the "SSE_weighted" '
@@ -270,16 +269,24 @@ def SSE_weighted(model):
 
     # check if all the values of the measurement error standard deviation
     # have been supplied
-    if all(
+    all_known_errors = all(
         model.measurement_error[y_hat] is not None for y_hat in model.experiment_outputs
-    ):
-        # calculate the weighted SSE between the prediction and observation of the
-        # measured variables
-        expr = (1 / 2) * sum(
-            ((y - y_hat) / model.measurement_error[y_hat]) ** 2
-            for y_hat, y in model.experiment_outputs.items()
-        )
-        return expr
+    )
+
+    if all_known_errors:
+        # calculate the weighted SSE between the prediction
+        # and observation of the measured variables
+        try:
+            expr = (1 / 2) * sum(
+                ((y - y_hat) / model.measurement_error[y_hat]) ** 2
+                for y_hat, y in model.experiment_outputs.items()
+            )
+            return expr
+        except ZeroDivisionError:
+            raise ValueError(
+                'Division by zero encountered in the "SSE_weighted" objective. '
+                'One or more values of the measurement error are zero.'
+            )
     else:
         raise ValueError(
             'One or more values are missing from "measurement_error". All values of '
@@ -582,36 +589,43 @@ def _finite_difference_FIM(
     # computing the condition number of the Jacobian matrix
     cond_number_jac = np.linalg.cond(J)
     if logging_level == logging.INFO:
-        logger.info(
-            f"The condition number of the Jacobian matrix " f"is {cond_number_jac}"
-        )
+        logger.info(f"The condition number of the Jacobian matrix is {cond_number_jac}")
 
     # grab the model
     model = _get_labeled_model_helper(experiment)
 
     # extract the measured variables and measurement errors
     y_hat_list = [y_hat for y_hat, y in model.experiment_outputs.items()]
 
-    # check if the model has a 'measurement_error' attribute and the measurement
-    # error standard deviation has been supplied
-    if hasattr(model, "measurement_error") and all(
+    # check if the model has a 'measurement_error' attribute and
+    # the measurement error standard deviation has been supplied
+    all_known_errors = all(
         model.measurement_error[y_hat] is not None for y_hat in model.experiment_outputs
-    ):
+    )
+
+    if hasattr(model, "measurement_error") and all_known_errors:
         error_list = [
             model.measurement_error[y_hat] for y_hat in model.experiment_outputs
         ]
 
-        # compute the matrix of the inverse of the measurement variance
-        # the following assumes independent measurement errors
-        W = np.diag([1 / (err**2) for err in error_list])
+        # compute the matrix of the inverse of the measurement error variance
+        # the following assumes independent and identically distributed
+        # measurement errors
+        try:
+            W = np.diag([1 / (err**2) for err in error_list])
+        except ZeroDivisionError:
+            raise ValueError(
+                'Division by zero encountered in computing the covariance matrix. '
+                'One or more values of the measurement error are zero.'
+            )
 
-        # check if error list is consistent
+        # check if the error list is consistent
         if len(error_list) == 0 or len(y_hat_list) == 0:
             raise ValueError(
                 "Experiment outputs and measurement errors cannot be empty."
             )
 
-        # check if the dimension of error_list is same with that of y_hat_list
+        # check if the dimension of error_list is the same with that of y_hat_list
         if len(error_list) != len(y_hat_list):
             raise ValueError(
                 "Experiment outputs and measurement errors are not the same length."
@@ -717,16 +731,24 @@ def _kaug_FIM(experiment, theta_vals, solver, tee, estimated_var=None):
     kaug_jac = np.array(jac).T
 
     # compute FIM
-    # compute matrix of the inverse of the measurement variance
-    # The following assumes independent measurement error.
+    # compute the matrix of the inverse of the measurement error variance
+    # the following assumes independent and identically distributed
+    # measurement errors
     W = np.zeros((len(model.measurement_error), len(model.measurement_error)))
     all_known_errors = all(
         model.measurement_error[y_hat] is not None for y_hat in model.experiment_outputs
     )
+
     count = 0
     for k, v in model.measurement_error.items():
         if all_known_errors:
-            W[count, count] = 1 / (v**2)
+            try:
+                W[count, count] = 1 / (v**2)
+            except ZeroDivisionError:
+                raise ValueError(
+                    'Division by zero encountered in computing the covariance matrix. '
+                    'One or more values of the measurement error are zero.'
+                )
         else:
             W[count, count] = 1 / estimated_var
         count += 1