code: def bfmi(energy):
      r
      energy_mat = np.atleast_2d(energy)
      num = np.square(np.diff(energy_mat, axis=1)).mean(axis=1)  # pylint: disable=no-member
      den = np.var(energy_mat, axis=1)
      return num / den
docs: r"""Calculate the estimated Bayesian fraction of missing information (BFMI).
  
      BFMI quantifies how well momentum resampling matches the marginal energy distribution. For more
      information on BFMI, see https://arxiv.org/pdf/1604.00695v1.pdf. The current advice is that
      values smaller than 0.3 indicate poor sampling. However, this threshold is provisional and may
      change. See http://mc-stan.org/users/documentation/case-studies/pystan_workflow.html for more
      information.
  
      Parameters
      ----------
      energy : NumPy array
          Should be extracted from a gradient based sampler, such as in Stan or PyMC3. Typically,
          after converting a trace or fit to InferenceData, the energy will be in
          `data.sample_stats.energy`.
  
      Returns
      -------
      z : array
          The Bayesian fraction of missing information of the model and trace. One element per
          chain in the trace.