def covariance(self, *, atol=1e-8) -> np.ndarray:
        """Compute the covariance matrix for the estimators of all settings.

        Like `variance`, this is the covariance of the sampling distribution
        of the sample mean. Practically, it is the 'normal' covariance
        divided by the number of observations (bitstrings).

        Args:
            atol: The absolute tolerance for asserting coefficients are real.

        Raises:
            ValueError: If there are no measurements.
        """
        if len(self.bitstrings) == 0:
            raise ValueError("No measurements")

        all_obs_vals = np.array(
            [
                _obs_vals_from_measurements(
                    bitstrings=self.bitstrings,
                    qubit_to_index=self._qubit_to_index,
                    observable=setting.observable,
                    atol=atol,
                )
                for setting in self._simul_settings
            ]
        )

        # https://github.com/numpy/numpy/issues/11502
        if all_obs_vals.shape[0] == 1:
            cov = np.array([[np.var(all_obs_vals[0], ddof=1)]])
            return cov

        cov = np.cov(all_obs_vals, ddof=1) / all_obs_vals.shape[1]
        return cov
Compute the covariance matrix for the estimators of all settings.

        Like `variance`, this is the covariance of the sampling distribution
        of the sample mean. Practically, it is the 'normal' covariance
        divided by the number of observations (bitstrings).

        Args:
            atol: The absolute tolerance for asserting coefficients are real.

        Raises:
            ValueError: If there are no measurements.
        
cirq-core/cirq/work/observable_measurement_data.py