Skip to content

Commit

Permalink
Moved 'calc_cov' to 'theta_est'. Added a paragraph in the documentation.
Browse files Browse the repository at this point in the history
  • Loading branch information
adowling2 committed Jun 1, 2020
1 parent 4a2549f commit 5f217d7
Show file tree
Hide file tree
Showing 5 changed files with 28 additions and 18 deletions.
10 changes: 10 additions & 0 deletions doc/OnlineDocs/contributed_packages/parmest/driver.rst
Original file line number Diff line number Diff line change
Expand Up @@ -144,3 +144,13 @@ include `model` and `data` and 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.

Covariance matrix estimation
----------------------------

If the optional argument ``calc_cov=True`` is specified for ``theta_est()``, parmest will calculate the covariance matrix of the fitted parameters 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`. `H^{-1}` is the inverse of the Hessian matrix for a 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`.
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ def SSE(model, data):

solver_options = {"max_iter": 6000} # not really needed in this case

pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE, solver_options, calc_cov=False)
pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE, solver_options)
obj, theta = pest.theta_est()
print(obj)
print(theta)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,8 @@ def SSE(model, data):

solver_options = {"max_iter": 6000} # not really needed in this case

pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE, solver_options, calc_cov=True)
obj, theta, cov = pest.theta_est()
pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE, solver_options)
obj, theta, cov = pest.theta_est(calc_cov=True)
print(obj)
print(theta)
print(cov)
26 changes: 13 additions & 13 deletions pyomo/contrib/parmest/parmest.py
Original file line number Diff line number Diff line change
Expand Up @@ -332,11 +332,9 @@ class Estimator(object):
If True, print diagnostics from the solver
solver_options: dict, optional
Provides options to the solver (also the name of an attribute)
calc_cov: bool, optional
Indicates if the covariance matrix should be computed and returned
"""
def __init__(self, model_function, data, theta_names, obj_function=None,
tee=False, diagnostic_mode=False, solver_options=None, calc_cov=False):
tee=False, diagnostic_mode=False, solver_options=None):

self.model_function = model_function
self.callback_data = data
Expand All @@ -353,9 +351,7 @@ def __init__(self, model_function, data, theta_names, obj_function=None,

self._second_stage_cost_exp = "SecondStageCost"
self._numbers_list = list(range(len(data)))

self.calc_cov = calc_cov



def _create_parmest_model(self, data):
"""
Expand Down Expand Up @@ -421,7 +417,7 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None):


def _Q_opt(self, ThetaVals=None, solver="ef_ipopt",
return_values=[], bootlist=None):
return_values=[], bootlist=None, calc_cov=False):
"""
Set up all thetas as first stage Vars, return resulting theta
values as well as the objective function value.
Expand Down Expand Up @@ -473,7 +469,7 @@ def _Q_opt(self, ThetaVals=None, solver="ef_ipopt",

assert not (need_gap and self.calc_cov), "Calculating both the gap and reduced hessian (covariance) is not currently supported."

if not self.calc_cov:
if not calc_cov:
# Do not calculate the reduced hessian

solver = SolverFactory('ipopt')
Expand Down Expand Up @@ -523,7 +519,7 @@ def _Q_opt(self, ThetaVals=None, solver="ef_ipopt",

objval = stsolver.root_E_obj()

if self.calc_cov:
if calc_cov:
# Calculate the covariance matrix

# Extract number of data points considered
Expand Down Expand Up @@ -553,12 +549,12 @@ def _Q_opt(self, ThetaVals=None, solver="ef_ipopt",
vals[var] = temp
var_values.append(vals)
var_values = pd.DataFrame(var_values)
if self.calc_cov:
if calc_cov:
return objval, thetavals, var_values, cov
else:
return objval, thetavals, var_values

if self.calc_cov:
if calc_cov:
return objval, thetavals, cov
else:
return objval, thetavals
Expand Down Expand Up @@ -750,7 +746,7 @@ def _get_sample_list(self, samplesize, num_samples, replacement=True):

return samplelist

def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None):
def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None, calc_cov=False):
"""
Parameter estimation using all scenarios in the data
Expand All @@ -762,6 +758,8 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None):
List of Variable names used to return values from the model
bootlist: list, optional
List of bootstrap sample numbers, used internally when calling theta_est_bootstrap
calc_cov: boolean, optional
If True, calculate and return the covariance matrix (only for "ef_ipopt" solver)
Returns
-------
Expand All @@ -774,13 +772,15 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None):
Hessian: dict
A dictionary of dictionaries for the Hessian.
The Hessian is not returned if the solver is ef_ipopt.
cov: numpy.array
Covariance matrix of the fitted parameters (only for ef_ipopt)
"""
assert isinstance(solver, str)
assert isinstance(return_values, list)
assert isinstance(bootlist, (type(None), list))

return self._Q_opt(solver=solver, return_values=return_values,
bootlist=bootlist)
bootlist=bootlist, calc_cov)


def theta_est_bootstrap(self, bootstrap_samples, samplesize=None,
Expand Down
4 changes: 2 additions & 2 deletions pyomo/contrib/parmest/tests/test_parmest.py
Original file line number Diff line number Diff line change
Expand Up @@ -226,10 +226,10 @@ def SSE(model, data):
expr = sum((data.y[i] - model.response_function[data.hour[i]])**2 for i in data.index)
return expr

self.pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE, calc_cov=True)
self.pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE)

def test_theta_est(self):
objval, thetavals, cov = self.pest.theta_est()
objval, thetavals, cov = self.pest.theta_est(calc_cov=True)

self.assertAlmostEqual(objval, 4.3317112, places=2)
self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) # 19.1426 from the paper
Expand Down

0 comments on commit 5f217d7

Please sign in to comment.