Return sampling error coveriance matrix Sigma eta

        Methodology source: Kjeldsen & Jones 2009, eq 9

        :param donor_catchments: Catchments to use as donors
        :type donor_catchments: list of :class:`Catchment`
        :return: 2-Dimensional, symmetric covariance matrix
        :rtype: :class:`numpy.ndarray`

def _matrix_sigma_eps(self, donor_catchments):
        """
        Return sampling error coveriance matrix Sigma eta

        Methodology source: Kjeldsen & Jones 2009, eq 9

        :param donor_catchments: Catchments to use as donors
        :type donor_catchments: list of :class:`Catchment`
        :return: 2-Dimensional, symmetric covariance matrix
        :rtype: :class:`numpy.ndarray`
        """
        p = len(donor_catchments)
        sigma = np.empty((p, p))
        for i in range(p):
            beta_i = self._beta(donor_catchments[i])
            n_i = donor_catchments[i].amax_records_end() - donor_catchments[i].amax_records_start() + 1
            for j in range(p):
                beta_j = self._beta(donor_catchments[j])
                n_j = donor_catchments[j].amax_records_end() - donor_catchments[j].amax_records_start() + 1
                rho_ij = self._lnqmed_corr(donor_catchments[i], donor_catchments[j])
                n_ij = min(donor_catchments[i].amax_records_end(), donor_catchments[j].amax_records_end()) - \
                       max(donor_catchments[i].amax_records_start(), donor_catchments[j].amax_records_start()) + 1
                sigma[i, j] = 4 * beta_i * beta_j * n_ij / n_i / n_j * rho_ij
        return sigma