First stab at revising Metropolis. Can now specify alternative proposals. Still need to add tuning.

Author: Christopher Fonnesbeck

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/step_methods/arraystep.py

@@ -1,32 +1,38 @@
-from ..core import * 
-import numpy as np 
+from ..core import *
+import numpy as np
 from numpy.random import uniform
 from numpy import log , isfinite
 
 __all__ = ['ArrayStep', 'metrop_select', 'SamplerHist']
 
+# TODO Add docstrings to ArrayStep
 class ArrayStep(object):
     def __init__(self, vars, fs, allvars = False):
         self.ordering = ArrayOrdering(vars)
         self.fs = fs
         self.allvars = allvars
-        
+
     def step(self, point):
         bij = DictToArrayBijection(self.ordering, point)
-        
-        inputs = map(bij.mapf, self.fs) 
+
+        inputs = map(bij.mapf, self.fs)
         if self.allvars:
             inputs += [point]
 
         apoint = self.astep(bij.map(point), *inputs)
         return bij.rmap(apoint)
 
 def metrop_select(mr, q, q0):
+    # Perform rejection/acceptance step for Metropolis class samplers
+
+    # Compare acceptance ratio to uniform random number
     if isfinite(mr) and log(uniform()) < mr:
+        # Accept proposed value
         return q
-    else: 
+    else:
+        # Reject proposed value
         return q0
-    
+
 
 class SamplerHist(object):
     def __init__(self):

pymc/step_methods/metropolis.py

@@ -1,31 +1,63 @@
-'''
-Created on Mar 7, 2011
-
-@author: johnsalvatier
-'''
 from numpy.linalg import cholesky
 
 from ..core import *
 from quadpotential import quad_potential
 
 from arraystep import *
+from numpy.random import normal, standard_cauchy, standard_exponential
 
 __all__ = ['Metropolis']
 
-# TODO Implement tuning for Metropolis step
+# Available proposal distributions for Metropolis
+def normal_proposal(s, n):
+    return lambda: normal(scale=s, size=n)
+
+def cauchy_proposal(s, n):
+    return lambda: standard_cauchy(size=n) * s
+
+def laplace_proposal(s, n):
+    return lambda: (standard_exponential(size=n) - standard_exponential(size=n)) * s
+
+
 class Metropolis(ArrayStep):
-    def __init__(self, vars, C, scaling=.25, is_cov = False, model = None):
+    """
+    Metropolis-Hastings sampling step
+
+    Parameters
+    ----------
+    vars : list
+        List of variables for sampler
+    S : standard deviation or covariance matrix
+        Some measure of variance to parameterize proposal distribution
+    proposal_dist : function
+        Function that returns zero-mean deviates when parameterized with
+        S (and n). Defaults to quad_potential.
+    scaling : scalar or array
+        Initial scale factor for proposal. Defaults to 1.
+    tune : bool
+        Flag for tuning. Defaults to True.
+    model : PyMC Model
+        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):
+
         model = modelcontext(model)
 
-        self.potential = quad_potential(C, not is_cov)
+        try:
+            self.proposal_dist = proposal_dist(S, n=len(vars))
+        except TypeError:
+            # quadpotential does not require n
+            self.proposal_dist = proposal_dist(S)
         self.scaling = scaling
+        self.tune = tune
         super(Metropolis,self).__init__(vars, [model.logpc])
-        
+
     def astep(self, q0, logp):
 
-        delta = self.potential.random() * self.scaling
-        
-        q = q0 + delta  
-        
-        return metrop_select(logp(q) - logp(q0),
-                                    q, q0)
+        delta = self.proposal_dist() * self.scaling
+
+        q = q0 + delta
+
+        return metrop_select(logp(q) - logp(q0), q, q0)