docs: def _computeAsymptoticCovarianceMatrix(self, W, N_k, method=None):
          """Compute estimate of the asymptotic covariance matrix.
  
          Parameters
          ----------
          W : np.ndarray, shape=(N, K), dtype='float'
              The normalized weight matrix for snapshots and states.
              W[n, k] is the weight of snapshot n in state k.
          N_k : np.ndarray, shape=(K), dtype='int'
              N_k[k] is the number of samples from state k.
          method : string, optional, default=None
              Method used to compute the asymptotic covariance matrix.
              Must be either "approximate", "svd", or "svd-ew".  If None,
              defaults to "svd-ew".
  
          Returns
          -------
          Theta: np.ndarray, shape=(K, K), dtype='float'
              Asymptotic covariance matrix
  
          Notes
          -----
          The computational costs of the various 'method' arguments varies:
            'svd' computes the generalized inverse using the singular value decomposition -- this should be efficient yet accurate (faster)
            'svd-ew' is the same as 'svd', but uses the eigenvalue decomposition of W'W to bypass the need to perform an SVD (fastest)
            'approximate' only requires multiplication of KxN and NxK matrices, but is an approximate underestimate of the uncertainty.
  
          svd and svd-ew are described in appendix D of Shirts, 2007 JCP, while
          "approximate" in Section 4 of Kong, 2003. J. R. Statist. Soc. B.
  
          We currently recommend 'svd-ew'.
          """
  
          # Set 'svd-ew' as default if uncertainty method specified as None.
          if method == None:
              method = 'svd-ew'
  
          # Get dimensions of weight matrix.
          [N, K] = W.shape
  
          # Check dimensions
          if(K != N_k.size):
              raise ParameterError(
                  'W must be NxK, where N_k is a K-dimensional array.')
          if(np.sum(N_k) != N):
              raise ParameterError('W must be NxK, where N = sum_k N_k.')
  
          check_w_normalized(W, N_k)
  
          # Compute estimate of asymptotic covariance matrix using specified method.
          if method == 'approximate':
              # Use fast approximate expression from Kong et al. -- this underestimates the true covariance, but may be a good approximation in some cases and requires no matrix inversions
              # Theta = P'P
  
              # Construct matrices
              W = np.matrix(W, dtype=np.float64)
  
              # Compute covariance
              Theta = W.T * W
  
          elif method == 'svd':
              # Use singular value decomposition based approach given in supplementary material to efficiently compute uncertainty
              # See Appendix D.1, Eq. D4 in [1].
  
              # Construct matrices
              Ndiag = np.matrix(np.diag(N_k), dtype=np.float64)
              W = np.matrix(W, dtype=np.float64)
              I = np.identity(K, dtype=np.float64)
  
              # Compute SVD of W
              [U, S, Vt] = linalg.svd(W, full_matrices=False)  # False Avoids O(N^2) memory allocation by only calculting the active subspace of U.
              Sigma = np.matrix(np.diag(S))
              V = np.matrix(Vt).T
  
              # Compute covariance
              Theta = V * Sigma * self._pseudoinverse(
                  I - Sigma * V.T * Ndiag * V * Sigma) * Sigma * V.T
  
          elif method == 'svd-ew':
              # Use singular value decomposition based approach given in supplementary material to efficiently compute uncertainty
              # The eigenvalue decomposition of W'W is used to forego computing the SVD.
              # See Appendix D.1, Eqs. D4 and D5 of [1].
  
              # Construct matrices
              Ndiag = np.matrix(np.diag(N_k), dtype=np.float64)
              W = np.matrix(W, dtype=np.float64)
              I = np.identity(K, dtype=np.float64)
  
              # Compute singular values and right singular vectors of W without using SVD
              # Instead, we compute eigenvalues and eigenvectors of W'W.
              # Note W'W = (U S V')'(U S V') = V S' U' U S V' = V (S'S) V'
              [S2, V] = linalg.eigh(W.T * W)
              # Set any slightly negative eigenvalues to zero.
              S2[np.where(S2 < 0.0)] = 0.0
              # Form matrix of singular values Sigma, and V.
              Sigma = np.matrix(np.diag(np.sqrt(S2)))
              V = np.matrix(V)
  
              # Compute covariance
              Theta = V * Sigma * self._pseudoinverse(
                  I - Sigma * V.T * Ndiag * V * Sigma) * Sigma * V.T
  
          else:
              # Raise an exception.
              raise ParameterError('Method ' + method + ' unrecognized.')
  
          return Theta