fixed examples and added deterministic traces

Author: jsalvatier

examples/ARM12-6.py

@@ -46,19 +46,19 @@
 
     lr = Normal('lr', floor*floor_m + means[group], sd**-2., observed = lradon)
 
-start = {'groupmean' : obs_means.mean(),
+    start = {'groupmean' : obs_means.mean(),
          'groupsd' : obs_means.std(), 
          'sd' : data.groupby('group').lradon.std().mean(), 
          'means' : np.array(obs_means),
          'floor_m' : 0.,
          }
 
-start = find_MAP(model, start, [groupmean,sd, floor_m])
-H = model.d2logpc()
-h = np.diag(H(start))
+    start = find_MAP(start, [groupmean,sd, floor_m])
+    H = model.d2logpc()
+    h = np.diag(H(start))
 
-step = HamiltonianMC(model, model.vars, h)
+    step = HamiltonianMC(model.vars, h)
+
+    trace = sample(3000, step, start)
 
-trace, state, t = sample(3000, step, start)
 
-print  " took: ", t

examples/ARM12-6uranium.py

@@ -48,20 +48,18 @@
 
     lr = Normal('lr', floor*floor_m + means[group] + ufull*u_m, sd**-2., observed = lradon)
 
-start = {'groupmean' : obs_means.mean(),
+    start = {'groupmean' : obs_means.mean(),
          'groupsd' : obs_means.std(), 
          'sd' : data.groupby('group').lradon.std().mean(), 
          'means' : np.array(obs_means),
          'u_m' : np.array([.72]),
          'floor_m' : 0.,
          }
 
-start = find_MAP(model, start, model.vars[:-1])
-H = model.d2logpc()
-h = np.diag(H(start))
+    start = find_MAP(start, model.vars[:-1])
+    H = model.d2logpc()
+    h = np.diag(H(start))
 
-step = HamiltonianMC(model, model.vars, h)
+    step = HamiltonianMC(model.vars, h)
 
-trace, state, t = sample(3000, step, start)
-
-print  " took: ", t
+    trace = sample(3000, step, start)

examples/ARM5-4.py

@@ -21,20 +21,18 @@
 
 Pa = np.array(P)
 
-model = Model()
-with model:
+with Model() as model:
     effects = Normal('effects', mu = 0, tau = 100.**-2, shape = len(P.columns))
     p = sigmoid(dot(Pa, effects))
 
     s = Bernoulli('s', p, observed = np.array(data.switch))
 
 
-#move the chain to the MAP which should be a good starting point
-start = find_MAP(model)
-H = model.d2logpc() #find a good orientation using the hessian at the MAP
-h = H(start)
+    #move the chain to the MAP which should be a good starting point
+    start = find_MAP()
+    H = model.d2logpc() #find a good orientation using the hessian at the MAP
+    h = H(start)
 
-step_method = HamiltonianMC(model, model.vars, h)
+    step = HamiltonianMC(model.vars, h)
 
-trace, state, t = sample(3e3, step_method, start)
-print "took :", t
+    trace = sample(3e3, step, start)

examples/banana.py

@@ -25,26 +25,24 @@
 http://arxiv.org/abs/1011.0057
 """
 N = 200
-model = Model()
-with model:
+with Model() as model:
 
 
     x = Normal('x', 0, 1)
     y = Normal('y', 0, 1)
     N = 200
     d = Normal('d', x + y**2, 1., observed = np.zeros(N))
 
-start = model.test_point
-h = np.ones(2)*np.diag(approx_hess(model, start))[0]
+    start = model.test_point
+    h = np.ones(2)*np.diag(approx_hess(start))[0]
 
 
-step = HamiltonianMC(model, model.vars, h, path_length = 4.)
+    step = HamiltonianMC(model.vars, h, path_length = 4.)
 
-history, state, t = sample(3e3, step, start)
+    trace = sample(3e3, step, start)
 
-print "took :", t
 pl.figure()
-pl.hexbin(history['x'], history['y'])
+pl.hexbin(trace['x'], trace['y'])
 
 
 

examples/dirichlet.py

@@ -14,11 +14,11 @@
                  simplextransform)
 
 
-H = model.d2logpc()
+    H = model.d2logpc()
 
-s = find_MAP(model)
+    s = find_MAP()
 
-step = HamiltonianMC(model, model.vars, H(s))
-trace, _,t = sample(1000, step, s) 
+    step = HamiltonianMC(model.vars, H(s))
+    trace = sample(1000, step, s) 
 
 

examples/hierarchical.py

@@ -49,10 +49,9 @@
 
 
 
-start = find_MAP(model)
-h = approx_hess(model, start) #find a good orientation using the hessian at the MAP
+    start = find_MAP()
+    h = approx_hess(start) #find a good orientation using the hessian at the MAP
 
-step = HamiltonianMC(model, model.vars, h, is_cov = False)
+    step = HamiltonianMC(model.vars, h, is_cov = False)
 
-
-print "took :", sample(3e3, step, start)
+    trace = sample(3e3, step, start)

examples/logistic.py

@@ -30,11 +30,12 @@ def tinvlogit(x):
 
     o = Bernoulli('o', p, observed = outcomes)
 
-#move the chain to the MAP which should be a good starting point
-start = find_MAP(model)
-h = np.diag(approx_hess(model, start)) #find a good orientation using the hessian at the MAP
+    
 
-step = HamiltonianMC(model, model.vars, h) 
+    #move the chain to the MAP which should be a good starting point
+    start = find_MAP()
+    h = np.diag(approx_hess(start)) #find a good orientation using the hessian at the MAP
 
-history, state, t = sample(3e2, step, start)
-print "took :", t
+    step = HamiltonianMC(model, model.vars, h) 
+
+    trace = psample(3e2, step, start)

examples/normal.py

@@ -13,14 +13,12 @@
 with model: 
     x = Normal('x', 0, 1., n)
 
-#start sampling at the MAP
-start = find_MAP(model)
-h = approx_hess(model, start) #find a good orientation using the hessian at the MAP
+    #start sampling at the MAP
+    start = find_MAP()
+    h = approx_hess(start) #find a good orientation using the hessian at the MAP
 
-#step_method = split_hmc_step(model, model.vars, hess, chain, hmc_cov)
-step = HamiltonianMC(model, model.vars, h, is_cov = False)
+    step = HamiltonianMC(model.vars, h, is_cov = False)
 
-ndraw = 3e3
-history, state, t = sample(ndraw, step, start)
+    ndraw = 3e3
+    trace = sample(ndraw, step, start)
 
-print "took :", t

examples/simpletest.py

@@ -18,23 +18,21 @@
 
     d = Normal('data', mu = x, tau = .75**-2, observed = data)
 
-start = model.test_point
+    start = model.test_point
 
-h = approx_hess(model, start)
-step = HamiltonianMC(model, model.vars, h)
+    h = approx_hess(start)
+    step = HamiltonianMC(model.vars, h)
 
-history,state, t = sample(1e3, step, start)
-
-print "took :", t
+    trace = sample(1e3, step, start)
 
 from pylab import * 
 subplot(2,2,1)
-plot(history['x'][:,0,0])
+plot(trace[x][:,0,0])
 subplot(2,2,2)
-hist(history['x'][:,0,0])
+hist(trace[x][:,0,0])
 
 subplot(2,2,3)
-plot(history['z'])
+plot(trace[x])
 subplot(2,2,4)
-hist(history['x'][:,0,0])
+hist(trace[x][:,0,0])
 show()

examples/speedtest_pymc3.py

@@ -46,16 +46,16 @@
 
 model = Model()
 with model:
-# Estimated occupancy
+    # Estimated occupancy
     p = Beta('b', 1,1)
 
-# Latent variable for occupancy
+    # Latent variable for occupancy
     z = Bernoulli('z', p , shape = y.shape)
 
-# Estimated mean count
+    # Estimated mean count
     theta = Uniform('theta', 0, 100)
 
-# Poisson likelihood
+    # Poisson likelihood
     z = ZeroInflatedPoisson('z', theta, z)