Re-implemented tuning to keep logic in step function

Author: Christopher Fonnesbeck

examples/lasso_block_update.py

@@ -4,15 +4,15 @@
 # <markdowncell>
 
 # Sometimes, it is very useful to update a set of parameters together. For example, variables that are highly correlated are often good to update together. In PyMC 3 block updating is simple, as example will demonstrate.
-# 
+#
 # Here we have a LASSO regression model where the two coefficients are strongly correlated. Normally, we would define the coefficient parameters as a single random variable, but here we define them separately to show how to do block updates.
-# 
+#
 # First we generate some fake data.
 
 # <codecell>
 from matplotlib.pylab import *
-from pymc import * 
-import numpy as np 
+from pymc import *
+import numpy as np
 
 d = np.random.normal(size=(3, 30))
 d1 = d[0] + 4
@@ -29,24 +29,24 @@
     s = Exponential('s', 1)
     m1 = Laplace('m1', 0, 100)
     m2 = Laplace('m2', 0, 100)
-    
+
     p = d1*m1 + d2*m2
-    
-    y = Normal('y', p, s**-2, observed = yd) 
+
+    y = Normal('y', p, s**-2, observed = yd)
 
 # <markdowncell>
 
-# For most samplers, including Metropolis and HamiltonianMC, simply pass a list of variables to sample as a block. This works with both scalar and array parameters. 
+# For most samplers, including Metropolis and HamiltonianMC, simply pass a list of variables to sample as a block. This works with both scalar and array parameters.
 
 # <codecell>
 
-with model: 
+with model:
     start = find_MAP()
-    
-    step1 = Metropolis([m1, m2], approx_hess(start, [m1,m2]))
-    
-    step2 = Metropolis([s], approx_hess(start, [s]))
-    
+
+    step1 = Metropolis([m1, m2])
+
+    step2 = Metropolis([s], laplace_proposal)
+
     trace = sample(5000, [step1, step2], start)
 
 # <codecell>

pymc/sample.py

@@ -8,8 +8,7 @@
 
 __all__ = ['sample', 'psample']
 
-def sample(draws, step, start=None, trace=None, track_progress=True,
-        tune_interval=100, model=None):
+def sample(draws, step, start=None, trace=None, track_progress=True, model=None):
     """
     Draw a number of samples using the given step method.
     Multiple step methods supported via compound step method
@@ -58,58 +57,11 @@ def sample(draws, step, start=None, trace=None, track_progress=True,
     for i in xrange(draws):
         point = step.step(point)
         trace = trace.record(point)
-        if i and not (i % tune_interval) and step.tune:
-            step = tune(step, tune_interval)
         if track_progress:
             progress.update(i)
 
     return trace
 
-
-def tune(step, tune_interval):
-    """
-    Tunes the scaling parameter for the proposal distribution
-    according to the acceptance rate over the last tune_interval:
-
-    Rate    Variance adaptation
-    ----    -------------------
-    <0.001        x 0.1
-    <0.05         x 0.5
-    <0.2          x 0.9
-    >0.5          x 1.1
-    >0.75         x 2
-    >0.95         x 10
-
-    """
-
-    # Calculate acceptance rate
-    acc_rate = step.accepted / float(tune_interval)
-
-    # Switch statement
-    if acc_rate<0.001:
-        # reduce by 90 percent
-        step.scaling *= 0.1
-    elif acc_rate<0.05:
-        # reduce by 50 percent
-        step.scaling *= 0.5
-    elif acc_rate<0.2:
-        # reduce by ten percent
-        step.scaling *= 0.9
-    elif acc_rate>0.95:
-        # increase by factor of ten
-        step.scaling *= 10.0
-    elif acc_rate>0.75:
-        # increase by double
-        step.scaling *= 2.0
-    elif acc_rate>0.5:
-        # increase by ten percent
-        step.scaling *= 1.1
-
-    # Re-initialize rejection count
-    step.accepted = 0
-
-    return step
-
 def argsample(args):
     """ defined at top level so it can be pickled"""
     return sample(*args)

pymc/step_methods/arraystep.py

@@ -7,11 +7,10 @@
 
 # TODO Add docstrings to ArrayStep
 class ArrayStep(object):
-    def __init__(self, vars, fs, allvars = False, tune=False):
+    def __init__(self, vars, fs, allvars = False):
         self.ordering = ArrayOrdering(vars)
         self.fs = fs
         self.allvars = allvars
-        self.tune = tune
 
     def step(self, point):
         bij = DictToArrayBijection(self.ordering, point)
@@ -29,10 +28,10 @@ def metrop_select(mr, q, q0):
     # Compare acceptance ratio to uniform random number
     if isfinite(mr) and log(uniform()) < mr:
         # Accept proposed value
-        return q, True
+        return q
     else:
         # Reject proposed value
-        return q0, False
+        return q0
 
 
 class SamplerHist(object):

pymc/step_methods/compound.py

@@ -8,7 +8,6 @@
 @quickclass(object)
 def CompoundStep(methods):
     methods = list(methods)
-    tune = False
     def step(point):
         for method in methods:
             point = method.step(point)

pymc/step_methods/metropolis.py

@@ -6,17 +6,17 @@
 from arraystep import *
 from numpy.random import normal, standard_cauchy, standard_exponential
 
-__all__ = ['Metropolis']
+__all__ = ['Metropolis', 'normal_proposal', 'cauchy_proposal', 'laplace_proposal']
 
 # Available proposal distributions for Metropolis
-def normal_proposal(s, n):
-    return lambda: normal(scale=s, size=n)
+def normal_proposal(s):
+    return lambda: normal(scale=s)
 
-def cauchy_proposal(s, n):
-    return lambda: standard_cauchy(size=n) * s
+def cauchy_proposal(s):
+    return lambda: standard_cauchy(size=np.size(s)) * s
 
-def laplace_proposal(s, n):
-    return lambda: (standard_exponential(size=n) - standard_exponential(size=n)) * s
+def laplace_proposal(s):
+    return lambda: (standard_exponential(size=np.size(s)) - standard_exponential(size=np.size(s))) * s
 
 
 class Metropolis(ArrayStep):
@@ -40,28 +40,77 @@ class Metropolis(ArrayStep):
         Optional model for sampling step. Defaults to None (taken from context).
 
     """
-    def __init__(self, vars, S, proposal_dist=quad_potential, scaling=1.,
-            tune=True, model=None):
+    def __init__(self, vars, S=None, proposal_dist=normal_proposal, scaling=1.,
+            tune=True, tune_interval=100, model=None):
 
         model = modelcontext(model)
 
-        try:
-            self.proposal_dist = proposal_dist(S, n=len(vars))
-        except TypeError:
-            # quadpotential does not require n
-            self.proposal_dist = proposal_dist(S)
+        if S is None:
+            S = [1]*len(vars)
+        self.proposal_dist = proposal_dist(S)
         self.scaling = scaling
+        self.tune = tune
+        self.tune_interval = tune_interval
+        self.steps_until_tune = tune_interval
         self.accepted = 0
-        super(Metropolis,self).__init__(vars, [model.logpc], tune=tune)
+        super(Metropolis,self).__init__(vars, [model.logpc])
 
     def astep(self, q0, logp):
 
+        if self.tune and not self.steps_until_tune:
+            # Tune scaling parameter
+            self.scaling = tune(self.scaling, self.accepted/float(self.tune_interval))
+            # Reset counter
+            self.steps_until_tune = self.tune_interval
+            self.accepted = 0
+
         delta = self.proposal_dist() * self.scaling
 
         q = q0 + delta
 
-        q_new, accepted = metrop_select(logp(q) - logp(q0), q, q0)
+        q_new = metrop_select(logp(q) - logp(q0), q, q0)
+
+        if (any(q_new!=q0) or all(q0==q)):
+            self.accepted += 1
 
-        self.accepted += int(accepted)
+        self.steps_until_tune -= 1
 
         return q_new
+
+def tune(scale, acc_rate):
+    """
+    Tunes the scaling parameter for the proposal distribution
+    according to the acceptance rate over the last tune_interval:
+
+    Rate    Variance adaptation
+    ----    -------------------
+    <0.001        x 0.1
+    <0.05         x 0.5
+    <0.2          x 0.9
+    >0.5          x 1.1
+    >0.75         x 2
+    >0.95         x 10
+
+    """
+
+    # Switch statement
+    if acc_rate<0.001:
+        # reduce by 90 percent
+        scale *= 0.1
+    elif acc_rate<0.05:
+        # reduce by 50 percent
+        scale *= 0.5
+    elif acc_rate<0.2:
+        # reduce by ten percent
+        scale *= 0.9
+    elif acc_rate>0.95:
+        # increase by factor of ten
+        scale *= 10.0
+    elif acc_rate>0.75:
+        # increase by double
+        scale *= 2.0
+    elif acc_rate>0.5:
+        # increase by ten percent
+        scale *= 1.1
+
+    return scale