docs: def computeMultipleExpectations(self, A_in, u_n, compute_uncertainty=True,
                                      uncertainty_method=None, warning_cutoff=1.0e-10, return_theta=False):
          """Compute the expectations of multiple observables of phase space functions.
  
          Compute the expectations of multiple observables of phase space functions.
          [A_0(x),A_1(x),...,A_i(x)] along with the covariances of their estimates.  The state is specified by
          the choice of u_n, which is the energy of the n samples evaluated at a the chosen state.
  
          Parameters
          ----------
          A_in : np.ndarray, float, shape=(I, k, N)
              A_in[i,n] = A_i(x_n), the value of phase observable i for configuration n at state of interest
          u_n : np.ndarray, float, shape=(N)
              u_n[n] is the reduced potential of configuration n at the state of interest
          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_i : np.ndarray, float, shape=(I)
              A_i[i] is the estimate for the expectation of A_i(x) at the state specified by u_kn
          d2A_ij : np.ndarray, float, shape=(I, I)
              d2A_ij[i,j] is the COVARIANCE in the estimates of A_i[i] and A_i[j],
              not the square root of the covariance
  
          Notes
          -----
          Not fully tested!
  
          Examples
          --------
  
          >>> from pymbar import testsystems
          >>> (x_n, u_kn, N_k, s_n) = testsystems.HarmonicOscillatorsTestCase().sample(mode='u_kn')
          >>> mbar = MBAR(u_kn, N_k)
          >>> A_in = np.array([x_n,x_n**2,x_n**3])
          >>> u_n = u_kn[0,:]
          >>> [A_i, d2A_ij] = mbar.computeMultipleExpectations(A_in, u_kn)
  
          """
  
          # Retrieve N and K for convenience.
          I = A_in.shape[0]  # number of observables
          K = self.K
          N = self.N  # N is total number of samples
  
          if len(np.shape(A_in)) == 3:
              A_in_old = A_in.copy()
              A_in = np.zeros([I, N], np.float64)
              for i in range(I):
                  A_in[i,:] = kn_to_n(A_in_old[i, :, :], N_k=self.N_k)
  
          A_min = np.zeros([I], dtype=np.float64)
          for i in range(I):
              A_min[i] = np.min(A_in[i, :]) #find the minimum
              A_in[i, :] -= (A_min[i]-1)  #all now values will be positive so that we can work in logarithmic scale
  
          if len(np.shape(u_n)) == 2:
              u_n = kn_to_n(u_n, N_k=self.N_k)
  
          # Augment W_nk, N_k, and c_k for q_A(x) for the observables, with one
          # row for the specified state and I rows for the observable at that
          # state.
          # log weight matrix
          Log_W_nk = np.zeros([N, K + 1 + I], np.float64)
          W_nk = np.zeros([N, K + 1 + I], np.float64)  # weight matrix
          N_k = np.zeros([K + 1 + I], np.int32)  # counts
          f_k = np.zeros([K + 1 + I], np.float64)  # free energies
  
          # Fill in first section of matrix with existing q_k(x) from states.
          Log_W_nk[:, 0:K] = self.Log_W_nk
          W_nk[:, 0:K] = np.exp(self.Log_W_nk)
          N_k[0:K] = self.N_k
          f_k[0:K] = self.f_k
  
          # Compute row of W matrix for the extra state corresponding to u_kn.
          Log_W_nk[:, K]  = self._computeUnnormalizedLogWeights(u_n)
          f_k[K] = -_logsum(Log_W_nk[:, K])
          Log_W_nk[:, K] += f_k[K]
  
          # Compute the remaining rows/columns of W_nk and c_k for the
          # observables.
          for i in range(I):
              Log_W_nk[:, K+1+i] = np.log(A_in[i, :]) + Log_W_nk[:, K]
              f_k[K + 1 + i] = -_logsum(Log_W_nk[:, K + 1 + i])
              Log_W_nk[:, K + 1 + i] += f_k[K + 1 + i]    # normalize this row
  
          # Compute estimates.
          A_i = np.zeros([I], np.float64)
          for i in range(I):
              A_i[i] = np.exp(-f_k[K + 1 + i])
  
          if compute_uncertainty or return_theta:
              # Compute augmented asymptotic covariance matrix.
              W_nk = np.exp(Log_W_nk)
              Theta_ij = self._computeAsymptoticCovarianceMatrix(
                  W_nk, N_k, method=uncertainty_method)
  
          if compute_uncertainty:
              # Compute estimates of statistical covariance
              # these variances will be the same whether or not we subtract a different constant from each A_i
              # todo: vectorize
              # compute the covariance component without doing the double loop
              d2A_ij = np.zeros([I, I], np.float64)
              for i in range(I):
                  for j in range(I):
                      d2A_ij[i, j] = A_i[i] * A_i[j] * (Theta_ij[K + 1 + i, K + 1 + j] - Theta_ij[
                                                        K + 1 + i, K] - Theta_ij[K, K + 1 + j] + Theta_ij[K, K])
  
          # Now that covariances are computed, add the constants back to A_i that
          # were required to enforce positivity
          A_i = A_i + (A_min - 1)
  
          returns = []
          returns.append(A_i)
  
          if compute_uncertainty:
              returns.append(d2A_ij)
  
          if return_theta:
              returns.append(Theta_ij)
  
          # Return expectations and uncertainties.
          return returns