many small fixes related to simple defaults

Author: jsalvatier

examples/logistic.py

@@ -36,6 +36,6 @@ def tinvlogit(x):
     start = find_MAP()
     h = np.diag(approx_hess(start)) #find a good orientation using the hessian at the MAP
 
-    step = HamiltonianMC(model, model.vars, h) 
+    step = HamiltonianMC(model.vars, h) 
 
     trace = sample(3e2, step, start)

examples/stochastic_volatility.ipynb

@@ -139,7 +139,7 @@
      "collapsed": false,
      "input": [
       "with model:\n",
-      "    start = find_MAP(model, vars = [s], fmin = optimize.fmin_l_bfgs_b)"
+      "    start = find_MAP(vars = [s], fmin = optimize.fmin_l_bfgs_b)"
      ],
      "language": "python",
      "metadata": {},

examples/stochastic_volatility.py

@@ -14,15 +14,15 @@
 # <markdowncell>
 
 # Asset prices have time-varying volatility (variance of day over day `returns`). In some periods, returns are highly vaiable, and in others very stable. Stochastic volatility models model this with a latent volatility variable, modeled as a stochastic process. The following model is similar to the one described in the No-U-Turn Sampler paper, Hoffman (2011) p21.
-#
+# 
 # $$ \sigma \sim Exponential(50) $$
-#
+# 
 # $$ \nu \sim Exponential(.1) $$
-#
+# 
 # $$ s_i \sim Normal(s_{i-1}, \sigma^{-2}) $$
-#
+# 
 # $$ log(\frac{y_i}{y_{i-1}}) \sim t(\nu, 0, exp(-2 s_i)) $$
-#
+# 
 # Here, $y$ is the daily return series and $s$ is the latent volatility process.
 
 # <markdowncell>
@@ -58,16 +58,16 @@
 
 # Fit Model
 # ------------
-# To get a decent scale for the hamiltonaian sampler, we find the hessian at a point. However, the 2nd derivatives for the degrees of freedom are negative and thus not very informative, so we make an educated guess. The interactions between `log_sigma`/`nu` and `s` are also not very useful, so we set them to zero.
-#
+# To get a decent scale for the hamiltonaian sampler, we find the hessian at a point. However, the 2nd derivatives for the degrees of freedom are negative and thus not very informative, so we make an educated guess. The interactions between `log_sigma`/`nu` and `s` are also not very useful, so we set them to zero. 
+# 
 # The hessian matrix is also very sparse, so we make it a sparse matrix for faster sampling.
 
 # <codecell>
 
 H = model.d2logpc()
 
-def hessian(point, nusd):
-    h = H(point)
+def hessian(point, nusd): 
+    h = H(Point(point))
     h[1,1] = nusd**-2
     h[:2,2:] = h[2:,:2] = 0
 
@@ -80,7 +80,7 @@ def hessian(point, nusd):
 # <codecell>
 
 with model:
-    start = find_MAP(model, vars = [s], fmin = optimize.fmin_l_bfgs_b)
+    start = find_MAP(vars = [s], fmin = optimize.fmin_l_bfgs_b)
 
 # <markdowncell>
 
@@ -90,7 +90,7 @@ def hessian(point, nusd):
 
 with model: 
     step = HamiltonianMC(model.vars, hessian(start, 6))
-    trace = sample(200, step, start) 
+    trace = sample(200, step, start, trace = model.vars + [sigma]) 
 
     start2 = trace.point(-1)
     step = HamiltonianMC(model.vars, hessian(start2, 6), path_length = 4.)
@@ -99,15 +99,19 @@ def hessian(point, nusd):
 # <codecell>
 
 #figsize(12,6)
-plt.title(str(s))
-plt.plot(trace[s][::10].T,'b', alpha = .01);
+title(str(s))
+plot(trace[s][::10].T,'b', alpha = .01);
 
 #figsize(12,6)
 traceplot(trace, model.vars[:-1]);
 
+# <codecell>
+
+trace.samples.keys()
+
 # <markdowncell>
 
 # References
 # -------------
-#     1. Hoffman & Gelman. (2011). The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. http://arxiv.org/abs/1111.4246
+#     1. Hoffman & Gelman. (2011). The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. http://arxiv.org/abs/1111.4246 
 

pymc/model.py

@@ -40,7 +40,7 @@ def get_context(cls):
 
 def modelcontext(model): 
     if model is None:
-       return Model.get_context() 
+        return Model.get_context() 
     return model
 
 class Model(Context):
@@ -87,7 +87,7 @@ def d2logpc(model, vars = None):
     @property
     def test_point(self):
         """Test point used to check that the model doesn't generate errors"""
-        return Point(self, ((var, var.tag.test_value) for var in self.vars))
+        return Point(((var, var.tag.test_value) for var in self.vars), model = self)
 
     @property
     def cont_vars(model):

pymc/sample.py

@@ -39,7 +39,7 @@ def sample(draws, step, start = None, trace = None, track_progress = True, model
     draws = int(draws)
     if start is None: 
         start = trace[-1]
-    point = Point(model, start)
+    point = Point(start, model = model)
 
     if not hasattr(trace, 'record'):
         if trace is None: 
@@ -66,7 +66,7 @@ def argsample(args):
     """ defined at top level so it can be pickled"""
     return sample(*args)
 
-def psample(draws, step, start, mtrace = None, threads = None, track = None, model = 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"""
 
@@ -86,7 +86,7 @@ def psample(draws, step, start, mtrace = None, threads = None, track = None, mod
 
     p = mp.Pool(threads)
 
-    argset = zip([model] *threads, [draws]*threads, [step]*threads, start, mtrace.traces)
+    argset = zip([draws]*threads, [step]*threads, start, mtrace.traces, [False]*threads, [model] *threads)
 
     traces = p.map(argsample, argset)
 

pymc/tuning/scaling.py

@@ -24,7 +24,7 @@ def approx_hess(point, vars=None, model = None):
     if vars is None :
         vars = model.cont_vars
 
-    point = Point(model, point)
+    point = Point(point, model = model)
 
     bij = DictToArrayBijection(ArrayOrdering(vars), point)
     dlogp = bij.mapf(model.dlogpc(vars))
@@ -53,7 +53,7 @@ def find_hessian(point, vars = None, model = None):
     """
     model = modelcontext(model)
     H = model.d2logpc(vars)
-    return H(Point(model, point))
+    return H(Point(point, model = model))
 
 def trace_cov(trace, vars = None):
     """

pymc/tuning/starting.py

@@ -11,7 +11,7 @@
 __all__ = ['find_MAP', 'scipyminimize']
 
 
-def find_MAP(start = None, vars=None, fmin = optimize.fmin_bfgs, return_raw = False, disp = False, model = None, model = None, *args, **kwargs):
+def find_MAP(start = None, vars=None, fmin = optimize.fmin_bfgs, return_raw = False, disp = False, model = None, *args, **kwargs):
     """
     Sets state to the local maximum a posteriori point given a model.
     Current default of fmin_Hessian does not deal well with optimizing close 
@@ -37,7 +37,7 @@ def find_MAP(start = None, vars=None, fmin = optimize.fmin_bfgs, return_raw = Fa
     if vars is None: 
         vars = model.cont_vars
 
-    start = Point(model, start)
+    start = Point(start, model = model)
     bij = DictToArrayBijection(ArrayOrdering(vars), start)
 
     logp = bij.mapf(model.logpc)

tests/test_distributions.py

@@ -115,7 +115,7 @@ def check_int_to_1(model, value, domains):
 
     for a in its.product(*domains):
         a = a + (value.tag.test_value,)
-        pt = Point(dict( (str(var), val) for var,val in zip(model.vars, a)))
+        pt = Point(dict( (str(var), val) for var,val in zip(model.vars, a)), model = model)
 
         bij = DictToVarBijection(value,() ,  pt)
 
@@ -145,6 +145,6 @@ def check_dlogp(model, value, domains):
     ndlogp = Gradient(logp)
 
     for a in its.product(*domains):
-        pt = Point(dict( (str(var), val) for var,val in zip(model.vars, a)))
+        pt = Point(dict( (str(var), val) for var,val in zip(model.vars, a)), model = model)
 
         pt = bij.map(pt)

tests/test_sampling.py

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