docs: 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