docs: def mbar_hessian(u_kn, N_k, f_k):
      """Hessian of MBAR objective function.
  
      Parameters
      ----------
      u_kn : np.ndarray, shape=(n_states, n_samples), dtype='float'
          The reduced potential energies, i.e. -log unnormalized probabilities
      N_k : np.ndarray, shape=(n_states), dtype='int'
          The number of samples in each state
      f_k : np.ndarray, shape=(n_states), dtype='float'
          The reduced free energies of each state
  
      Returns
      -------
      H : np.ndarray, dtype=float, shape=(n_states, n_states)
          Hessian of mbar objective function.
  
      Notes
      -----
      Equation (C9) in JCP MBAR paper.
      """
      u_kn, N_k, f_k = validate_inputs(u_kn, N_k, f_k)
  
      W = mbar_W_nk(u_kn, N_k, f_k)
  
      H = W.T.dot(W)
      H *= N_k
      H *= N_k[:, np.newaxis]
      H -= np.diag(W.sum(0) * N_k)
  
      return -1.0 * H