docs: def computeMultipleExpectations(self, A_in, u_n, compute_uncertainty=True, compute_covariance=False,
                                      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)] at a single state,
          along with the error in the estimates and the uncertainty in
          the 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
          compute_uncertainty : bool, optional, default=True
              If True, calculate the uncertainty
          compute_covariance : bool, optional, default=False
              If True, calculate the covariance
          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)
  
          Returns
          -------------
          result_vals : dictionary
  
          Possible keys in the result_vals dictionary:
  
          'mu' : np.ndarray, float, shape=(I)
              result_vals['mu'] is the estimate for the expectation of A_i(x) at the state specified by u_kn
          'sigma' : np.ndarray, float, shape = (I)
              result_vals['sigma'] is the uncertainty in the expectation of A_state_map[i](x) at the state specified by u_n[state_map[i],:]
              or None if compute_uncertainty is False
          'covariances' : np.ndarray, float, shape=(I, I)
              result_vals['covariances'] is the COVARIANCE in the estimates of A_i[i] and A_i[j]: we can't actually take a square root
              or None if compute_covariance is False
          'Theta': np.ndarray, float, shape=(I, I), covariances of the log weights, useful for some additional calculations.
  
          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,:]
          >>> results = 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()  # convert to k by n format
              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)
  
          if len(np.shape(u_n)) == 2:
              u_n = kn_to_n(u_n, N_k = self.N_k)
  
          state_map = np.zeros([2,I],int)
          state_map[1,:] = np.arange(I)  # same (first) state for all variables.
  
          inner_results = self.computeExpectationsInner(A_in,u_n,state_map,
                                                        return_theta=(compute_uncertainty or compute_covariance),
                                                        uncertainty_method=uncertainty_method,
                                                        warning_cutoff=warning_cutoff)
          result_vals = dict()
          expectations, uncertainties, covariances = None, None, None
          result_vals['mu'] = inner_results['observables']
  
          if compute_uncertainty or compute_covariance or return_theta:
              Adiag = np.zeros([2*I,2*I],dtype=np.float64)
              diag = np.ones(2*I,dtype=np.float64)
              diag[0:I] = diag[I:2*I] = inner_results['observables']-inner_results['Amin']
              np.fill_diagonal(Adiag,diag)
              Theta = Adiag*inner_results['Theta']*Adiag
              if return_theta:
                  result_vals['Theta'] = Theta
  
              if compute_uncertainty:
                  covA_ij = np.array(Theta[0:I,0:I]+Theta[I:2*I,I:2*I]-Theta[0:I,I:2*I]-Theta[I:2*I,0:I])
                  result_vals['sigma'] = np.sqrt(covA_ij[0:I,0:I].diagonal())
  
              if compute_covariance:
                  # compute estimate of statistical covariance of the observables
                  result_vals['covariances'] = inner_results['Theta'][0:I,0:I]
  
          return result_vals