docs: def estimate_multivariate_ess(samples, batch_size_generator=None, full_output=False):
      r"""Compute the multivariate Effective Sample Size of your (single instance set of) samples.
  
      This multivariate ESS is defined in Vats et al. (2016) and is given by:
  
      .. math::
  
          ESS = n \bigg(\frac{|\Lambda|}{|\Sigma|}\bigg)^{1/p}
  
      Where :math:`n` is the number of samples, :math:`p` the number of parameters, :math:`\Lambda` is the covariance
      matrix of the parameters and :math:`\Sigma` captures the covariance structure in the target together with
      the covariance due to correlated samples. :math:`\Sigma` is estimated using
      :func:`estimate_multivariate_ess_sigma`.
  
      In the case of NaN in any part of the computation the ESS is set to 0.
  
      To compute the multivariate ESS for multiple problems, please use :func:`multivariate_ess`.
  
      Args:
          samples (ndarray): an pxn matrix with for p parameters and n samples.
          batch_size_generator (MultiVariateESSBatchSizeGenerator): the batch size generator, tells us how many
              batches and of which size we use for estimating the minimum ESS. Defaults to :class:`SquareRootSingleBatch`
          full_output (boolean): set to True to return the estimated :math:`\Sigma` and the optimal batch size.
  
      Returns:
          float or tuple: when full_output is set to True we return a tuple with the estimated multivariate ESS,
              the estimated :math:`\Sigma` matrix and the optimal batch size. When full_output is False (the default)
              we only return the ESS.
  
      References:
          Vats D, Flegal J, Jones G (2016). Multivariate Output Analysis for Markov Chain Monte Carlo.
          arXiv:1512.07713v2 [math.ST]
      """
      batch_size_generator = batch_size_generator or SquareRootSingleBatch()
  
      batch_sizes = batch_size_generator.get_multivariate_ess_batch_sizes(*samples.shape)
  
      nmr_params, chain_length = samples.shape
      nmr_batches = len(batch_sizes)
  
      det_lambda = det(np.cov(samples))
  
      ess_estimates = np.zeros(nmr_batches)
      sigma_estimates = np.zeros((nmr_params, nmr_params, nmr_batches))
  
      for i in range(0, nmr_batches):
          sigma = estimate_multivariate_ess_sigma(samples, int(batch_sizes[i]))
          ess = chain_length * (det_lambda**(1.0 / nmr_params) / det(sigma)**(1.0 / nmr_params))
  
          ess_estimates[i] = ess
          sigma_estimates[..., i] = sigma
  
      ess_estimates = np.nan_to_num(ess_estimates)
  
      if nmr_batches > 1:
          idx = np.argmin(ess_estimates)
      else:
          idx = 0
  
      if full_output:
          return ess_estimates[idx], sigma_estimates[..., idx], batch_sizes[idx]
      return ess_estimates[idx]