docs: def computePerturbedExpectation(self, u_n, A_n, compute_uncertainty=True, uncertainty_method=None, warning_cutoff=1.0e-10, return_theta=False):
          """Compute the expectation of an observable of phase space function A(x) for a single new state.
  
          Parameters
          ----------
          u_n : np.ndarray, float, shape=(K, N_max)
              u_n[n] = u(x_n) - the energy of the new state at all N samples previously sampled.
          A_n : np.ndarray, float, shape=(K, N_max)
              A_n[n] = A(x_n) - the phase space function of the new state at all N samples previously sampled.  If this does NOT depend on state (e.g. position), it's simply the value of the observation.  If it DOES depend on the current state, then the observables from the previous states need to be reevaluated at THIS state.
          compute_uncertainty : bool, optional
              If False, the uncertainties will not be computed (default: True)
          uncertainty_method : string, optional
              Choice of method used to compute asymptotic covariance method, or None to use default
              See help for computeAsymptoticCovarianceMatrix() for more information on various methods. (default: None)
          warning_cutoff : float, optional
              Warn if squared-uncertainty is negative and larger in magnitude than this number (default: 1.0e-10)
          return_theta : bool, optional
              Whether or not to return the theta matrix.  Can be useful for complicated differences.
  
  
          Returns
          -------
          A : float
              A is the estimate for the expectation of A(x) for the specified state
          dA : float
              dA is uncertainty estimate for A
  
          Notes
          -----
          See Section IV of [1].
          # Compute estimators and uncertainty.
          #A = sum(W_n[:,K] * A_n[:]) # Eq. 15 of [1]
          #dA = abs(A) * np.sqrt(Theta_ij[K,K] + Theta_ij[K+1,K+1] - 2.0 * Theta_ij[K,K+1]) # Eq. 16 of [1]
          """
  
          if len(np.shape(u_n)) == 2:
              u_n = kn_to_n(u_n, N_k=self.N_k)
  
          if len(np.shape(A_n)) == 2:
              A_n = kn_to_n(A_n, N_k=self.N_k)
  
          # Convert to np matrix.
          A_n = np.array(A_n, dtype=np.float64)
  
          # Retrieve N and K for convenience.
          N = self.N
          K = self.K
  
          # Make A_k all positive so we can operate logarithmically for
          # robustness
          A_min = np.min(A_n)
          A_n = A_n - (A_min - 1)
  
          # Augment W_nk, N_k, and c_k for q_A(x) for the observable, with one
          # extra row/column for the specified state (Eq. 13 of [1]).
          # weight matrix
          Log_W_nk = np.zeros([N, K + 2], dtype=np.float64)
          N_k = np.zeros([K + 2], dtype=np.int32)  # counts
          f_k = np.zeros([K + 2], dtype=np.float64)  # free energies
  
          # Fill in first K states with existing q_k(x) from states.
          Log_W_nk[:, 0:K] = self.Log_W_nk
          N_k[0:K] = self.N_k
  
          # compute the free energy of the additional state
          log_w_n = self._computeUnnormalizedLogWeights(u_n)
          # Compute free energies
          f_k[K] = -_logsum(log_w_n)
          Log_W_nk[:, K] = log_w_n + f_k[K]
  
          # compute the observable at this state
          Log_W_nk[:, K + 1] = np.log(A_n) + Log_W_nk[:, K]
          f_k[K + 1] = -_logsum(Log_W_nk[:, K + 1])
          Log_W_nk[:, K + 1] += f_k[K + 1]              # normalize the row
          A = np.exp(-f_k[K + 1])
  
          if (compute_uncertainty or return_theta):
              # Compute augmented asymptotic covariance matrix.
              Theta_ij = self._computeAsymptoticCovarianceMatrix(
                  np.exp(Log_W_nk), N_k, method=uncertainty_method)
  
          if (compute_uncertainty):
              dA = np.abs(A) * np.sqrt(
                  Theta_ij[K + 1, K + 1] + Theta_ij[K, K] - 2.0 * Theta_ij[K, K + 1])  # Eq. 16 of [1]
  
          # shift answers back with the offset now that variances are computed
          A += (A_min - 1)
  
          returns = []
          returns.append(A)
  
          if (compute_uncertainty):
              returns.append(dA)
  
          if (return_theta):
              returns.append(Theta_ij)
  
          # Return expectations and uncertainties.
          return returns