Merge branch 'master' into sample_SMC

Author: aloctavodia

RELEASE-NOTES.md

@@ -13,6 +13,11 @@
 - Improve error message `NaN occurred in optimization.` during ADVI
 - Save and load traces without `pickle` using `pm.save_trace` and `pm.load_trace`
 - Add `Kumaraswamy` distribution
+- Rewrite parallel sampling of multiple chains on py3. This resolves
+  long standing issues when tranferring large traces to the main process,
+  avoids pickleing issues on UNIX, and allows us to show a progress bar
+  for all chains. If parallel sampling is interrupted, we now return
+  partial results.
 
 ### Fixes
 

benchmarks/benchmarks/benchmarks.py

@@ -1,4 +1,5 @@
 import time
+import timeit
 
 import numpy as np
 import pandas as pd
@@ -63,6 +64,7 @@ class OverheadSuite(object):
     samplers
     """
     params = [pm.NUTS, pm.HamiltonianMC, pm.Metropolis, pm.Slice]
+    timer = timeit.default_timer
 
     def setup(self, step):
         self.n_steps = 10000
@@ -78,6 +80,7 @@ def time_overhead_sample(self, step):
 class ExampleSuite(object):
     """Implements examples to keep up with benchmarking them."""
     timeout = 360.0  # give it a few minutes
+    timer = timeit.default_timer
 
     def time_drug_evaluation(self):
         drug = np.array([101, 100, 102, 104, 102, 97, 105, 105, 98, 101,

docs/source/notebooks/dependent_density_regression.ipynb

@@ -99,8 +99,9 @@ def record(self, point):
         self._fh.write(','.join(columns) + '\n')
 
     def close(self):
-        self._fh.close()
-        self._fh = None  # Avoid serialization issue.
+        if self._fh is not None:
+            self._fh.close()
+            self._fh = None  # Avoid serialization issue.
 
     # Selection methods
 

pymc3/backends/text.py

@@ -158,18 +158,33 @@ def random_choice(*args, **kwargs):
                 return np.random.choice(k, p=w, *args, **kwargs)
 
         w = draw_values([self.w], point=point)[0]
-
+        comp_tmp = self._comp_samples(point=point, size=None)
+        if self.shape.size == 0:
+            distshape = np.asarray(np.broadcast(w, comp_tmp).shape)[..., :-1]
+        else:
+            distshape = self.shape
         w_samples = generate_samples(random_choice,
                                      w=w,
                                      broadcast_shape=w.shape[:-1] or (1,),
-                                     dist_shape=self.shape,
+                                     dist_shape=distshape,
                                      size=size).squeeze()
-        comp_samples = self._comp_samples(point=point, size=size)
-
-        if comp_samples.ndim > 1:
-            return np.squeeze(comp_samples[np.arange(w_samples.size), w_samples])
+        if (size is None) or (distshape.size == 0):
+            comp_samples = self._comp_samples(point=point, size=size)
+            if comp_samples.ndim > 1:
+                samples = np.squeeze(comp_samples[np.arange(w_samples.size), ..., w_samples])
+            else:
+                samples = np.squeeze(comp_samples[w_samples])
         else:
-            return np.squeeze(comp_samples[w_samples])
+            samples = np.zeros((size,)+tuple(distshape))
+            for i in range(size):
+                w_tmp = w_samples[i, :]
+                comp_tmp = self._comp_samples(point=point, size=None)
+                if comp_tmp.ndim > 1:
+                    samples[i, :] = np.squeeze(comp_tmp[np.arange(w_tmp.size), ..., w_tmp])
+                else:
+                    samples[i, :] = np.squeeze(comp_tmp[w_tmp])
+
+        return samples
 
 
 class NormalMixture(Mixture):
@@ -197,22 +212,22 @@ class NormalMixture(Mixture):
         the component standard deviations
     tau : array of floats
         the component precisions
+    comp_shape : shape of the Normal component
+        notice that it should be different than the shape
+        of the mixture distribution, with one axis being
+        the number of components.
 
     Note: You only have to pass in sd or tau, but not both.
     """
 
-    def __init__(self, w, mu, *args, **kwargs):
+    def __init__(self, w, mu, comp_shape=(), *args, **kwargs):
         _, sd = get_tau_sd(tau=kwargs.pop('tau', None),
                            sd=kwargs.pop('sd', None))
 
-        distshape = np.broadcast(mu, sd).shape
         self.mu = mu = tt.as_tensor_variable(mu)
         self.sd = sd = tt.as_tensor_variable(sd)
 
-        if not distshape:
-            distshape = np.broadcast(mu.tag.test_value, sd.tag.test_value).shape
-
-        super(NormalMixture, self).__init__(w, Normal.dist(mu, sd=sd, shape=distshape),
+        super(NormalMixture, self).__init__(w, Normal.dist(mu, sd=sd, shape=comp_shape),
                                             *args, **kwargs)
 
     def _repr_latex_(self, name=None, dist=None):

pymc3/distributions/mixture.py

@@ -69,10 +69,42 @@ def setter(self, val):
     return gp_wrapper
 
 
-def plot_gp_dist(ax, samples, x, plot_samples=True, palette="Reds"):
-    """ A helper function for plotting 1D GP posteriors from trace """
+def plot_gp_dist(ax, samples, x, plot_samples=True, palette="Reds", fill_alpha=0.8, samples_alpha=0.1, fill_kwargs=None, samples_kwargs=None):
+    """ A helper function for plotting 1D GP posteriors from trace 
+    
+        Parameters
+    ----------
+    ax : axes
+        Matplotlib axes.
+    samples : trace or list of traces
+        Trace(s) or posterior predictive sample from a GP.
+    x : array
+        Grid of X values corresponding to the samples. 
+    plot_samples: bool
+        Plot the GP samples along with posterior (defaults True).
+    palette: str
+        Palette for coloring output (defaults to "Reds").
+    fill_alpha : float
+        Alpha value for the posterior interval fill (defaults to 0.8).
+    samples_alpha : float
+        Alpha value for the sample lines (defaults to 0.1).
+    fill_kwargs : dict
+        Additional arguments for posterior interval fill (fill_between).
+    samples_kwargs : dict
+        Additional keyword arguments for samples plot.
+
+    Returns
+    -------
+
+    ax : Matplotlib axes
+    """
     import matplotlib.pyplot as plt
 
+    if fill_kwargs is None:
+        fill_kwargs = {}
+    if samples_kwargs is None:
+        samples_kwargs = {}
+
     cmap = plt.get_cmap(palette)
     percs = np.linspace(51, 99, 40)
     colors = (percs - np.min(percs)) / (np.max(percs) - np.min(percs))
@@ -82,8 +114,11 @@ def plot_gp_dist(ax, samples, x, plot_samples=True, palette="Reds"):
         upper = np.percentile(samples, p, axis=1)
         lower = np.percentile(samples, 100-p, axis=1)
         color_val = colors[i]
-        ax.fill_between(x, upper, lower, color=cmap(color_val), alpha=0.8)
+        ax.fill_between(x, upper, lower, color=cmap(color_val), alpha=fill_alpha, **fill_kwargs)
     if plot_samples:
         # plot a few samples
         idx = np.random.randint(0, samples.shape[1], 30)
-        ax.plot(x, samples[:,idx], color=cmap(0.9), lw=1, alpha=0.1)
+        ax.plot(x, samples[:,idx], color=cmap(0.9), lw=1, alpha=samples_alpha,
+                **samples_kwargs)
+
+    return ax