Merge pull request #214 from pymc-devs/metropolis

Add self-tuning to metropolis step method.

Author: Chris 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/__init__.py

@@ -1,12 +1,14 @@
-from core import * 
+from core import *
 from distributions import *
-from math import * 
+from math import *
 
 from trace import *
-from sample import * 
+from sample import *
 from step_methods import *
 from tuning import *
 
-from debug import * 
+from debug import *
 
 from plots import *
+
+from .tests import test

pymc/model.py

@@ -1,23 +1,18 @@
-'''
-Created on Mar 7, 2011
-
-@author: johnsalvatier
-'''
 from point import *
 from types import *
 
 from theano import theano, tensor as t, function
 from theano.gof.graph import inputs
 
-import numpy as np 
+import numpy as np
 from functools import wraps
 
-__all__ = ['Model', 'compilef', 'gradient', 'hessian', 'modelcontext', 'Point'] 
+__all__ = ['Model', 'compilef', 'gradient', 'hessian', 'modelcontext', 'Point']
 
 
 
-class Context(object): 
-    def __enter__(self): 
+class Context(object):
+    def __enter__(self):
         type(self).get_contexts().append(self)
         return self
 
@@ -28,7 +23,7 @@ def __exit__(self, typ, value, traceback):
     def get_contexts(cls):
         if not hasattr(cls, "contexts"):
             cls.contexts = []
-            
+
         return cls.contexts
 
     @classmethod
@@ -38,45 +33,45 @@ def get_context(cls):
         except IndexError:
             raise TypeError("No context on context stack")
 
-def modelcontext(model): 
+def modelcontext(model):
     if model is None:
-        return Model.get_context() 
+        return Model.get_context()
     return model
 
 class Model(Context):
     """
-    Base class for encapsulation of the variables and 
+    Base class for encapsulation of the variables and
     likelihood factors of a model.
     """
 
     def __init__(self):
         self.vars = []
         self.factors = []
-    
+
     @property
     def logp(model):
         """
         log-probability of the model
-            
+
         Parameters
         ----------
-            
-        model : Model  
+
+        model : Model
 
         Returns
         -------
 
         logp : Theano scalar
-            
+
         """
         return t.add(*map(t.sum, model.factors))
 
     @property
-    def logpc(model): 
+    def logpc(model):
         """Compiled log probability density function"""
         return compilef(model.logp)
 
-    def dlogpc(model, vars = None): 
+    def dlogpc(model, vars = None):
         """Compiled log probability density gradient function"""
         return compilef(gradient(model.logp, vars))
 
@@ -92,7 +87,7 @@ def test_point(self):
     @property
     def cont_vars(model):
         """All the continuous variables in the model"""
-        return typefilter(model.vars, continuous_types) 
+        return typefilter(model.vars, continuous_types)
 
     """
     these functions add random variables
@@ -108,8 +103,8 @@ def Var(model, name, dist):
         model.factors.append(dist.logp(var))
         return var
 
-    def TransformedVar(model, name, dist, trans): 
-        tvar = model.Var(trans.name + '_' + name, trans.apply(dist)) 
+    def TransformedVar(model, name, dist, trans):
+        tvar = model.Var(trans.name + '_' + name, trans.apply(dist))
 
         return named(trans.backward(tvar),name), tvar
 
@@ -118,26 +113,26 @@ def AddPotential(model, potential):
 
 
 def Point(*args,**kwargs):
-    """ 
-    Build a point. Uses same args as dict() does. 
-    Filters out variables not in the model. All keys are strings.  
+    """
+    Build a point. Uses same args as dict() does.
+    Filters out variables not in the model. All keys are strings.
 
     Parameters
     ----------
-        *args, **kwargs 
+        *args, **kwargs
             arguments to build a dict
     """
     if 'model' in kwargs :
         model = kwargs['model']
         del kwargs['model']
-    else: 
+    else:
         model = Model.get_context()
 
     d = dict(*args, **kwargs)
     varnames = map(str, model.vars)
-    return dict((str(k),np.array(v)) 
-            for (k,v) in d.iteritems() 
-            if str(k) in varnames) 
+    return dict((str(k),np.array(v))
+            for (k,v) in d.iteritems()
+            if str(k) in varnames)
 
 
 def compilef(outs, mode = None):
@@ -148,42 +143,42 @@ def compilef(outs, mode = None):
     ----------
     outs : Theano variable or iterable of Theano variables
     mode : Theano compilation mode
-    
+
     Returns
     -------
     Compiled Theano function
     """
     return PointFunc(
-                function(inputvars(outs), outs, 
-                         allow_input_downcast = True, 
+                function(inputvars(outs), outs,
+                         allow_input_downcast = True,
                          on_unused_input = 'ignore',
                          mode = mode)
            )
 
 def named(var, name):
-    var.name = name 
+    var.name = name
     return var
 
 def as_iterargs(data):
-    if isinstance(data, tuple): 
+    if isinstance(data, tuple):
         return data
     if hasattr(data, 'columns'): #data frames
-        return [np.asarray(data[c]) for c in data.columns] 
+        return [np.asarray(data[c]) for c in data.columns]
     else:
         return [data]
 
-def makeiter(a): 
+def makeiter(a):
     if isinstance(a, (tuple, list)):
         return a
     else :
         return [a]
 
-def inputvars(a): 
+def inputvars(a):
     return [v for v in inputs(makeiter(a)) if isinstance(v, t.TensorVariable)]
 
 """
-Theano derivative functions 
-""" 
+Theano derivative functions
+"""
 
 def cont_inputs(f):
     return typefilter(inputvars(f), continuous_types)
@@ -193,7 +188,7 @@ def gradient1(f, v):
     return t.flatten(t.grad(f, v, disconnected_inputs='warn'))
 
 def gradient(f, vars = None):
-    if not vars: 
+    if not vars:
         vars = cont_inputs(f)
 
     return t.concatenate([gradient1(f, v) for v in vars], axis = 0)
@@ -202,14 +197,14 @@ def jacobian1(f, v):
     """jacobian of f wrt v"""
     f = t.flatten(f)
     idx = t.arange(f.shape[0])
-    
-    def grad_i(i): 
+
+    def grad_i(i):
         return gradient1(f[i], v)
 
     return theano.map(grad_i, idx)[0]
 
 def jacobian(f, vars = None):
-    if not vars: 
+    if not vars:
         vars = cont_inputs(f)
 
     return t.concatenate([jacobian1(f, v) for v in vars], axis = 1)
@@ -218,6 +213,6 @@ def hessian(f, vars = None):
     return -jacobian(gradient(f, vars), vars)
 
 
-#theano stuff 
+#theano stuff
 theano.config.warn.sum_div_dimshuffle_bug = False
 theano.config.compute_test_value = 'raise'

pymc/sample.py

@@ -8,41 +8,41 @@
 
 __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, 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))
 
@@ -65,8 +65,8 @@ def sample(draws, step, start = None, trace = None, track_progress = True, model
 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 +78,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 +87,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)