Parmest Quick Start Guide

This quick start guide shows how to use parmest to estimate model parameters from experimental data as well as compute their uncertainty. The model and data used in this guide were taken from [RB01].

The mathematical model of interest is:

\[y_i(\theta_1, \theta_2, t_i) = \theta_1 \left(1 - e^{-\theta_2 t_i} \right), \quad \forall \quad i \, \in \, {1, \ldots, n}\]

Where \(y\) is the observation of the measured variable, \(t\) is the time, \(\theta_1\) is the asymptote, and \(\theta_2\) is the rate constant.

The experimental data is given in the table below:

Table 6 Data
hour	y
1	8.3
2	10.3
3	19.0
4	16.0
5	15.6
7	19.8

To use parmest to estimate \(\theta_1\) and \(\theta_2\) from the data, we provide the following detailed steps:

Step 0: Import parmest and Pandas

Before solving the parameter estimation problem, the following code must be executed to import the required packages for parameter estimation in parmest:

>>> import pyomo.contrib.parmest.parmest as parmest
>>> import pandas as pd

Step 1: Create an Experiment Class

Parmest requires that the user create an Experiment class that builds an annotated Pyomo model denoting experiment outputs, unknown parameters, and measurement errors using Pyomo Suffix components.

m.experiment_outputs maps the experiment output (or measurement) terms in the model (Pyomo Param, Var, or Expression) to their associated data values (float, int).
m.unknown_parameters maps the model parameters to estimate (Pyomo Param or Var) to their component unique identifier (Pyomo ComponentUID) which is used to identify equivalent parameters across multiple experiments. Within parmest, any parameters that are to be estimated are converted to unfixed variables. Variables that are to be estimated are also unfixed.
m.measurement_error maps the measurement error (float, int) of the experiment output, or measurement (Pyomo Param, Var, or Expression) defined in the model.

The experiment class has one required method:

get_labeled_model which returns the labeled Pyomo model.

An example Experiment class is shown below.

Listing 2 RooneyBieglerExperiment class from the parmest example

class RooneyBieglerExperiment(Experiment):
    """Experiment class for Rooney-Biegler parameter estimation and design of experiments.

    This class wraps the Rooney-Biegler exponential model for use with Pyomo's
    parmest (parameter estimation) and DoE (Design of Experiments) tools.
    """

    def __init__(self, data, measure_error=None, theta=None):
        """Initialize a Rooney-Biegler experiment instance.

        Parameters
        ----------
        data : pandas.Series or dict
            Experiment data containing 'hour' (time input) and 'y' (response) values.
        measure_error : float, optional
            Standard deviation of measurement error for the response variable.
            Required for DoE and covariance estimation. Default is None.
        theta : dict, optional
            Initial parameter values with 'asymptote' and 'rate_constant' keys.
            Default is None, which uses {'asymptote': 15, 'rate_constant': 0.5}.
        """
        self.data = data
        self.model = None
        self.measure_error = measure_error
        self.theta = theta

    def create_model(self):
        # rooney_biegler_model expects a dataframe
        if hasattr(self.data, 'to_frame'):
            # self.data is a pandas Series
            data_df = self.data.to_frame().transpose()
        else:
            # self.data is a dict
            data_df = pd.DataFrame([self.data])
        self.model = rooney_biegler_model(data_df, theta=self.theta)

    def label_model(self):

        m = self.model

        # Add experiment outputs as a suffix
        # Experiment outputs suffix is required for parmest
        m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.experiment_outputs.update([(m.y, self.data['y'])])

        # Add unknown parameters as a suffix
        # Unknown parameters suffix is required for both Pyomo.DoE and parmest
        m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.unknown_parameters.update(
            (k, pyo.value(k)) for k in [m.asymptote, m.rate_constant]
        )

        # Add measurement error as a suffix
        # Measurement error suffix is required for Pyomo.DoE and
        #  `cov` estimation in parmest
        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.measurement_error.update([(m.y, self.measure_error)])

        # Add hour as an experiment input
        # Experiment inputs suffix is required for Pyomo.DoE
        m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.experiment_inputs.update([(m.hour, self.data['hour'])])

    def get_labeled_model(self):
        if self.model is None:
            self.create_model()
            self.label_model()
        return self.model

Step 2: Load the Data and Create a List of Experiments

Load the experimental data into Python and create an instance of your Experiment class for each set of experimental data. In this example, each measurement of y is treated as a separate experiment.

>>> 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"])
>>> exp_list = []
>>> for i in range(data.shape[0]):
...     exp_list.append(RooneyBieglerExperiment(data.loc[i, :]))

Step 3: Create the Estimator Object

To use parmest, the user creates an Estimator object which includes the following methods:

`theta_est`	Parameter estimation using all scenarios in the data
`cov_est`	Covariance matrix calculation using all scenarios in the data
`theta_est_bootstrap`	Parameter estimation using bootstrap resampling of the data
`theta_est_leaveNout`	Parameter estimation where N data points are left out of each sample
`objective_at_theta`	Objective value for each theta
`confidence_region_test`	Confidence region test to determine if theta values are within a rectangular, multivariate normal, or Gaussian kernel density distribution for a range of alpha values
`likelihood_ratio_test`	Likelihood ratio test to identify theta values within a confidence region using the \(\chi^2\) distribution
`leaveNout_bootstrap_test`	Leave-N-out bootstrap test to compare theta values where N data points are left out to a bootstrap analysis using the remaining data, results indicate if theta is within a confidence region determined by the bootstrap analysis

Additional functions are available in parmest to plot results and fit distributions to theta values.

`pairwise_plot`	Plot pairwise relationship for theta values, and optionally alpha-level confidence intervals and objective value contours
`grouped_boxplot`	Plot a grouped boxplot to compare two datasets
`grouped_violinplot`	Plot a grouped violinplot to compare two datasets
`fit_rect_dist`	Fit an alpha-level rectangular distribution to theta values
`fit_mvn_dist`	Fit a multivariate normal distribution to theta values
`fit_kde_dist`	Fit a Gaussian kernel-density distribution to theta values

A Estimator object can be created using the following code. A description of the arguments are listed below.

>>> pest = parmest.Estimator(exp_list, obj_function="SSE")

Alternatively, the weighted sum of squared errors objective can be used.

>>> pest = parmest.Estimator(exp_list, obj_function="SSE_weighted")

Optionally, solver options can be supplied, e.g.,

>>> solver_options = {"max_iter": 6000}
>>> pest = parmest.Estimator(exp_list, obj_function="SSE", solver_options=solver_options)

Objective function

The obj_function keyword argument is used to specify the objective function to use for parameter estimation if the user has not implemented their own custom objective function. 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. If the user wants to use an objective that is different from the built-in options, a custom objective function can be specified in the user’s model, however, covariance matrix estimation (see Covariance Matrix Estimation Section) is not supported for custom objective functions.

When declaring a custom objective function, parmest assumes the model has the structure of a two-stage stochastic programming problem so the objective function should be implemented using Pyomo Expressions for the first stage cost (named “FirstStageCost”) and the second stage cost (named “SecondStageCost”). For parameter estimation problems the first stage cost is usually set to zero and the second stage cost is usually defined as the deviation between the model and the observations.

Step 4: Estimate the Parameters

After creating the Estimator object with the desired objective function, solve the parameter estimation problem by calling theta_est, e.g.,

>>> pest = parmest.Estimator(exp_list, obj_function="SSE")
>>> obj_val, theta_val = pest.theta_est()

More Examples Beyond this Quick Guide

More detailed examples, such as parameter estimation of reaction kinetics are provided in the Examples Section.