Merge pull request #221 from pymc-devs/slice

Slice sampler step function

Author: jsalvatier

pymc/sample.py

@@ -17,7 +17,6 @@ def sample(draws, step, start=None, trace=None, track_progress=True, model=None)
     Parameters
     ----------
 
-    model : Model (optional if in `with` context)
     draws : int
         The number of samples to draw
     step : function
@@ -26,8 +25,9 @@ def sample(draws, step, start=None, trace=None, track_progress=True, model=None)
         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
+    track_progress : bool
+        Flag for progress bar
+    model : Model (optional if in `with` context)
 
     Examples
     --------

pymc/step_methods/__init__.py

@@ -1,4 +1,5 @@
 from compound import *
 from hmc import *
-from metropolis import *    
+from metropolis import *
 from gibbs import *
+from slicer import *

pymc/step_methods/slicer.py

@@ -0,0 +1,64 @@
+# Modified from original implementation by Dominik Wabersich (2013)
+
+from ..core import *
+from arraystep import *
+from numpy import floor, abs, atleast_1d
+from numpy.random import standard_exponential, random, uniform
+
+__all__ = ['Slice']
+
+class Slice(ArrayStep):
+    """Slice sampler"""
+    def __init__(self, vars, w=1, m=20, n_tune=0,
+            tune=True, tune_interval=100, model=None):
+
+        model = modelcontext(model)
+
+        self.vars = vars
+        self.w = w
+        self.m = m
+        self.n_tune = n_tune
+        self.w_tune = []
+        self.tune = tune
+        self.tune_interval = tune_interval
+        self.model = model
+
+        super(Slice,self).__init__(vars, [model.logpc])
+
+    def astep(self, q0, logp):
+
+        y = logp(q0) - standard_exponential()
+
+        # Stepping out procedure
+        L = q0 - self.w*random()
+        R = L + self.w
+        J = floor(self.m*random())
+        K = (self.m - 1) - J
+
+        while(J>0 and y<logp(L)):
+            L = L - self.w
+            J = J - 1
+
+        while(K>0 and y<logp(R)):
+            R = R + self.w
+            K = K - 1
+
+        # Sample uniformly from slice
+        q_new = atleast_1d(uniform(L, R))
+        y_new = logp(q_new)
+
+        while(y_new<y):
+            # Shrink bounds of uniform
+            if (q_new < q0):
+                L = q_new
+            else:
+                R = q_new
+            q_new = atleast_1d(uniform(L,R))
+            y_new = logp(q_new)
+
+        if (len(self.w_tune) > self.n_tune):
+            # Tune sampler parameters
+            self.w_tune.append(abs(q0 - q_new))
+            self.w = 2 * sum(self.w_tune)/len(self.w_tune)
+
+        return q_new

pymc/tests/models.py

@@ -1,14 +1,7 @@
 from pymc import Model, Normal, Metropolis
-import numpy as np 
+import numpy as np
 import pymc as pm
 
-def simple_init(): 
-    start, model, moments = simple_model()
-
-    step = Metropolis(model.vars, np.diag([1.]), model = model)
-    return model, start, step, moments
-
-
 def simple_model():
     mu = -2.1
     tau = 1.3
@@ -17,6 +10,12 @@ def simple_model():
 
     return model.test_point, model, (mu, tau**-1)
 
+def simple_init():
+    start, model, moments = simple_model()
+
+    step = Metropolis(model.vars, np.diag([1.]), model = model)
+    return model, start, step, moments
+
 def simple_2model():
     mu = -2.1
     tau = 1.3
@@ -34,7 +33,7 @@ def mv_simple():
         [.05 , .1,  0],
         [1. ,-0.05,5.5]])
 
-    tau = np.dot(p,p.T) 
+    tau = np.dot(p,p.T)
 
 
     with pm.Model() as model:

pymc/tests/test_sampling.py

@@ -17,9 +17,9 @@ def test_sample():
         test_samplers.append(psample)
 
     with model:
-        for     samplr  in test_samplers: 
-            for n       in [0, 10, 1000]:
-                yield samplr, n, step, start 
+        for samplr in test_samplers:
+            for n in [0, 10, 1000]:
+                yield samplr, n, step, start
 
 
 

pymc/tests/test_step.py

@@ -1,5 +1,5 @@
-from checks import *
-from models import simple_model, mv_simple
+from .checks import *
+from .models import simple_model, mv_simple
 from theano.tensor import constant
 from scipy.stats.mstats import moment
 
@@ -11,26 +11,23 @@ def check_stat(name, trace, var, stat, value, bound):
 def test_step_continuous():
     start, model, (mu, C) = mv_simple()
 
-    hmc = pm.HamiltonianMC(model.vars,  C, is_cov = True, model = model)
-    mh = pm.Metropolis(model.vars , C, is_cov = True, scaling = 2, model = model)
+    hmc = pm.HamiltonianMC(model.vars,  C, is_cov=True, model=model)
+    mh = pm.Metropolis(model.vars, model=model)
+    slicer = pm.Slice(model.vars, model=model)
     compound = pm.CompoundStep([hmc, mh])
 
-    steps = [mh, hmc, compound]
+    steps = [mh, hmc, compound, slicer]
+
 
     unc = np.diag(C)**.5
     check = [('x', np.mean, mu, unc/10.),
              ('x', np.std , unc, unc/10.)]
 
     for st in steps:
         np.random.seed(1)
-        h = sample(8000, st, start, track_progress=False, model = model)
+        h = sample(8000, st, start, model=model)
         for (var, stat, val, bound) in check:
             np.random.seed(1)
-            h = sample(8000, st, start, track_progess=False, model = model)
-
-            yield check_stat,repr(st), h, var, stat, val, bound
-
-
-
-
+            h = sample(8000, st, start, model=model)
 
+            yield check_stat,repr(st), h, var, stat, val, bound