Repository: arviz-devs/arviz

Issue:
Incorrect denominator in bfmi
The original equation for `bfmi` is at the bottom of page 7 in [this paper](https://arxiv.org/pdf/1604.00695v1.pdf). The numerator in numpy for a vector `E` is `(length(E) - 1) * np.mean(np.square(np.diff(E)))`. The denominator is `(length(E) - 1) * np.var(E, ddof = 1)`. The two scaling factors cancel and we get
```python
def bfmi(E):
    return np.mean(np.square(np.diff(E))) / np.var(E, ddof = 1)
```
The current implementation uses `ddof=0`, so it's off from the equation in the paper by `length(E) / (length(E) - 1)`

**Additional context**
PyStan's implementation [here](https://github.com/stan-dev/pystan/blob/cbc5c0fc8597504776c08a138adc9ae63329e6a1/pystan/diagnostics.py#L213-L217) is equivalent to the paper's. RStan's implementation of `get_bfmi` [here](https://rdrr.io/cran/rstan/src/R/check_hmc_diagnostics.R) has a different implementation than either ArviZ's current one or PyStan's. Theirs amounts to the original paper's equation multiplied by `(length(E) - 1) / length(E)`.

Discussion:
nitishp25: can someone tell me why this fails?
percygautam: Hi @nitishp25 , I think you have to change the test_bfmi() function of test_diagnostics.py. Run the code in test function with updated ddof value and replace the result with the previous value.
nitishp25: Got it, thanks @percygautam!
nitishp25: Ready to merge.
OriolAbril: Also, this is a little unrelated but it would be great if you could fix it. I have just seen that [bfmi docs](https://arviz-devs.github.io/arviz/generated/arviz.bfmi.html#arviz.bfmi) are not rendered properly (see [rhat docs](https://arviz-devs.github.io/arviz/generated/arviz.rhat.html#arviz.rhat) to see how they should look like) there is no `Out` section even though there should be and the `In []` is not properly highlighted. A white line between [lines 48 and 49](https://github.com/arviz-devs/arviz/blob/c3158613f9e793902eddd0e5ecd715595c31c824/arviz/stats/diagnostics.py#L48) should fix this.
nitishp25: > Also, this is a little unrelated but it would be great if you could fix it. I have just seen that [bfmi docs](https://arviz-devs.github.io/arviz/generated/arviz.bfmi.html#arviz.bfmi) are not rendered properly (see [rhat docs](https://arviz-devs.github.io/arviz/generated/arviz.rhat.html#arviz.rhat) to see how they should look like) there is no `Out` section even though there should be and the `In []` is not properly highlighted. A white line between [lines 48 and 49](https://github.com/arviz-devs/arviz/blob/c3158613f9e793902eddd0e5ecd715595c31c824/arviz/stats/diagnostics.py#L48) should fix this.

Got it. Working on it.
sethaxen: This looks like the right fix. With the fix to the else statement suggested by @OriolAbril, this looks good to me. Thanks for the PR!
nitishp25: Please review @OriolAbril 

Patch:
diff --git a/arviz/stats/diagnostics.py b/arviz/stats/diagnostics.py
index eae8692..21577b6 100644
--- a/arviz/stats/diagnostics.py
+++ b/arviz/stats/diagnostics.py
@@ -46,6 +46,7 @@ def bfmi(data):
     Examples
     --------
     Compute the BFMI of an InferenceData object
+
     .. ipython::
 
         In [1]: import arviz as az
@@ -526,9 +527,9 @@ def _bfmi(energy):
     energy_mat = np.atleast_2d(energy)
     num = np.square(np.diff(energy_mat, axis=1)).mean(axis=1)  # pylint: disable=no-member
     if energy_mat.ndim == 2:
-        den = _numba_var(svar, np.var, energy_mat, axis=1, ddof=0)
+        den = _numba_var(svar, np.var, energy_mat, axis=1, ddof=1)
     else:
-        den = np.var(energy, axis=1)
+        den = np.var(energy, axis=1, ddof=1)
     return num / den
 
 


Test patch:
diff --git a/arviz/tests/test_diagnostics.py b/arviz/tests/test_diagnostics.py
index 1ba514a..4100cdb 100644
--- a/arviz/tests/test_diagnostics.py
+++ b/arviz/tests/test_diagnostics.py
@@ -45,7 +45,7 @@ def data():
 class TestDiagnostics:
     def test_bfmi(self):
         energy = np.array([1, 2, 3, 4])
-        assert_almost_equal(bfmi(energy), 0.8)
+        assert_almost_equal(bfmi(energy), 0.6)
 
     def test_bfmi_dataset(self):
         data = load_arviz_data("centered_eight")