Revert divergences PR

Author: eigenfoo

README.md

@@ -16,11 +16,11 @@ Currently Implemented
 - Leapfrog integrator
 - Hamiltonian Monte Carlo
 - Some log probabilities (normal, multivariate normal, mixtures, funnel)
-- Divergences
 
 Roadmap
 -------
 
+- Divergences
 - Mass matrix adaptation
 - Diagnostics
 - [NUTS](https://arxiv.org/abs/1111.4246)

minimc/minimc.py

@@ -76,18 +76,14 @@ def hamiltonian_monte_carlo(
         start_log_p = np.sum(momentum.logpdf(p0)) - initial_potential
         new_log_p = np.sum(momentum.logpdf(p_new)) - final_V
         energy_change = new_log_p - start_log_p
-        p_accept = min(1, np.exp(energy_change))
 
-        if abs(energy_change) < max_energy_change:
-            # Check Metropolis acceptance criterion
-            if np.random.rand() < p_accept:
-                samples.append(q_new)
-                initial_potential = final_V
-                initial_potential_grad = final_dVdq
-            else:
-                samples.append(np.copy(samples[-1]))
+        # Check Metropolis acceptance criterion
+        p_accept = min(1, np.exp(energy_change))
+        if np.random.rand() < p_accept:
+            samples.append(q_new)
+            initial_potential = final_V
+            initial_potential_grad = final_dVdq
         else:
-            # Divergence encountered
             samples.append(np.copy(samples[-1]))
 
         if idx < tune - 1:

minimc/minimc_slow.py

@@ -76,18 +76,14 @@ def hamiltonian_monte_carlo(
         energy_change = start_log_p - new_log_p
         p_accept = np.exp(energy_change)
 
-        if abs(energy_change) < max_energy_change:
-            if np.random.rand() < p_accept:
-                samples.append(q_new)
-                accepted.append(True)
-            else:
-                if do_reject:
-                    samples.append(np.copy(samples[-1]))
-                else:
-                    samples.append(q_new)
-                accepted.append(False)
+        if np.random.rand() < p_accept:
+            samples.append(q_new)
+            accepted.append(True)
         else:
-            samples.append(np.copy(samples[-1]))
+            if do_reject:
+                samples.append(np.copy(samples[-1]))
+            else:
+                samples.append(q_new)
             accepted.append(False)
         p_accepts.append(p_accept)
 

test/test_minimc.py

@@ -79,4 +79,4 @@ def test_hamiltonian_monte_carlo_mv():
     )
     assert samples.shape[0] == 100
     assert_allclose(mu, np.mean(samples, axis=0), atol=0.3)
-    assert_allclose(cov, np.cov(samples.T), atol=0.3)
+    assert_allclose(cov, np.cov(samples.T), atol=0.5)