diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index 6643485fc9b..7f05b19b9ba 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -1074,16 +1074,13 @@ def __init__(
         assert isinstance(experiment_list, list)
         self.exp_list = experiment_list
 
-        # get the number of experiments
-        self.number_exp = _count_total_experiments(self.exp_list)
-
         # check if the experiment has a ``get_labeled_model`` function
         model = _get_labeled_model(self.exp_list[0])
 
         # check if the model has all the required suffixes
         _check_model_labels(model)
 
-        # populate keyword argument options
+        # check if the objective function is supplied correctly
         if isinstance(obj_function, str):
             try:
                 self.obj_function = ObjectiveType(obj_function)
@@ -1113,6 +1110,18 @@ def __init__(
                 f"{[e.value for e in RegularizationType]}."
             )
 
+        # get the number of data points
+        # this is used to compute the covariance matrix for the case
+        # when the measurement errors are not supplied by the user
+        get_measurement_error = getattr(model, "measurement_error", None)
+
+        all_unknown_errors = get_measurement_error is None or all(
+            model.measurement_error[y_hat] is None for y_hat in model.experiment_outputs
+        )
+
+        if self.obj_function == ObjectiveType.SSE and all_unknown_errors:
+            self.number_exp = _count_total_experiments(self.exp_list)
+
         self.tee = tee
         self.diagnostic_mode = diagnostic_mode
         self.solver_options = solver_options
@@ -1691,7 +1700,11 @@ def _cov_at_theta(self, method, solver, step):
 
                 # fix the value of the unknown parameters to the estimated values
                 for param in model.unknown_parameters:
-                    param.fix(self.estimated_theta[param.name])
+                    if param.is_indexed():
+                        for idx in param:
+                            param[idx].fix(self.estimated_theta[param[idx].name])
+                    else:
+                        param.fix(self.estimated_theta[param.name])
 
                 # re-solve the model with the estimated parameters
                 results = pyo.SolverFactory(solver).solve(model, tee=self.tee)
@@ -1720,12 +1733,6 @@ def _cov_at_theta(self, method, solver, step):
                 f"The sum of squared errors at the estimated parameter(s) is: {sse}"
             )
 
-        # Number of data points considered
-        n = self.number_exp
-
-        # Extract the number of fitted parameters
-        l = len(self.estimated_theta)
-
         """Calculate covariance assuming experimental observation errors are
         independent and follow a Gaussian distribution with constant variance.
 
@@ -1763,6 +1770,17 @@ def _cov_at_theta(self, method, solver, step):
 
                 # check if the user supplied the values of the measurement errors
                 if all(item is None for item in meas_error):
+                    # Number of data points considered
+                    n = self.number_exp
+
+                    # Extract the number of fitted parameters
+                    l = len(self.estimated_theta)
+
+                    assert n > l, (
+                        "The number of datapoints must be greater than the "
+                        "number of parameters to estimate."
+                    )
+
                     if cov_method == CovarianceMethod.reduced_hessian:
                         # in the "reduced_hessian" method, use the objective value
                         # to calculate the measurement error variance because this
@@ -2027,18 +2045,6 @@ def cov_est(self, method="finite_difference", solver="ipopt", step=1e-3):
         if not isinstance(step, float):
             raise TypeError("Expected a float for the step, e.g., 1e-2")
 
-        # number of unknown parameters
-        num_unknowns = max(
-            [
-                len(_expanded_unknown_parameter_info(experiment.get_labeled_model())[0])
-                for experiment in self.exp_list
-            ]
-        )
-        assert self.number_exp > num_unknowns, (
-            "The number of datapoints must be greater than the "
-            "number of parameters to estimate."
-        )
-
         return self._cov_at_theta(method=method, solver=solver, step=step)
 
     def theta_est_bootstrap(
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index f44331a86ac..af8fc9ca757 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -2272,9 +2272,12 @@ def test_indexed_unknown_parameter_names_are_expanded_consistently(self):
 
         self.assertEqual(pest._return_theta_names(), ["theta[a]", "theta[b]"])
 
+    @unittest.skipUnless(ipopt_available, "Test requires ipopt")
     def test_cov_est_counts_expanded_indexed_unknown_parameters(self):
         pest = _build_indexed_theta_estimator([(1.0, 2.0), (2.0, 4.0)])
 
+        obj, theta = pest.theta_est()
+
         with self.assertRaisesRegex(
             AssertionError,
             "The number of datapoints must be greater than the number of parameters to estimate.",