Merge pull request #518 from CUQI-DTU/fix_invalid_proposal_accepted

Fix issue of invalid proposals being accepted

cuqi/experimental/mcmc/_cwmh.py

@@ -158,7 +158,9 @@ def step(self):
 
             # accept/reject
             u_theta = np.log(np.random.rand())
-            if (u_theta <= alpha): # accept
+            if (u_theta <= alpha) and \
+               (not np.isnan(target_eval_star)) and \
+               (not np.isinf(target_eval_star)):
                 x_t[j] = x_all_components[j]
                 target_eval_t = target_eval_star
                 acc[j] = 1

cuqi/experimental/mcmc/_hmc.py

@@ -256,7 +256,10 @@ def step(self):
 
             # Metropolis step
             alpha2 = min(1, (n_prime/n)) #min(0, np.log(n_p) - np.log(n))
-            if (s_prime == 1) and (np.random.rand() <= alpha2):
+            if (s_prime == 1) and \
+                (np.random.rand() <= alpha2) and \
+                (not np.isnan(logd_prime)) and \
+                (not np.isinf(logd_prime)):
                 self.current_point = point_prime
                 self.current_target_logd = logd_prime
                 self.current_target_grad = np.copy(grad_prime)

cuqi/experimental/mcmc/_langevin_algorithm.py

@@ -102,10 +102,14 @@ def _accept_or_reject(self, x_star, target_eval_star, target_grad_star):
         scalar
             1 (accepted)
         """
-        self.current_point = x_star
-        self.current_target_logd = target_eval_star
-        self.current_target_grad = target_grad_star
-        acc = 1
+        acc = 0
+        if (not np.isnan(target_eval_star)) and \
+           (not np.isinf(target_eval_star)):
+            self.current_point = x_star
+            self.current_target_logd = target_eval_star
+            self.current_target_grad = target_grad_star
+            acc = 1
+
         return acc
 
     def step(self):
@@ -211,7 +215,9 @@ def _accept_or_reject(self, x_star, target_eval_star, target_grad_star):
         # accept/reject with Metropolis
         acc = 0
         log_u = np.log(np.random.rand())
-        if (log_u <= log_alpha) and (np.isnan(target_eval_star) == False):
+        if (log_u <= log_alpha) and \
+           (not np.isnan(target_eval_star)) and \
+           (not np.isinf(target_eval_star)):
             self.current_point = x_star
             self.current_target_logd = target_eval_star
             self.current_target_grad = target_grad_star

cuqi/experimental/mcmc/_mh.py

@@ -57,7 +57,9 @@ def step(self):
         # accept/reject
         u_theta = np.log(np.random.rand())
         acc = 0
-        if (u_theta <= alpha):
+        if (u_theta <= alpha) and \
+           (not np.isnan(target_eval_star)) and \
+           (not np.isinf(target_eval_star)):
             self.current_point = x_star
             self.current_target_logd = target_eval_star
             acc = 1

tests/zexperimental/test_mcmc.py

@@ -903,4 +903,59 @@ def test_HybridGibbs_updates_state_only_after_accepting_sample():
         # In case sampler has accepted samples, the state should have been updated
         else:
             assert sampler_states[key] != new_state, f"Sampler {key} state was erroneously not updated in Gibbs scheme, even when new samples were accepted. State: \n {new_state}"
-
+
+sampler_instances_for_bounded_distribution = [
+    cuqi.experimental.mcmc.MH(
+        target=cuqi.distribution.Beta(0.5, 0.5),
+        initial_point=np.array([0.1]),
+        scale=0.1,
+    ),
+    cuqi.experimental.mcmc.ULA(
+        cuqi.distribution.Beta(0.5, 0.5),
+        initial_point=np.array([0.1]),
+        scale=0.1,
+    ),
+    cuqi.experimental.mcmc.MALA(
+        cuqi.distribution.Beta(0.5, 0.5),
+        initial_point=np.array([0.1]),
+        scale=0.1,
+    ),
+    cuqi.experimental.mcmc.NUTS(
+        cuqi.distribution.Beta(0.5, 0.5), initial_point=np.array([0.1])
+    ),
+    cuqi.experimental.mcmc.MH(
+        target=cuqi.distribution.Beta(np.array([0.5, 0.5]), np.array([0.5, 0.5])),
+        initial_point=np.array([0.1, 0.1]),
+        scale=0.1,
+    ),
+    cuqi.experimental.mcmc.CWMH(
+        target=cuqi.distribution.Beta(
+            np.array([0.5, 0.5]), np.array([0.5, 0.5])
+        ),
+        initial_point=np.array([0.1, 0.1]),
+        scale=0.1,
+    ),
+    cuqi.experimental.mcmc.ULA(
+        cuqi.distribution.Beta(np.array([0.5, 0.5]), np.array([0.5, 0.5])),
+        initial_point=np.array([0.1, 0.1]),
+        scale=0.1,
+    ),
+    cuqi.experimental.mcmc.MALA(
+        cuqi.distribution.Beta(np.array([0.5, 0.5]), np.array([0.5, 0.5])),
+        initial_point=np.array([0.1, 0.1]),
+        scale=0.1,
+    ),
+    cuqi.experimental.mcmc.NUTS(
+        cuqi.distribution.Beta(np.array([0.5, 0.5]), np.array([0.5, 0.5])),
+        initial_point=np.array([0.1, 0.1]),
+    ),
+]
+
+@pytest.mark.parametrize("sampler", sampler_instances_for_bounded_distribution)
+def test_if_invalid_sample_accepted(sampler: cuqi.experimental.mcmc.Sampler):
+    sampler.sample(50)
+    samples = sampler.get_samples().samples
+    tol = 1e-8
+    assert (
+        samples.min() > 0.0 - tol and samples.max() < 1.0 + tol
+    ), f"Invalid samples accepted for sampler {sampler.__class__.__name__}."