Upstream rebase

Author: lucianopaz

README.rst

@@ -53,7 +53,7 @@ Learn Bayesian statistics with a book together with PyMC3:
 -  `PyMC3 port of the book "Doing Bayesian Data Analysis" by John Kruschke <https://github.com/aloctavodia/Doing_bayesian_data_analysis>`__ as well as the `second edition <https://github.com/JWarmenhoven/DBDA-python>`__: Principled introduction to Bayesian data analysis.
 -  `PyMC3 port of the book "Statistical Rethinking A Bayesian Course with Examples in R and Stan" by Richard McElreath <https://github.com/pymc-devs/resources/tree/master/Rethinking>`__
 -  `PyMC3 port of the book "Bayesian Cognitive Modeling" by Michael Lee and EJ Wagenmakers <https://github.com/pymc-devs/resources/tree/master/BCM>`__: Focused on using Bayesian statistics in cognitive modeling.
--  `Bayesian Analysis with Python by Osvaldo Martin <https://www.packtpub.com/big-data-and-business-intelligence/bayesian-analysis-python>`__ (and `errata <https://github.com/aloctavodia/BAP>`__): Great introductory book.
+-  `Bayesian Analysis with Python  <https://www.packtpub.com/big-data-and-business-intelligence/bayesian-analysis-python-second-edition>`__ (second edition) by Osvaldo Martin: Great introductory book. (`code <https://github.com/aloctavodia/BAP>`__ and errata).
 
 PyMC3 talks
 -----------

pymc3/distributions/multivariate.py

@@ -775,7 +775,7 @@ def WishartBartlett(name, S, nu, is_cholesky=False, return_cholesky=False, testv
     -----
     This is not a standard Distribution class but follows a similar
     interface. Besides the Wishart distribution, it will add RVs
-    c and z to your model which make up the matrix.
+    name_c and name_z to your model which make up the matrix.
 
     This distribution is usually a bad idea to use as a prior for multivariate
     normal. You should instead use LKJCholeskyCov or LKJCorr.
@@ -797,11 +797,11 @@ def WishartBartlett(name, S, nu, is_cholesky=False, return_cholesky=False, testv
         diag_testval = None
         tril_testval = None
 
-    c = tt.sqrt(ChiSquared('c', nu - np.arange(2, 2 + n_diag), shape=n_diag,
+    c = tt.sqrt(ChiSquared('%s_c' % name, nu - np.arange(2, 2 + n_diag), shape=n_diag,
                            testval=diag_testval))
-    pm._log.info('Added new variable c to model diagonal of Wishart.')
-    z = Normal('z', 0., 1., shape=n_tril, testval=tril_testval)
-    pm._log.info('Added new variable z to model off-diagonals of Wishart.')
+    pm._log.info('Added new variable %s_c to model diagonal of Wishart.' % name)
+    z = Normal('%s_z' % name, 0., 1., shape=n_tril, testval=tril_testval)
+    pm._log.info('Added new variable %s_z to model off-diagonals of Wishart.' % name)
     # Construct A matrix
     A = tt.zeros(S.shape, dtype=np.float32)
     A = tt.set_subtensor(A[diag_idx], c)

pymc3/sampling.py

@@ -6,7 +6,6 @@
 
 from six import integer_types
 from joblib import Parallel, delayed
-from tempfile import mkdtemp
 import numpy as np
 import theano.gradient as tg
 
@@ -245,7 +244,8 @@ def sample(draws=500, step=None, init='auto', n_init=200000, start=None, trace=N
     chains : int
         The number of chains to sample. Running independent chains is important for some
         convergence statistics and can also reveal multiple modes in the posterior. If `None`,
-        then set to either `cores` or 2, whichever is larger. For SMC the default value is 100.
+        then set to either `cores` or 2, whichever is larger. For SMC the number of chains is the
+        number of draws.
     cores : int
         The number of chains to run in parallel. If `None`, set to the number of CPUs in the
         system, but at most 4 (for 'SMC' defaults to 1). Keep in mind that some chains might
@@ -323,7 +323,6 @@ def sample(draws=500, step=None, init='auto', n_init=200000, start=None, trace=N
     if isinstance(step, pm.step_methods.smc.SMC):
         if step_kwargs is None:
             step_kwargs = {}
-        test_folder = mkdtemp(prefix='SMC_TEST')
         trace = smc.sample_smc(draws=draws,
                                step=step,
                                progressbar=progressbar,
@@ -1523,4 +1522,4 @@ def init_nuts(init='auto', chains=1, n_init=500000, model=None,
 
     step = pm.NUTS(potential=potential, model=model, **kwargs)
 
-    return start, step
+    return start, step

pymc3/step_methods/smc.py

@@ -5,7 +5,6 @@
 import theano
 import pymc3 as pm
 from tqdm import tqdm
-import warnings
 
 from .arraystep import metrop_select
 from .metropolis import MultivariateNormalProposal
@@ -17,8 +16,6 @@
 
 __all__ = ["SMC", "sample_smc"]
 
-proposal_dists = {"MultivariateNormal": MultivariateNormalProposal}
-
 
 class SMC:
     """
@@ -27,20 +24,24 @@ class SMC:
     Parameters
     ----------
     n_steps : int
-        The number of steps of a Markov Chain. If `tune == True` `n_steps` will be used for
-        the first stage, and the number of steps of the other states will be determined
+        The number of steps of a Markov Chain. If `tune_steps == True` `n_steps` will be used for
+        the first stage and the number of steps of the other stages will be determined
         automatically based on the acceptance rate and `p_acc_rate`.
+        The number of steps will never be larger than `n_steps`.
     scaling : float
         Factor applied to the proposal distribution i.e. the step size of the Markov Chain. Only
-        works if `tune == False` otherwise is determined automatically
+        works if `tune_scaling == False` otherwise is determined automatically.
     p_acc_rate : float
-        Probability of not accepting a Markov Chain proposal. Used to compute `n_steps` when
-        `tune == True`. It should be between 0 and 1.
-    proposal_name :
-        Type of proposal distribution. Currently the only valid option is `MultivariateNormal`.
+        Used to compute `n_steps` when `tune_steps == True`. The higher the value of `p_acc_rate`
+        the higher the number of steps computed automatically. Defaults to 0.99. It should be
+        between 0 and 1.
+    tune_scaling : bool
+        Whether to compute the scaling automatically or not. Defaults to True
+    tune_steps : bool
+        Whether to compute the number of steps automatically or not. Defaults to True
     threshold : float
         Determines the change of beta from stage to stage, i.e.indirectly the number of stages,
-        the higher the value of threshold the higher the number of stages. Defaults to 0.5.
+        the higher the value of `threshold` the higher the number of stages. Defaults to 0.5.
         It should be between 0 and 1.
     model : :class:`pymc3.Model`
         Optional model for sampling step. Defaults to None (taken from context).
@@ -63,17 +64,18 @@ def __init__(
         self,
         n_steps=25,
         scaling=1.0,
-        p_acc_rate=0.01,
-        tune=True,
-        proposal_name="MultivariateNormal",
+        p_acc_rate=0.99,
+        tune_scaling=True,
+        tune_steps=True,
         threshold=0.5,
     ):
 
         self.n_steps = n_steps
+        self.max_steps = n_steps
         self.scaling = scaling
-        self.p_acc_rate = p_acc_rate
-        self.tune = tune
-        self.proposal = proposal_dists[proposal_name]
+        self.p_acc_rate = 1 - p_acc_rate
+        self.tune_scaling = tune_scaling
+        self.tune_steps = tune_steps
         self.threshold = threshold
 
 
@@ -95,15 +97,15 @@ def sample_smc(draws=5000, step=None, progressbar=False, model=None, random_seed
     random_seed : int
         random seed
     """
-    warnings.warn("Warning: SMC is experimental, hopefully it will be ready for PyMC 3.6")
     model = modelcontext(model)
 
     if random_seed != -1:
         np.random.seed(random_seed)
 
-    beta = 0
+    beta = 0.
     stage = 0
-    acc_rate = 1
+    acc_rate = 1.
+    proposed = draws * step.n_steps
     model.marginal_likelihood = 1
     variables = inputvars(model.vars)
     discrete = np.concatenate([[v.dtype in pm.discrete_types] * (v.dsize or 1) for v in variables])
@@ -128,31 +130,32 @@ def sample_smc(draws=5000, step=None, progressbar=False, model=None, random_seed
 
         # compute proposal distribution based on weights
         covariance = _calc_covariance(posterior, weights)
-        proposal = step.proposal(covariance)
-
-        # compute scaling and number of Markov chains steps (optional), based on previous
-        # acceptance rate
-        if step.tune and stage > 0:
-            if acc_rate == 0:
-                acc_rate = 1.0 / step.n_steps
-            step.scaling = _tune(acc_rate)
-            step.n_steps = 1 + int(np.log(step.p_acc_rate) / np.log(1 - acc_rate))
+        proposal = MultivariateNormalProposal(covariance)
+
+        # compute scaling (optional) and number of Markov chains steps (optional), based on the
+        # acceptance rate of the previous stage
+        if (step.tune_scaling or step.tune_steps) and stage > 0:
+            if step.tune_scaling:
+                step.scaling = _tune(acc_rate)
+            if step.tune_steps:
+                acc_rate = max(1. / proposed, acc_rate)
+                step.n_steps = min(
+                    step.max_steps, 1 + int(np.log(step.p_acc_rate) / np.log(1 - acc_rate))
+                )
 
         pm._log.info(
-            "Stage: {:d} Beta: {:f} Steps: {:d} Acc: {:f}".format(
-                stage, beta, step.n_steps, acc_rate
-            )
+            "Stage: {:d} Beta: {:f} Steps: {:d}".format(stage, beta, step.n_steps, acc_rate)
         )
         # Apply Metropolis kernel (mutation)
-        proposed = 0.0
-        accepted = 0.0
+        proposed = draws * step.n_steps
+        accepted = 0.
         priors = np.array([prior_logp(sample) for sample in posterior]).squeeze()
         tempered_post = priors + likelihoods * beta
         for draw in tqdm(range(draws), disable=not progressbar):
             old_tempered_post = tempered_post[draw]
             q_old = posterior[draw]
             deltas = np.squeeze(proposal(step.n_steps) * step.scaling)
-            for n_step in range(0, step.n_steps):
+            for n_step in range(step.n_steps):
                 delta = deltas[n_step]
 
                 if any_discrete:
@@ -170,10 +173,9 @@ def sample_smc(draws=5000, step=None, progressbar=False, model=None, random_seed
 
                 q_old, accept = metrop_select(new_tempered_post - old_tempered_post, q_new, q_old)
                 if accept:
-                    accepted += accept
+                    accepted += 1
                     posterior[draw] = q_old
                     old_tempered_post = new_tempered_post
-                proposed += 1.0
 
         acc_rate = accepted / proposed
         stage += 1
@@ -258,7 +260,7 @@ def _calc_covariance(posterior, weights):
     """
     cov = np.cov(posterior, aweights=weights.ravel(), bias=False, rowvar=0)
     if np.isnan(cov).any() or np.isinf(cov).any():
-        raise ValueError('Sample covariances not valid! Likely "chains" is too small!')
+        raise ValueError('Sample covariances not valid! Likely "draws" is too small!')
     return np.atleast_2d(cov)