Merge branch 'aloctavodia-sample_ppc_ma'

Author: aloctavodia

.travis.yml

@@ -1,23 +1,25 @@
 language: python
+sudo: false
 
 before_install:
   - . ./scripts/install_miniconda.sh
-  - "export DISPLAY=:99.0"
-  - "sh -e /etc/init.d/xvfb start"
-  - if [ ${MATPLOTLIB} = "1.2" ]; then mkdir $HOME/.matplotlib; fi
-  - if [ ${MATPLOTLIB} = "1.2" ]; then cp ${SRCDIR}/tools/matplotlibrc $HOME/.matplotlib/matplotlibrc; fi
+  - sh -e /etc/init.d/xvfb start
+  - export DISPLAY=":99.0"
 
 install:
   - . ./scripts/create_testenv.sh
   - pip install coveralls pylint
 
 env:
-  - PYTHON_VERSION=2.7 TESTCMD="--durations=10 --ignore=pymc3/tests/test_examples.py --cov-append --ignore=pymc3/tests/test_distributions_random.py --ignore=pymc3/tests/test_variational_inference.py --ignore=pymc3/tests/test_shared.py --ignore=pymc3/tests/test_smc.py --ignore=pymc3/tests/test_updates.py"
-  - PYTHON_VERSION=2.7 RUN_PYLINT="true" TESTCMD="--durations=10 --cov-append pymc3/tests/test_distributions_random.py pymc3/tests/test_shared.py pymc3/tests/test_smc.py"
-  - PYTHON_VERSION=2.7 TESTCMD="--durations=10 --cov-append pymc3/tests/test_examples.py pymc3/tests/test_variational_inference.py pymc3/tests/test_updates.py"
-  - PYTHON_VERSION=3.6 TESTCMD="--durations=10 --cov-append --ignore=pymc3/tests/test_examples.py --ignore=pymc3/tests/test_distributions_random.py --ignore=pymc3/tests/test_variational_inference.py --ignore=pymc3/tests/test_shared.py --ignore=pymc3/tests/test_smc.py --ignore=pymc3/tests/test_updates.py"
-  - PYTHON_VERSION=3.6 TESTCMD="--durations=10 --cov-append pymc3/tests/test_distributions_random.py pymc3/tests/test_shared.py pymc3/tests/test_smc.py"
-  - PYTHON_VERSION=3.6 TESTCMD="--durations=10 --cov-append pymc3/tests/test_examples.py pymc3/tests/test_variational_inference.py pymc3/tests/test_updates.py"
+  - PYTHON_VERSION=2.7 FLOATX='float32' TESTCMD="--durations=10 --ignore=pymc3/tests/test_examples.py --cov-append --ignore=pymc3/tests/test_distributions_random.py --ignore=pymc3/tests/test_variational_inference.py --ignore=pymc3/tests/test_shared.py --ignore=pymc3/tests/test_smc.py --ignore=pymc3/tests/test_updates.py"
+  - PYTHON_VERSION=2.7 FLOATX='float32' RUN_PYLINT="true" TESTCMD="--durations=10 --cov-append pymc3/tests/test_distributions_random.py pymc3/tests/test_shared.py pymc3/tests/test_smc.py"
+  - PYTHON_VERSION=2.7 FLOATX='float32' TESTCMD="--durations=10 --cov-append pymc3/tests/test_examples.py pymc3/tests/test_variational_inference.py pymc3/tests/test_updates.py"
+  - PYTHON_VERSION=2.7 FLOATX='float64' TESTCMD="--durations=10 --ignore=pymc3/tests/test_examples.py --cov-append --ignore=pymc3/tests/test_distributions_random.py --ignore=pymc3/tests/test_variational_inference.py --ignore=pymc3/tests/test_shared.py --ignore=pymc3/tests/test_smc.py --ignore=pymc3/tests/test_updates.py"
+  - PYTHON_VERSION=2.7 FLOATX='float64' RUN_PYLINT="true" TESTCMD="--durations=10 --cov-append pymc3/tests/test_distributions_random.py pymc3/tests/test_shared.py pymc3/tests/test_smc.py"
+  - PYTHON_VERSION=2.7 FLOATX='float64' TESTCMD="--durations=10 --cov-append pymc3/tests/test_examples.py pymc3/tests/test_variational_inference.py pymc3/tests/test_updates.py"
+  - PYTHON_VERSION=3.6 FLOATX='float64' TESTCMD="--durations=10 --cov-append --ignore=pymc3/tests/test_examples.py --ignore=pymc3/tests/test_distributions_random.py --ignore=pymc3/tests/test_variational_inference.py --ignore=pymc3/tests/test_shared.py --ignore=pymc3/tests/test_smc.py --ignore=pymc3/tests/test_updates.py"
+  - PYTHON_VERSION=3.6 FLOATX='float64' TESTCMD="--durations=10 --cov-append pymc3/tests/test_distributions_random.py pymc3/tests/test_shared.py pymc3/tests/test_smc.py"
+  - PYTHON_VERSION=3.6 FLOATX='float64' TESTCMD="--durations=10 --cov-append pymc3/tests/test_examples.py pymc3/tests/test_variational_inference.py pymc3/tests/test_updates.py"
 script:
   - . ./scripts/test.sh $TESTCMD
 

docs/source/examples.rst

@@ -12,6 +12,8 @@ Howto
    notebooks/sampler-stats.ipynb
    notebooks/Diagnosing_biased_Inference_with_Divergences.ipynb
    notebooks/posterior_predictive.ipynb
+   notebooks/model_comparison.ipynb
+   notebooks/model_averaging.ipynb
    notebooks/howto_debugging.ipynb
    notebooks/PyMC3_tips_and_heuristic.ipynb
    notebooks/LKJ.ipynb

docs/source/notebooks/model_averaging.ipynb

@@ -177,13 +177,13 @@ def highest_sampled_stage(self):
         -------
         stage number : int
         """
-        return max(self.stage_number(s) for s in glob(self.path('*')))
+        return max(self.stage_number(s) for s in glob(self.stage_path('*')))
 
     def atmip_path(self, stage_number):
         """Consistent naming for atmip params."""
         return os.path.join(self.stage_path(stage_number), 'atmip.params.pkl')
 
-    def load_atmip_params(self, stage_number):
+    def load_atmip_params(self, stage_number, model):
         """Load saved parameters from last sampled ATMIP stage.
 
         Parameters
@@ -196,8 +196,14 @@ def load_atmip_params(self, stage_number):
         else:
             prev = stage_number - 1
         pm._log.info('Loading parameters from completed stage {}'.format(prev))
-        with open(self.atmip_path(prev), 'rb') as buff:
-            return pickle.load(buff)
+
+        with model:
+            with open(self.atmip_path(prev), 'rb') as buff:
+                step = pickle.load(buff)
+
+        # update step stage to current stage
+        step.stage = stage_number
+        return step
 
     def dump_atmip_params(self, step):
         """Save atmip params to file."""
@@ -278,7 +284,7 @@ def recover_existing_results(self, stage, draws, step, n_jobs, model=None):
             # load incomplete stage results
             pm._log.info('Reloading existing results ...')
             mtrace = self.load_multitrace(stage, model=model)
-            if len(mtrace) > 0:
+            if len(mtrace.chains) > 0:
                 # continue sampling if traces exist
                 pm._log.info('Checking for corrupted files ...')
                 return self.check_multitrace(mtrace, draws=draws, n_chains=step.n_chains)

docs/source/notebooks/Model Comparison.ipynb renamed to docs/source/notebooks/model_comparison.ipynb

@@ -87,7 +87,7 @@ def getattr_value(self, val):
 
     def _repr_latex_(self, name=None, dist=None):
         return None
-                                                                     
+
 
 def TensorType(dtype, shape):
     return tt.TensorType(str(dtype), np.atleast_1d(shape) == 1)
@@ -123,6 +123,11 @@ def __init__(self, shape=(), dtype=None, defaults=('mode', ),
                 dtype = 'int64'
         if dtype != 'int16' and dtype != 'int64':
             raise TypeError('Discrete classes expect dtype to be int16 or int64.')
+
+        if kwargs.get('transform', None) is not None:
+            raise ValueError("Transformations for discrete distributions "
+                             "are not allowed.")
+
         super(Discrete, self).__init__(
             shape, dtype, defaults=defaults, *args, **kwargs)
 

pymc3/backends/smc_text.py

@@ -70,13 +70,6 @@ def __init__(self, dist, transform, *args, **kwargs):
             # force the last dim not broadcastable
             self.type = tt.TensorType(v.dtype, b)
 
-    def _repr_latex_(self, name=None, dist=None):
-        if name is None:
-            name = self.name
-        if dist is None:
-            dist = self.dist
-        return dist._repr_latex_(self, name=name, dist=dist)
-
     def logp(self, x):
         return (self.dist.logp(self.transform_used.backward(x)) +
                 self.transform_used.jacobian_det(x))

pymc3/distributions/distribution.py

@@ -0,0 +1,18 @@
+kcal.per.g,neocortex,log_mass
+0.490,0.552,0.668
+0.470,0.645,1.658
+0.560,0.645,1.681
+0.890,0.676,0.920
+0.920,0.688,-0.386
+0.800,0.589,-2.120
+0.460,0.617,-0.755
+0.710,0.603,-1.139
+0.680,0.700,0.438
+0.970,0.704,1.176
+0.840,0.734,2.510
+0.620,0.675,1.681
+0.540,0.713,3.569
+0.490,0.726,4.375
+0.480,0.702,3.707
+0.550,0.763,3.500
+0.710,0.755,4.006

pymc3/distributions/transforms.py

@@ -539,6 +539,7 @@ def Var(self, name, dist, data=None, total_size=None):
                                 name=name,
                                 orig_name=get_transformed_name(name, dist.transform)))
                 self.deterministics.append(var)
+                self.add_random_variable(var)
                 return var
         elif isinstance(data, dict):
             with self:
@@ -985,7 +986,7 @@ def Deterministic(name, var, model=None):
 
     Returns
     -------
-    n : var but with name name
+    var : var, with name attribute
     """
     model = modelcontext(model)
     var.name = model.name_for(name)
@@ -1009,6 +1010,7 @@ def Potential(name, var, model=None):
     model = modelcontext(model)
     var.name = model.name_for(name)
     model.potentials.append(var)
+    model.add_random_variable(var)
     return var
 
 

pymc3/examples/data/milk.csv

@@ -18,7 +18,7 @@
 import sys
 sys.setrecursionlimit(10000)
 
-__all__ = ['sample', 'iter_sample', 'sample_ppc', 'init_nuts']
+__all__ = ['sample', 'iter_sample', 'sample_ppc', 'sample_ppc_w', 'init_nuts']
 
 STEP_METHODS = (NUTS, HamiltonianMC, Metropolis, BinaryMetropolis,
                 BinaryGibbsMetropolis, Slice, CategoricalGibbsMetropolis)
@@ -484,14 +484,15 @@ def _update_start_vals(a, b, model):
 
     a.update({k: v for k, v in b.items() if k not in a})
 
+
 def sample_ppc(trace, samples=None, model=None, vars=None, size=None,
                random_seed=None, progressbar=True):
     """Generate posterior predictive samples from a model given a trace.
 
     Parameters
     ----------
     trace : backend, list, or MultiTrace
-        Trace generated from MCMC sampling
+        Trace generated from MCMC sampling.
     samples : int
         Number of posterior predictive samples to generate. Defaults to the
         length of `trace`
@@ -503,12 +504,19 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None,
     size : int
         The number of random draws from the distribution specified by the
         parameters in each sample of the trace.
+    random_seed : int
+        Seed for the random number generator.
+    progressbar : bool
+        Whether or not to display a progress bar in the command line. The
+        bar shows the percentage of completion, the sampling speed in
+        samples per second (SPS), and the estimated remaining time until
+        completion ("expected time of arrival"; ETA).
 
     Returns
     -------
     samples : dict
-        Dictionary with the variables as keys. The values corresponding
-        to the posterior predictive samples.
+        Dictionary with the variables as keys. The values corresponding to the
+        posterior predictive samples.
     """
     if samples is None:
         samples = len(trace)
@@ -521,18 +529,128 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None,
 
     seed(random_seed)
 
+    indices = randint(0, len(trace), samples)
     if progressbar:
-        indices = tqdm(randint(0, len(trace), samples), total=samples)
-    else:
-        indices = randint(0, len(trace), samples)
+        indices = tqdm(indices, total=samples)
+
+    try:
+        ppc = defaultdict(list)
+        for idx in indices:
+            param = trace[idx]
+            for var in vars:
+                ppc[var.name].append(var.distribution.random(point=param,
+                                                             size=size))
+
+    except KeyboardInterrupt:
+        pass
+
+    finally:
+        if progressbar:
+            indices.close()
+
+    return {k: np.asarray(v) for k, v in ppc.items()}
+
+
+def sample_ppc_w(traces, samples=None, models=None, size=None, weights=None,
+                 random_seed=None, progressbar=True):
+    """Generate weighted posterior predictive samples from a list of models and
+    a list of traces according to a set of weights.
+
+    Parameters
+    ----------
+    traces : list
+        List of traces generated from MCMC sampling. The number of traces should
+        be equal to the number of weights.
+    samples : int
+        Number of posterior predictive samples to generate. Defaults to the
+        length of the shorter trace in traces.
+    models : list
+        List of models used to generate the list of traces. The number of models
+        should be equal to the number of weights and the number of observed RVs
+        should be the same for all models.
+        By default a single model will be inferred from `with` context, in this
+        case results will only be meaningful if all models share the same
+        distributions for the observed RVs.
+    size : int
+        The number of random draws from the distributions specified by the
+        parameters in each sample of the trace.
+    weights: array-like
+        Individual weights for each trace. Default, same weight for each model.
+    random_seed : int
+        Seed for the random number generator.
+    progressbar : bool
+        Whether or not to display a progress bar in the command line. The
+        bar shows the percentage of completion, the sampling speed in
+        samples per second (SPS), and the estimated remaining time until
+        completion ("expected time of arrival"; ETA).
+
+    Returns
+    -------
+    samples : dict
+        Dictionary with the variables as keys. The values corresponding to the
+        posterior predictive samples from the weighted models.
+    """
+    seed(random_seed)
+
+    if models is None:
+        models = [modelcontext(models)] * len(traces)
+
+    if weights is None:
+        weights = [1] * len(traces)
+
+    if len(traces) != len(weights):
+        raise ValueError('The number of traces and weights should be the same')
+
+    if len(models) != len(weights):
+        raise ValueError('The number of models and weights should be the same')
+
+    lenght_morv = len(models[0].observed_RVs)
+    if not all(len(i.observed_RVs) == lenght_morv for i in models):
+        raise ValueError(
+            'The number of observed RVs should be the same for all models')
+
+    weights = np.asarray(weights)
+    p = weights / np.sum(weights)
+
+    min_tr = min([len(i) for i in traces])
+
+    n = (min_tr * p).astype('int')
+    # ensure n sum up to min_tr
+    idx = np.argmax(n)
+    n[idx] = n[idx] + min_tr - np.sum(n)
 
-    ppc = defaultdict(list)
-    for idx in indices:
-        param = trace[idx]
-        for var in vars:
+    trace = np.concatenate([np.random.choice(traces[i], j)
+                            for i, j in enumerate(n)])
+
+    variables = []
+    for i, m in enumerate(models):
+        variables.extend(m.observed_RVs * n[i])
+
+    len_trace = len(trace)
+
+    if samples is None:
+        samples = len_trace
+
+    indices = randint(0, len_trace, samples)
+
+    if progressbar:
+        indices = tqdm(indices, total=samples)
+
+    try:
+        ppc = defaultdict(list)
+        for idx in indices:
+            param = trace[idx]
+            var = variables[idx]
             ppc[var.name].append(var.distribution.random(point=param,
                                                          size=size))
 
+    except KeyboardInterrupt:
+        pass
+
+    finally:
+        if progressbar:
+            indices.close()
+
     return {k: np.asarray(v) for k, v in ppc.items()}
 
 

pymc3/model.py

@@ -520,7 +520,7 @@ def ATMIP_sample(n_steps, step=None, start=None, homepath=None, chain=0, stage=0
         step.stage = stage
         draws = 1
     else:
-        step = stage_handler.load_atmip_params(stage)
+        step = stage_handler.load_atmip_params(stage, model=model)
         draws = step.n_steps
 
     stage_handler.clean_directory(stage, None, rm_flag)

pymc3/sampling.py

@@ -7,6 +7,7 @@
 from pymc3.tests import models
 from pymc3.backends import base
 import pytest
+import theano
 
 
 class ModelBackendSetupTestCase(object):
@@ -227,6 +228,7 @@ def record_point(self, val):
         else:
             self.strace.record(point=point)
 
+    @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
     def test_standard_close(self):
         for idx in range(self.draws):
             self.record_point(idx)
@@ -266,13 +268,15 @@ class SelectionTestCase(ModelBackendSampledTestCase):
     - shape
     """
 
+    @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
     def test_get_values_default(self):
         for varname in self.test_point.keys():
             expected = np.concatenate([self.expected[chain][varname]
                                        for chain in [0, 1]])
             result = self.mtrace.get_values(varname)
             npt.assert_equal(result, expected)
 
+    @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
     def test_get_values_nocombine_burn_keyword(self):
         burn = 2
         for varname in self.test_point.keys():