def estimate_multivariate_ess_sigma(samples, batch_size):
    r"""Calculates the Sigma matrix which is part of the multivariate ESS calculation.

    This implementation is based on the Matlab implementation found at: https://github.com/lacerbi/multiESS

    The Sigma matrix is defined as:

    .. math::

        \Sigma = \Lambda + 2 * \sum_{k=1}^{\infty}{Cov(Y_{1}, Y_{1+k})}

    Where :math:`Y` are our samples and :math:`\Lambda` is the covariance matrix of the samples.

    This implementation computes the :math:`\Sigma` matrix using a Batch Mean estimator using the given batch size.
    The batch size has to be :math:`1 \le b_n \le n` and a typical value is either :math:`\lfloor n^{1/2} \rfloor`
    for slow mixing chains or :math:`\lfloor n^{1/3} \rfloor` for reasonable mixing chains.

    If the length of the chain is longer than the sum of the length of all the batches, this implementation
    calculates :math:`\Sigma` for every offset and returns the average of those offsets.

    Args:
        samples (ndarray): the samples for which we compute the sigma matrix. Expects an (p, n) array with
            p the number of parameters and n the sample size
        batch_size (int): the batch size used in the approximation of the correlation covariance

    Returns:
        ndarray: an pxp array with p the number of parameters in the samples.

    References:
        Vats D, Flegal J, Jones G (2016). Multivariate Output Analysis for Markov Chain Monte Carlo.
        arXiv:1512.07713v2 [math.ST]
    """
    sample_means = np.mean(samples, axis=1, dtype=np.float64)
    nmr_params, chain_length = samples.shape

    nmr_batches = int(np.floor(chain_length / batch_size))
    sigma = np.zeros((nmr_params, nmr_params))

    nmr_offsets = chain_length - nmr_batches * batch_size + 1

    for offset in range(nmr_offsets):
        batches = np.reshape(samples[:, np.array(offset + np.arange(0, nmr_batches * batch_size), dtype=np.int)].T,
                             [batch_size, nmr_batches, nmr_params], order='F')

        batch_means = np.squeeze(np.mean(batches, axis=0, dtype=np.float64))

        Z = batch_means - sample_means

        for x, y in itertools.product(range(nmr_params), range(nmr_params)):
            sigma[x, y] += np.sum(Z[:, x] * Z[:, y])

    return sigma * batch_size / (nmr_batches - 1) / nmr_offsets
r"""Calculates the Sigma matrix which is part of the multivariate ESS calculation.

    This implementation is based on the Matlab implementation found at: https://github.com/lacerbi/multiESS

    The Sigma matrix is defined as:

    .. math::

        \Sigma = \Lambda + 2 * \sum_{k=1}^{\infty}{Cov(Y_{1}, Y_{1+k})}

    Where :math:`Y` are our samples and :math:`\Lambda` is the covariance matrix of the samples.

    This implementation computes the :math:`\Sigma` matrix using a Batch Mean estimator using the given batch size.
    The batch size has to be :math:`1 \le b_n \le n` and a typical value is either :math:`\lfloor n^{1/2} \rfloor`
    for slow mixing chains or :math:`\lfloor n^{1/3} \rfloor` for reasonable mixing chains.

    If the length of the chain is longer than the sum of the length of all the batches, this implementation
    calculates :math:`\Sigma` for every offset and returns the average of those offsets.

    Args:
        samples (ndarray): the samples for which we compute the sigma matrix. Expects an (p, n) array with
            p the number of parameters and n the sample size
        batch_size (int): the batch size used in the approximation of the correlation covariance

    Returns:
        ndarray: an pxp array with p the number of parameters in the samples.

    References:
        Vats D, Flegal J, Jones G (2016). Multivariate Output Analysis for Markov Chain Monte Carlo.
        arXiv:1512.07713v2 [math.ST]