Merge pull request #3 from kaklise/parmest-covariance

Minor changes to parmest for covariance estimation

doc/OnlineDocs/contributed_packages/parmest/api.rst

@@ -3,11 +3,22 @@
 API
 ============
 
+parmest
+---------
 .. automodule:: pyomo.contrib.parmest.parmest
     :members:
     :undoc-members:
     :show-inheritance:
 
+scenariocreator
+------------------
+.. automodule:: pyomo.contrib.parmest.scenariocreator
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+graphics
+---------
 .. automodule:: pyomo.contrib.parmest.graphics
     :members:
     :undoc-members:

doc/OnlineDocs/contributed_packages/parmest/covariance.rst

@@ -0,0 +1,9 @@
+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 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`.

doc/OnlineDocs/contributed_packages/parmest/datarec.rst

@@ -1,6 +1,6 @@
 .. _datarecsection:
 
-Data Reconciliation using parmest
+Data Reconciliation 
 =================================
 
 The method :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`

doc/OnlineDocs/contributed_packages/parmest/driver.rst

@@ -1,6 +1,6 @@
 .. _driversection:
 
-Parameter Estimation using parmest
+Parameter Estimation 
 ==================================
 
 Parameter Estimation using parmest requires a Pyomo model, experimental
@@ -145,12 +145,3 @@ 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 :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 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`.

doc/OnlineDocs/contributed_packages/parmest/index.rst

@@ -12,11 +12,12 @@ confidence regions and subsequent creation of scenarios for PySP.
    installation.rst
    driver.rst
    datarec.rst
+   covariance.rst
+   scencreate.rst
    graphics.rst
    examples.rst
    parallel.rst
    api.rst
-   scencreate.rst
 
 Indices and Tables
 ------------------

doc/OnlineDocs/contributed_packages/parmest/scencreate.rst

@@ -8,9 +8,6 @@ object, which has methods to add ``ParmestScen`` scenario objects to a
 ``ScenarioSet`` object, which can write them to a csv file or output them
 via an iterator method.
 
-Example
--------
-
 This example is in the semibatch subdirectory of the examples directory in
 the file ``scencreate.py``. It creates a csv file with scenarios that
 correspond one-to-one with the experiments used as input data. It also
@@ -23,12 +20,3 @@ scenarios to the screen, accessing them via the ``ScensItator`` a ``print``
 .. note::
    This example may produce an error message your version of Ipopt is not based
    on a good linear solver.
-
-
-API
----
-
-.. automodule:: pyomo.contrib.parmest.scenariocreator
-    :members:
-    :undoc-members:
-    :show-inheritance:

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

@@ -44,9 +44,8 @@ def SSE(model, data):
 parmest.pairwise_plot(bootstrap_theta, title='Bootstrap theta')
 parmest.pairwise_plot(bootstrap_theta, theta, 0.8, ['MVN', 'KDE', 'Rect'], 
                       title='Bootstrap theta with confidence regions')
-'''
+
 ### Likelihood ratio test
-'''
 asym = np.arange(10, 30, 2)
 rate = np.arange(0, 1.5, 0.1)
 theta_vals = pd.DataFrame(list(product(asym, rate)), columns=theta_names)

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

@@ -8,14 +8,10 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
-import numpy as np
 import pandas as pd
-from itertools import product
 import pyomo.contrib.parmest.parmest as parmest
 from rooney_biegler import rooney_biegler_model
 
-### Parameter estimation
-
 # Vars to estimate
 theta_names = ['asymptote', 'rate_constant']
 

pyomo/contrib/parmest/parmest.py

@@ -780,7 +780,7 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None, calc_cov
         assert isinstance(bootlist, (type(None), list))
 
         return self._Q_opt(solver=solver, return_values=return_values,
-                           bootlist=bootlist, calc_cov)
+                           bootlist=bootlist, calc_cov=calc_cov)
 
 
     def theta_est_bootstrap(self, bootstrap_samples, samplesize=None,