Merge pull request #2339 from hvasbath/smc_model_llk

Remove necessity to define likelihood variable inside model, is done …

Author: Junpeng Lao

docs/source/notebooks/SMC2_gaussians.ipynb

@@ -50,8 +50,12 @@ def __init__(self, s):
         self.chol = scipy.linalg.cholesky(s, lower=True)
 
     def __call__(self, num_draws=None):
-        b = np.random.randn(self.n)
-        return np.dot(self.chol, b)
+        if num_draws is not None:
+            b = np.random.randn(self.n, num_draws)
+            return np.dot(self.chol, b).T
+        else:
+            b = np.random.randn(self.n)
+            return np.dot(self.chol, b)
 
 
 class Metropolis(ArrayStepShared):

pymc3/step_methods/metropolis.py

@@ -25,33 +25,15 @@
 from ..theanof import inputvars, make_shared_replacements, join_nonshared_inputs
 import numpy.random as nr
 
+from .metropolis import MultivariateNormalProposal
 from .arraystep import metrop_select
 from ..backends import smc_text as atext
 
-__all__ = ['SMC', 'ATMIP_sample']
+__all__ = ['SMC', 'sample_smc']
 
 EXPERIMENTAL_WARNING = "Warning: SMC is an experimental step method, and not yet"\
     " recommended for use in PyMC3!"
 
-
-class Proposal(object):
-    """Proposal distributions modified from pymc3 to initially create all the
-    Proposal steps without repeated execution of the RNG - significant speedup!
-
-    Parameters
-    ----------
-    s : :class:`numpy.ndarray`
-    """
-    def __init__(self, s):
-        self.s = np.atleast_1d(s)
-
-
-class MultivariateNormalProposal(Proposal):
-    def __call__(self, num_draws=None):
-        return np.random.multivariate_normal(
-                mean=np.zeros(self.s.shape[0]), cov=self.s, size=num_draws)
-
-
 proposal_dists = {
     'MultivariateNormal': MultivariateNormalProposal,
         }
@@ -147,6 +129,13 @@ def __init__(self, vars=None, out_vars=None, n_chains=100, scaling=1., covarianc
         vars = inputvars(vars)
 
         if out_vars is None:
+            if not any(likelihood_name == RV.name for RV in model.unobserved_RVs):
+                with model:
+                    llk = pm.Deterministic(likelihood_name, model.logpt)
+            else:
+                raise ValueError(
+                    'The model likelihood name is already being used by a RV!')
+
             out_vars = model.unobserved_RVs
 
         out_varnames = [out_var.name for out_var in out_vars]
@@ -419,9 +408,9 @@ def resample(self):
         return outindx
 
 
-def ATMIP_sample(n_steps, step=None, start=None, homepath=None, chain=0, stage=0, n_jobs=1,
+def sample_smc(n_steps, step=None, start=None, homepath=None, chain=0, stage=0, n_jobs=1,
                  tune=None, progressbar=False, model=None, random_seed=-1, rm_flag=False):
-    """(C)ATMIP sampling algorithm (Cascading - (C) not always relevant)
+    """Sequential Monte Carlo sampling
 
     Samples the solution space with n_chains of Metropolis chains, where each
     chain has n_steps iterations. Once finished, the sampled traces are
@@ -524,25 +513,8 @@ def ATMIP_sample(n_steps, step=None, start=None, homepath=None, chain=0, stage=0
         draws = step.n_steps
 
     stage_handler.clean_directory(stage, None, rm_flag)
-    with model:
-        chains = stage_handler.recover_existing_results(stage, draws, step, n_jobs)
-        if chains is not None:
-            rest = len(chains) % n_jobs
-            if rest > 0:
-                pm._log.info('Fixing %i chains ...' % rest)
-                chains, rest_chains = chains[:-rest], chains[-rest:]
-                # process traces that are not a multiple of n_jobs
-                sample_args = {
-                    'draws': draws,
-                    'step': step,
-                    'stage_path': stage_handler.stage_path(stage),
-                    'progressbar': progressbar,
-                    'model': model,
-                    'n_jobs': rest,
-                    'chains': rest_chains}
 
-                _iter_parallel_chains(**sample_args)
-                pm._log.info('Back to normal!')
+    chains = stage_handler.recover_existing_results(stage, draws, step, n_jobs)
 
     with model:
         while step.beta < 1:
@@ -556,7 +528,7 @@ def ATMIP_sample(n_steps, step=None, start=None, homepath=None, chain=0, stage=0
             pm._log.info('Beta: %f Stage: %i' % (step.beta, step.stage))
 
             # Metropolis sampling intermediate stages
-            chains = stage_handler.clean_directory(stage, chains, rm_flag)
+            chains = stage_handler.clean_directory(step.stage, chains, rm_flag)
             sample_args = {
                     'draws': draws,
                     'step': step,

pymc3/step_methods/smc.py

@@ -51,26 +51,25 @@ def two_gaussians(x):
                            upper=2. * np.ones_like(mu1),
                            testval=-1. * np.ones_like(mu1),
                            transform=None)
-            like = pm.Deterministic('like', two_gaussians(X))
-            llk = pm.Potential('like_potential', like)
+            llk = pm.Potential('muh', two_gaussians(X))
+
+        self.step = smc.SMC(
+            n_chains=self.n_chains,
+            tune_interval=self.tune_interval,
+            model=self.ATMIP_test)
 
         self.muref = mu1
 
-    @pytest.mark.parametrize('n_jobs', [1, 2])
-    def test_sample_n_core(self, n_jobs):
+    @pytest.mark.parametrize(['n_jobs', 'stage'], [[1, 0], [2, 6]])
+    def test_sample_n_core(self, n_jobs, stage):
 
         def last_sample(x):
             return x[(self.n_steps - 1)::self.n_steps]
 
-        step = smc.SMC(
-            n_chains=self.n_chains,
-            tune_interval=self.tune_interval,
-            model=self.ATMIP_test,
-            likelihood_name=self.ATMIP_test.deterministics[0].name)
-
-        mtrace = smc.ATMIP_sample(
+        mtrace = smc.sample_smc(
             n_steps=self.n_steps,
-            step=step,
+            step=self.step,
+            stage=stage,
             n_jobs=n_jobs,
             progressbar=True,
             homepath=self.test_folder,

pymc3/tests/test_smc.py

@@ -5,7 +5,7 @@
 from .checks import close_to
 from .models import simple_categorical, mv_simple, mv_simple_discrete, simple_2model, mv_prior_simple
 from pymc3.sampling import assign_step_methods, sample
-from pymc3.model import Model, Deterministic
+from pymc3.model import Model
 from pymc3.step_methods import (NUTS, BinaryGibbsMetropolis, CategoricalGibbsMetropolis,
                                 Metropolis, Slice, CompoundStep, NormalProposal,
                                 MultivariateNormalProposal, HamiltonianMC,
@@ -163,8 +163,7 @@ def check_trace(self, step_method):
         with Model():
             x = Normal('x', mu=0, sd=1)
             if step_method.__name__ == 'SMC':
-                Deterministic('like', - 0.5 * tt.log(2 * np.pi) - 0.5 * x.T.dot(x))
-                trace = smc.ATMIP_sample(n_steps=n_steps, step=step_method(random_seed=1),
+                trace = smc.sample_smc(n_steps=n_steps, step=step_method(random_seed=1),
                                          n_jobs=1, progressbar=False,
                                          homepath=self.temp_dir)
             else: