Merge pull request #2867 from junpenglao/fix_2866

fix BinaryGibbsMetropolis issue for p=.5

Author: Junpeng Lao

pymc3/step_methods/metropolis.py

@@ -297,13 +297,31 @@ def competence(var):
 
 
 class BinaryGibbsMetropolis(ArrayStep):
-    """A Metropolis-within-Gibbs step method optimized for binary variables"""
+    """A Metropolis-within-Gibbs step method optimized for binary variables
+
+    Parameters
+    ----------
+    vars : list
+        List of variables for sampler
+    order : list or 'random'
+        List of integers indicating the Gibbs update order
+        e.g., [0, 2, 1, ...]. Default is random
+    transit_p : float
+        The diagonal of the transition kernel. A value > .5 gives anticorrelated proposals,
+        which resulting in more efficient antithetical sampling.
+    model : PyMC Model
+        Optional model for sampling step. Defaults to None (taken from context).
+
+    """
     name = 'binary_gibbs_metropolis'
 
-    def __init__(self, vars, order='random', model=None):
+    def __init__(self, vars, order='random', transit_p=.8, model=None):
 
         model = pm.modelcontext(model)
 
+        # transition probabilities
+        self.transit_p = transit_p
+
         self.dim = sum(v.dsize for v in vars)
 
         if order == 'random':
@@ -330,11 +348,14 @@ def astep(self, q0, logp):
         logp_curr = logp(q)
 
         for idx in order:
-            curr_val, q[idx] = q[idx], True - q[idx]
-            logp_prop = logp(q)
-            q[idx], accepted = metrop_select(logp_prop - logp_curr, q[idx], curr_val)
-            if accepted:
-                logp_curr = logp_prop
+            # No need to do metropolis update if the same value is proposed,
+            # as you will get the same value regardless of accepted or reject
+            if nr.rand() < self.transit_p:
+                curr_val, q[idx] = q[idx], True - q[idx]
+                logp_prop = logp(q)
+                q[idx], accepted = metrop_select(logp_prop - logp_curr, q[idx], curr_val)
+                if accepted:
+                    logp_curr = logp_prop
 
         return q