@classmethod
def _inverse_covariance(cls, Phi, covariance, EPS_Sigma=1e-07, EPS_residual
    =0.001):
    """Computes the unbiased inverse of an estimated covariance matrix, using the `Hartlap, Simon & Schneider (2007) <https://arxiv.org/abs/astro-ph/0608064>`__ correction factor. See equations (13) and (15) in |Leclercqetal2019|_.

        Parameters
        ----------
        Phi : array, double, dimension=(N,P)
            a set of N simulations
        covariance : array, double, dimension=(P,P)
            estimated covariance matrix
        EPS_Sigma : double, optional, default=1e-7
            epsilon to be added to the diagonal of the covariance matrix, if necessary
        EPS_residual : double, optional, default=1e-3
            maximum residual in covariance*covariance^{-1} before attempting regularization

        Returns
        -------
        Sigma^hat_theta^{-1}' : array, double, dimension=(P,P)
            alpha*(Sigma^hat_theta')^{-1} where (Sigma^hat_theta')^{-1} is the inverse sample covariance

        """
    from pyselfi.utils import regular_inv
    N = Phi.shape[0]
    P = Phi.shape[1]
    alpha = (N - P - 2.0) / (N - 1.0) if P < N - 2 else 1.0
    inv_covariance = regular_inv(covariance, EPS_Sigma, EPS_residual)
    return alpha * inv_covariance
