@classmethod
def _covariance(cls, Phi):
    """Returns the estimated covariance of a set of summaries, with the prefactor for a virtual signal field. See equations (13)-(14) in |Leclercqetal2019|_.

        Parameters
        ----------
        Phi : array, double, dimension=(N,P)
            a set of N simulations

        Returns
        -------
        Sigma^hat_theta' : array, double, dimension=(P,P)
            (N+1)/N * Sigma^hat_theta where Sigma^hat_theta is the unbiased covariance estimator

        """
    import numpy as np
    from ctypes import c_double
    N = Phi.shape[0]
    covariance_prefactor = (N + 1) / N
    return covariance_prefactor * np.cov(Phi.T, ddof=1).astype(c_double)
