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
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]