Added tuning function for step functions

Author: Christopher Fonnesbeck

pymc/sample.py

@@ -8,41 +8,42 @@
 
 __all__ = ['sample', 'psample']
 
-def sample(draws, step, start = None, trace = None, track_progress = True, model = None): 
+def sample(draws, step, start=None, trace=None, track_progress=True,
+        tune_interval=100, model=None):
     """
-    Draw a number of samples using the given step method. 
-    Multiple step methods supported via compound step method 
+    Draw a number of samples using the given step method.
+    Multiple step methods supported via compound step method
     returns the amount of time taken.
-        
+
     Parameters
     ----------
-        
+
     model : Model (optional if in `with` context)
-    draws : int  
+    draws : int
         The number of samples to draw
     step : function
         A step function
-    start : dict 
+    start : dict
         Starting point in parameter space (Defaults to trace.point(-1))
     trace : NpTrace
         A trace of past values (defaults to None)
     track : list of vars
         The variables to follow
-        
+
     Examples
     --------
-        
+
     >>> an example
-        
+
     """
     model = modelcontext(model)
     draws = int(draws)
-    if start is None: 
+    if start is None:
         start = trace[-1]
     point = Point(start, model = model)
 
     if not hasattr(trace, 'record'):
-        if trace is None: 
+        if trace is None:
             trace = model.vars
         trace = NpTrace(list(trace))
 
@@ -57,16 +58,63 @@ def sample(draws, step, start = None, trace = None, track_progress = True, model
     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)
-  
-def psample(draws, step, start, mtrace = None, track = None, model = None, threads = None):
+
+def psample(draws, step, start, mtrace=None, track=None, model=None, threads=None):
     """draw a number of samples using the given step method. Multiple step methods supported via compound step method
     returns the amount of time taken"""
 
@@ -78,7 +126,7 @@ def psample(draws, step, start, mtrace = None, track = None, model = None, threa
     if isinstance(start, dict) :
         start = threads * [start]
 
-    if track is None: 
+    if track is None:
         track = model.vars
 
     if not mtrace:
@@ -87,7 +135,7 @@ def psample(draws, step, start, mtrace = None, track = None, model = None, threa
     p = mp.Pool(threads)
 
     argset = zip([draws]*threads, [step]*threads, start, mtrace.traces, [False]*threads, [model] *threads)
-    
+
     traces = p.map(argsample, argset)
-        
+
     return MultiTrace(traces)

pymc/step_methods/arraystep.py

@@ -7,10 +7,11 @@
 
 # TODO Add docstrings to ArrayStep
 class ArrayStep(object):
-    def __init__(self, vars, fs, allvars = False):
+    def __init__(self, vars, fs, allvars = False, tune=False):
         self.ordering = ArrayOrdering(vars)
         self.fs = fs
         self.allvars = allvars
+        self.tune = tune
 
     def step(self, point):
         bij = DictToArrayBijection(self.ordering, point)
@@ -28,10 +29,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
+        return q, True
     else:
         # Reject proposed value
-        return q0
+        return q0, False
 
 
 class SamplerHist(object):

pymc/step_methods/compound.py

@@ -2,12 +2,13 @@
 Created on Mar 7, 2011
 
 @author: johnsalvatier
-''' 
-from ..quickclass import * 
+'''
+from ..quickclass import *
 
 @quickclass(object)
 def CompoundStep(methods):
-    methods = list(methods) 
+    methods = list(methods)
+    tune = False
     def step(point):
         for method in methods:
             point = method.step(point)

pymc/step_methods/hmc.py

@@ -2,16 +2,16 @@
 Created on Mar 7, 2011
 
 @author: johnsalvatier
-''' 
+'''
 from numpy import floor
 from quadpotential import *
 from arraystep import *
-from ..core import * 
+from ..core import *
 import numpy as np
 
 __all__ = ['HamiltonianMC']
 
-#TODO: 
+#TODO:
 #add constraint handling via page 37 of Radford's http://www.cs.utoronto.ca/~radford/ham-mcmc.abstract.html
 
 def unif(step_size, elow = .85, ehigh = 1.15):
@@ -30,20 +30,20 @@ def __init__(self, vars, C, step_scale = .25, path_length = 2., is_cov = False,
             step_scale : float, default=.25
                 Size of steps to take, automatically scaled down by 1/n**(1/4) (defaults to .25)
             path_length : float, default=2
-                total length to travel 
+                total length to travel
             is_cov : bool, default=False
                 Treat C as a covariance matrix/vector if True, else treat it as a precision matrix/vector
-            step_rand : function float -> float, default=unif 
-                A function which takes the step size and returns an new one used to randomize the step size at each iteration. 
-            state 
+            step_rand : function float -> float, default=unif
+                A function which takes the step size and returns an new one used to randomize the step size at each iteration.
+            state
                 State object
             model : Model
         """
         model = modelcontext(model)
         n = C.shape[0]
-        
+
         self.step_size = step_scale / n**(1/4.)
-        
+
         self.potential = quad_potential(C, is_cov)
         self.path_length = path_length
         self.step_rand = step_rand
@@ -52,37 +52,37 @@ def __init__(self, vars, C, step_scale = .25, path_length = 2., is_cov = False,
             state = SamplerHist()
         self.state = state
 
-        ArrayStep.__init__(self, 
+        ArrayStep.__init__(self,
                 vars, [model.logpc, model.dlogpc(vars)]
                 )
 
     def astep(self, q0, logp, dlogp):
-        
-            
+
+
         #randomize step size
-        e = self.step_rand(self.step_size) 
+        e = self.step_rand(self.step_size)
         nstep = int(floor(self.path_length / self.step_size))
-        
-        q = q0 
+
+        q = q0
         p = p0 = self.potential.random()
-        
+
         #use the leapfrog method
         p = p - (e/2) * -dlogp(q) # half momentum update
-        
-        for i in range(nstep): 
+
+        for i in range(nstep):
             #alternate full variable and momentum updates
             q = q + e * self.potential.velocity(p)
             if i != nstep - 1:
                 p = p - e * -dlogp(q)
-             
+
         p = p - (e/2) * -dlogp(q)  # do a half step momentum update to finish off
-        
-        p = -p 
-            
+
+        p = -p
+
         # - H(q*, p*) + H(q, p) = -H(q, p) + H(q0, p0) = -(- logp(q) + K(p)) + (-logp(q0) + K(p0))
         mr = (-logp(q0)) + self.potential.energy(p0) - ((-logp(q))  + self.potential.energy(p))
 
         self.state.metrops.append(mr)
-        
-        return metrop_select(mr, q, q0)
-        
+
+        return metrop_select(mr, q, q0)[0]
+

pymc/step_methods/metropolis.py

@@ -51,13 +51,17 @@ def __init__(self, vars, S, proposal_dist=quad_potential, scaling=1.,
             # quadpotential does not require n
             self.proposal_dist = proposal_dist(S)
         self.scaling = scaling
-        self.tune = tune
-        super(Metropolis,self).__init__(vars, [model.logpc])
+        self.accepted = 0
+        super(Metropolis,self).__init__(vars, [model.logpc], tune=tune)
 
     def astep(self, q0, logp):
 
         delta = self.proposal_dist() * self.scaling
 
         q = q0 + delta
 
-        return metrop_select(logp(q) - logp(q0), q, q0)
+        q_new, accepted = metrop_select(logp(q) - logp(q0), q, q0)
+
+        self.accepted += int(accepted)
+
+        return q_new