Extended Parmest Capability for weighted SSE objective

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

@@ -1,16 +1,121 @@
 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:
+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.
 
-.. math::
-   V_{\theta} = 2 \sigma^2 H^{-1} 
+1. Reduced Hessian Method
 
-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`.
+    When the objective function is the sum of squared errors (SSE):
+    :math:`\text{SSE} = \sum_{i = 1}^n \left(y_{i} - \hat{y}_{i}\right)^2`,
+    the covariance matrix is:
+
+    .. 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}^*}
+
+    When the objective function is the weighted SSE (WSSE):
+    :math:`\text{WSSE} = \frac{1}{2} \left(\mathbf{y} - f(\mathbf{x};\boldsymbol{\theta})\right)^\text{T}
+    \mathbf{W} \left(\mathbf{y} - f(\mathbf{x};\boldsymbol{\theta})\right)`,
+    the covariance matrix is:
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left(\frac{\partial^2 \text{WSSE}}
+        {\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}\right)^{-1}_{\boldsymbol{\theta}
+        = \boldsymbol{\theta}^*}
+
+    Where :math:`V_{\boldsymbol{\theta}}` is the covariance matrix of the estimated
+    parameters, :math:`y` are the observed measured variables, :math:`\hat{y}` are the
+    predicted measured variables, :math:`n` is the number of data points,
+    :math:`\boldsymbol{\theta}` are the unknown parameters, :math:`\boldsymbol{\theta^*}`
+    are the estimates of the unknown parameters, :math:`\mathbf{x}` are the decision
+    variables, and :math:`\mathbf{W}` is a diagonal matrix containing the inverse of the
+    variance of the measurement error, :math:`\sigma^2`. 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:`l` is the number of
+    fitted parameters, and :math:`e_i` is the residual for experiment :math:`i`.
+
+    In parmest, this method computes the inverse of the Hessian by scaling the
+    objective function (SSE or WSSE) with a constant probability factor.
+
+2. Finite Difference Method
+
+    In this method, the covariance matrix, :math:`V_{\boldsymbol{\theta}}`, is
+    calculated by applying the Gauss-Newton approximation to the Hessian,
+    :math:`\frac{\partial^2 \text{SSE}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}`
+    or
+    :math:`\frac{\partial^2 \text{WSSE}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}`,
+    leading to:
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left(\sum_{i = 1}^n \mathbf{G}_{i}^{\mathrm{T}} \mathbf{W}
+        \mathbf{G}_{i} \right)^{-1}
+
+    This method uses central finite difference to compute the Jacobian matrix,
+    :math:`\mathbf{G}_{i}`, for experiment :math:`i`, which is the sensitivity of
+    the measured variables with respect to the parameters, :math:`\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
+
+    Similar to the finite difference method, the covariance matrix is calculated as:
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left( \sum_{i = 1}^n \mathbf{G}_{\text{kaug},\, i}^{\mathrm{T}}
+        \mathbf{W} \mathbf{G}_{\text{kaug},\, i} \right)^{-1}
+
+    However, this method uses the model optimality (KKT) condition to compute the
+    Jacobian matrix, :math:`\mathbf{G}_{\text{kaug},\, i}`, for experiment :math:`i`.
+
+The covariance matrix calculation is only supported with the built-in objective
+functions "SSE" or "SSE_weighted".
+
+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, with the default method being "finite_difference", call
+the :class:`~pyomo.contrib.parmest.parmest.Estimator.cov_est` function, 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'],
+    )
+
+    # Create the Experiment class
+    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()
+
+Optionally, one of the three methods; "reduced_hessian", "finite_difference",
+and "automatic_differentiation_kaug" can be supplied for the covariance calculation,
+e.g.,
+
+.. doctest::
+    :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
+
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
+    >>> obj_val, theta_val = pest.theta_est()
+    >>> cov_method = "reduced_hessian"
+    >>> cov = pest.cov_est(method=cov_method)

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

@@ -44,8 +44,8 @@ The following example returns model values from a Pyomo Expression.
 
    >>> # Define objective
    >>> def SSE(model):
-   ...     expr = (model.experiment_outputs[model.y]
-   ...             - model.response_function[model.experiment_outputs[model.hour]]
+   ...     expr = (model.experiment_outputs[model.y[model.hour]]
+   ...             - model.y[model.hour]
    ...            ) ** 2
    ...     return expr
 

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

@@ -16,19 +16,21 @@ the model and the observations (typically defined as the sum of squared
 deviation between model values and observed values).
 
 If the Pyomo model is not formatted as a two-stage stochastic
-programming problem in this format, the user can supply a custom
-function to use as the second stage cost and the Pyomo model will be
+programming problem in this format, the user can choose either the
+built-in "SSE" or "SSE_weighted" objective functions, or supply a custom
+objective function to use as the second stage cost. The Pyomo model will then be
 modified within parmest to match the required specifications.
-The stochastic programming callback function is also defined within parmest.  The callback
-function returns a populated and initialized model for each scenario.
+The stochastic programming callback function is also defined within parmest.
+The callback function returns a populated and initialized model for each scenario.
 
-To use parmest, the user creates a :class:`~pyomo.contrib.parmest.parmest.Estimator` object 
-which includes the following methods:
+To use parmest, the user creates a :class:`~pyomo.contrib.parmest.parmest.Estimator`
+object which includes the following methods:
 
 .. autosummary::
    :nosignatures:
 
    ~pyomo.contrib.parmest.parmest.Estimator.theta_est
+   ~pyomo.contrib.parmest.parmest.Estimator.cov_est
    ~pyomo.contrib.parmest.parmest.Estimator.theta_est_bootstrap
    ~pyomo.contrib.parmest.parmest.Estimator.theta_est_leaveNout
    ~pyomo.contrib.parmest.parmest.Estimator.objective_at_theta
@@ -65,16 +67,9 @@ Section.
         columns=['hour', 'y'],
     )
 
-    # Sum of squared error function
-    def SSE(model):
-        expr = (
-            model.experiment_outputs[model.y]
-            - model.response_function[model.experiment_outputs[model.hour]]
-        ) ** 2
-        return expr
-
     # 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, :]))
@@ -83,15 +78,15 @@ Section.
     :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)
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
 
 Optionally, solver options can be supplied, e.g.,
 
 .. doctest::
     :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
 
     >>> solver_options = {"max_iter": 6000}
-    >>> pest = parmest.Estimator(exp_list, obj_function=SSE, solver_options=solver_options)
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE", solver_options=solver_options)
 
 
 List of experiment objects
@@ -137,17 +132,20 @@ expressions that are used to build an objective for the two-stage
 stochastic programming problem.  
 
 If the Pyomo model is not written as a two-stage stochastic programming problem in
-this format, and/or if the user wants to use an objective that is
-different than the original model, a custom objective function can be
-defined for parameter estimation.  The objective function has a single argument, 
-which is the model from a single experiment.
+this format, the user can select the "SSE" or "SSE_weighted" built-in objective
+functions. If the user wants to use an objective that is different from the built-in
+options, a custom objective function can be defined for parameter estimation. However,
+covariance matrix estimation will not support this custom objective function. The objective
+function (built-in or custom) has a single argument, which is the model from a single
+experiment.
 The objective function returns a Pyomo
 expression which is used to define "SecondStageCost".  The objective
 function can be used to customize data points and weights that are used
 in parameter estimation.
 
-Parmest includes one built in objective function to compute the sum of squared errors ("SSE") between the 
-``m.experiment_outputs`` model values and data values.
+Parmest includes two built-in objective functions ("SSE" and "SSE_weighted") to compute
+the sum of squared errors between the ``m.experiment_outputs`` model values and
+data values.
 
 Suggested initialization procedure for parameter estimation problems
 --------------------------------------------------------------------
@@ -162,4 +160,11 @@ estimation solve from the square problem solution, set optional argument ``solve
 argument ``(initialize_parmest_model=True)``. Different initial guess values for the fitted 
 parameters can be provided using optional argument `theta_values` (**Pandas Dataframe**)
 
-3. Solve parameter estimation problem by calling :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`
+3. Solve parameter estimation problem by calling
+:class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`, e.g.,
+
+.. doctest::
+    :skipif: not parmest_available
+
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
+    >>> obj_val, theta_val = pest.theta_est()

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

@@ -41,23 +41,32 @@ that for most experiments, only small parts of :math:`x` will change
 from one experiment to the next.
 
 The following least squares objective can be used to estimate parameter
-values, where data points are indexed by :math:`s=1,\ldots,S`
+values assuming Gaussian independent and identically distributed measurement
+errors, where data points are indexed by :math:`s=1,\ldots,S`
 
 .. math::
 
    \min_{{\theta}} Q({\theta};{\tilde{x}}, {\tilde{y}}) \equiv \sum_{s=1}^{S}q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) \;\;
 
-where
+where :math:`q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s})` can be:
 
-.. math::
+1. Sum of squared errors
+
+    .. math::
+
+       q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) =
+        \sum_{i=1}^{m}\left({\tilde{y}}_{s,i} - g_{i}({\tilde{x}}_{s};{\theta})\right)^{2}
+
+2. Weighted sum of squared errors
+
+    .. math::
 
-   q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) = \sum_{i=1}^{m}w_{i}\left[{\tilde{y}}_{si} - g_{i}({\tilde{x}}_{s};{\theta})\right]^{2}, 
+       q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) =
+        \sum_{i=1}^{m}\left(\frac{{\tilde{y}}_{s,i} - g_{i}({\tilde{x}}_{s};{\theta})}{w_i}\right)^{2}
 
 i.e., the contribution of sample :math:`s` to :math:`Q`, where :math:`w
-\in \Re^{m}` is a vector of weights for the responses. For
-multi-dimensional :math:`y`, this is the squared weighted :math:`L_{2}`
-norm and for univariate :math:`y` the weighted squared deviation.
-Custom objectives can also be defined for parameter estimation.
+\in \Re^{m}` is a vector containing the standard deviation of the measurement
+errors of :math:`y`. Custom objectives can also be defined for parameter estimation.
 
 In the applications of interest to us, the function :math:`g(\cdot)` is
 usually defined as an optimization problem with a large number of

pyomo/contrib/parmest/examples/rooney_biegler/bootstrap_example.py

@@ -27,8 +27,7 @@ def main():
     # Sum of squared error function
     def SSE(model):
         expr = (
-            model.experiment_outputs[model.y]
-            - model.response_function[model.experiment_outputs[model.hour]]
+            model.experiment_outputs[model.y[model.hour]] - model.y[model.hour]
         ) ** 2
         return expr
 

pyomo/contrib/parmest/examples/rooney_biegler/likelihood_ratio_example.py

@@ -28,8 +28,7 @@ def main():
     # Sum of squared error function
     def SSE(model):
         expr = (
-            model.experiment_outputs[model.y]
-            - model.response_function[model.experiment_outputs[model.hour]]
+            model.experiment_outputs[model.y[model.hour]] - model.y[model.hour]
         ) ** 2
         return expr
 

pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py

@@ -27,8 +27,7 @@ def main():
     # Sum of squared error function
     def SSE(model):
         expr = (
-            model.experiment_outputs[model.y]
-            - model.response_function[model.experiment_outputs[model.hour]]
+            model.experiment_outputs[model.y[model.hour]] - model.y[model.hour]
         ) ** 2
         return expr
 

pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py

@@ -26,19 +26,15 @@ def rooney_biegler_model(data):
     model.asymptote = pyo.Var(initialize=15)
     model.rate_constant = pyo.Var(initialize=0.5)
 
-    model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True)
-    model.y = pyo.Param(within=pyo.PositiveReals, mutable=True)
+    model.y = pyo.Var(data.hour, within=pyo.PositiveReals, initialize=5)
 
     def response_rule(m, h):
-        expr = m.asymptote * (1 - pyo.exp(-m.rate_constant * h))
-        return expr
+        return m.y[h] == m.asymptote * (1 - pyo.exp(-m.rate_constant * h))
 
-    model.response_function = pyo.Expression(data.hour, rule=response_rule)
+    model.response_function = pyo.Constraint(data.hour, rule=response_rule)
 
     def SSE_rule(m):
-        return sum(
-            (data.y[i] - m.response_function[data.hour[i]]) ** 2 for i in data.index
-        )
+        return sum((data.y[i] - m.y[data.hour[i]]) ** 2 for i in data.index)
 
     model.SSE = pyo.Objective(rule=SSE_rule, sense=pyo.minimize)
 
@@ -47,9 +43,10 @@ def SSE_rule(m):
 
 class RooneyBieglerExperiment(Experiment):
 
-    def __init__(self, data):
+    def __init__(self, data, measure_error=None):
         self.data = data
         self.model = None
+        self.measure_error = measure_error
 
     def create_model(self):
         # rooney_biegler_model expects a dataframe
@@ -61,22 +58,22 @@ def label_model(self):
         m = self.model
 
         m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.experiment_outputs.update(
-            [(m.hour, self.data['hour']), (m.y, self.data['y'])]
-        )
+        m.experiment_outputs.update([(m.y[self.data['hour']], self.data['y'])])
 
         m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.unknown_parameters.update(
             (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant]
         )
 
+        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.measurement_error.update([(m.y[self.data['hour']], self.measure_error)])
+
     def finalize_model(self):
 
         m = self.model
 
-        # Experiment output values
+        # Experiment input values
         m.hour = self.data['hour']
-        m.y = self.data['y']
 
     def get_labeled_model(self):
         self.create_model()