Fixing and documenting how the Estimator calculates stds

qiskit/primitives/backend_estimator_v2.py

@@ -72,7 +72,7 @@ class _PreprocessedData:
 
 
 class BackendEstimatorV2(BaseEstimatorV2):
-    """Evaluates expectation values for provided quantum circuit and observable combinations
+    r"""Evaluates expectation values for provided quantum circuit and observable combinations
 
     The :class:`~.BackendEstimatorV2` class is a generic implementation of the
     :class:`~.BaseEstimatorV2` interface that is used to wrap a :class:`~.BackendV2`
@@ -87,7 +87,18 @@ class BackendEstimatorV2(BaseEstimatorV2):
     precludes doing any provider- or backend-specific optimizations.
 
     This class does not perform any measurement or gate mitigation, and, presently, is only
-    compatible with Pauli-based observables.
+    compatible with Pauli-based observables. More formally, given an observable of the type
+    :math:`O=\sum_{i=1}^Na_iP_i`, where :math:`a_i` is a complex number and :math:`P_i` is a
+    Pauli operator, the estimator calculates the expectation :math:`E(P_i)` of each :math:`P_i`
+    and finally calculates the expectation value of :math:`O` as
+    :math:`E(O)==\sum_{i=1}^Na_iE(P_i)`. The reported ``std`` is calculated as
+
+    .. math::
+
+        \frac{\sum_{i=1}^{n}a_i\sqrt{\textrm{Var}\big(P_i\big)}}{\sqrt{N}}\:,
+
+    where :math:`\textrm{Var}(P_i)` is the variance of :math:`P_i`, :math:`N=O(\epsilon^{-2})` is
+    the number of shots, and :math:`\epsilon` is the target precision.
 
     Each tuple of ``(circuit, observables, <optional> parameter values, <optional> precision)``,
     called an estimator primitive unified bloc (PUB), produces its own array-based result. The
@@ -254,8 +265,8 @@ def _postprocess_pub(
             for pauli, coeff in bc_obs[index].items():
                 expval, variance = expval_map[param_index, pauli]
                 evs[index] += expval * coeff
-                variances[index] += variance * coeff**2
-        stds = np.sqrt(variances / shots)
+                variances[index] += np.abs(coeff) * variance**0.5
+        stds = variances / np.sqrt(shots)
         data_bin = DataBin(evs=evs, stds=stds, shape=evs.shape)
         return PubResult(data_bin, metadata={"target_precision": pub.precision})