added deterministic traces

Author: jsalvatier

examples/stochastic_volatility.ipynb

@@ -12,7 +12,7 @@
 import numpy as np 
 from functools import wraps
 
-__all__ = ['Model', 'compilef', 'gradient', 'hessian', 'withmodel'] 
+__all__ = ['Model', 'compilef', 'gradient', 'hessian', 'withmodel', 'Point'] 
 
 
 
@@ -109,13 +109,33 @@ def Var(model, name, dist):
 
     def TransformedVar(model, name, dist, trans): 
         tvar = model.Var(trans.name + '_' + name, trans.apply(dist)) 
-        return trans.backward(tvar), tvar
+
+        return named(trans.backward(tvar),name), tvar
 
     def AddPotential(model, potential):
         model.factors.append(potential)
 
 withmodel = withcontext(Model, 'model')
 
+@withmodel
+def Point(model, *args,**kwargs):
+    """ 
+    Build a point. Uses same args as dict() does. 
+    Filters out variables not in the model. All keys are strings.  
+
+    Parameters
+    ----------
+        model : Model (in context) 
+        *args, **kwargs 
+            arguments to build a dict
+    """
+
+    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) 
+
 
 def compilef(outs, mode = None):
     return PointFunc(
@@ -125,6 +145,9 @@ def compilef(outs, mode = None):
                          mode = mode)
            )
 
+def named(var, name):
+    var.name = name 
+    return var
 
 def as_iterargs(data):
     if isinstance(data, tuple): 

pymc/model.py

@@ -1,9 +1,5 @@
 import numpy as np
 
-def Point(d) : 
-    d = dict(d)
-    return dict((str(k),np.array(v)) for (k,v) in d.iteritems()) 
- 
 class PointFunc(object): 
     def __init__(self, f):
         self.f = f

pymc/point.py

@@ -4,11 +4,12 @@
 from time import time
 from core import *
 import step_methods
+from progressbar import ProgressBar
 
 __all__ = ['sample', 'psample']
 
 @withmodel
-def sample(model, draws, step, start = None, trace = None, vars = None): 
+def sample(model, draws, step, start = None, trace = None, progress_bar = True): 
     """
     Draw a number of samples using the given step method. 
     Multiple step methods supported via compound step method 
@@ -25,33 +26,40 @@ def sample(model, draws, step, start = None, trace = None, vars = None):
         Starting point in parameter space (Defaults to trace.point(-1))
     trace : NpTrace
         A trace of past values (defaults to None)
-    state : 
-        The current state of the sampler (defaults to None)
+    track : list of vars
+        The variables to follow
         
     Examples
     --------
         
     >>> an example
         
     """
+    draws = int(draws)
     if start is None: 
         start = trace[-1]
     point = Point(start)
 
-    if vars is None: 
-        vars = model.vars
+    if not hasattr(trace, 'record'):
+        if trace is None: 
+            trace = model.vars
+        trace = NpTrace(list(trace))
 
-    if trace is None: 
-        trace = NpTrace(vars)
 
     try:
         step = step_methods.CompoundStep(step)
     except TypeError:
         pass
 
-    for _ in xrange(int(draws)):
+
+    
+    progress = ProgressBar(draws)
+
+    for i in xrange(draws):
         point = step.step(point)
         trace = trace.record(point)
+        if progress_bar:
+            progress.animate(i)
 
     return trace
 
@@ -60,7 +68,7 @@ def argsample(args):
     return sample(*args)
 
 @withmodel
-def psample(model, draws, step, start, mtrace = None, threads = None, vars = None):
+def psample(model, draws, step, start, mtrace = None, threads = None, track = 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"""
 
@@ -70,11 +78,11 @@ def psample(model, draws, step, start, mtrace = None, threads = None, vars = Non
     if isinstance(start, dict) :
         start = threads * [start]
 
-    if vars is None: 
-        vars = model.vars
+    if track is None: 
+        track = model.vars
 
     if not mtrace:
-        mtrace = MultiTrace(threads, vars)
+        mtrace = MultiTrace(threads, track)
 
     p = mp.Pool(threads)
 

pymc/sample.py

@@ -36,7 +36,7 @@ def __getitem__(self, key):
         return self.samples[str(key)].value
 
     def point(self, index):
-        return Point((k, v.value[index]) for (k,v) in self.samples.iteritems())
+        return dict((k, v.value[index]) for (k,v) in self.samples.iteritems())
 
 class ListArray(object):
     def __init__(self):

pymc/trace.py

@@ -1,2 +1,2 @@
-from starting import find_MAP
-from scaling import approx_hess
+from starting import *
+from scaling import *

pymc/tuning/__init__.py

@@ -7,7 +7,7 @@
 import numpy as np 
 from ..core import *
 
-__all__ = ['approx_hess']
+__all__ = ['approx_hess', 'find_hessian', 'trace_cov']
 
 @withmodel
 def approx_hess(model, start, vars=None):
@@ -24,7 +24,7 @@ def approx_hess(model, start, vars=None):
     if vars is None :
         vars = model.cont_vars
 
-    start = Point(start)
+    start = Point(model, start)
 
     bij = DictToArrayBijection(ArrayOrdering(vars), start)
     dlogp = bij.mapf(model.dlogpc(vars))
@@ -40,6 +40,10 @@ def grad_logp(point):
     '''
     return -nd.Jacobian(grad_logp)(bij.map(start))
 
+@withmodel
+def find_hessian(model, point, vars = None): 
+    H = model.d2logpc(vars)
+    return H(Point(model, point))
 
 def trace_cov(trace, vars = None):
     """

pymc/tuning/scaling.py

@@ -37,7 +37,7 @@ def find_MAP(model, start = None, vars=None, fmin = optimize.fmin_bfgs, return_r
     if vars is None: 
         vars = model.cont_vars
 
-    start = Point(start)
+    start = Point(model, start)
     bij = DictToArrayBijection(ArrayOrdering(vars), start)
 
     logp = bij.mapf(model.logpc)