NegativeBinomial distribution sample generator wrapper

RELEASE-NOTES.md

@@ -18,6 +18,8 @@
 - Deprecated `sd` in version 3.7 has been replaced by `sigma` now raises `DepreciationWarning` on using `sd` in continuous, mixed and timeseries distributions. (see #3837 and #3688).
 - In named models, `pm.Data` objects now get model-relative names (see [#3843](https://github.com/pymc-devs/pymc3/pull/3843)).
 - `pm.sample` now takes 1000 draws and 1000 tuning samples by default, instead of 500 previously (see [#3855](https://github.com/pymc-devs/pymc3/pull/3855)).
+- Dropped the outdated 'nuts' initialization method for `pm.sample` (see [#3863](https://github.com/pymc-devs/pymc3/pull/3863)).
+- Moved argument division out of `NegativeBinomial` `random` method. Fixes [#3864](https://github.com/pymc-devs/pymc3/issues/3864) in the style of [#3509](https://github.com/pymc-devs/pymc3/pull/3509).
 
 ## PyMC3 3.8 (November 29 2019)
 

pymc3/distributions/discrete.py

@@ -673,12 +673,24 @@ def random(self, point=None, size=None):
         array
         """
         mu, alpha = draw_values([self.mu, self.alpha], point=point, size=size)
-        g = generate_samples(stats.gamma.rvs, alpha, scale=mu / alpha,
+        g = generate_samples(self._random, mu=mu, alpha=alpha,
                              dist_shape=self.shape,
                              size=size)
         g[g == 0] = np.finfo(float).eps  # Just in case
         return np.asarray(stats.poisson.rvs(g)).reshape(g.shape)
 
+    def _random(self, mu, alpha, size):
+        """ Wrapper around stats.gamma.rvs that converts NegativeBinomial's
+        parametrization to scipy.gamma. All parameter arrays should have
+        been broadcasted properly by generate_samples at this point and size is
+        the scipy.rvs representation.
+        """
+        return stats.gamma.rvs(
+            a=alpha,
+            scale=mu / alpha,
+            size=size,
+        )
+
     def logp(self, value):
         """
         Calculate log-probability of NegativeBinomial distribution at specified value.

pymc3/sampling.py

@@ -274,16 +274,14 @@ def sample(
         * advi: Run ADVI to estimate posterior mean and diagonal mass matrix.
         * advi_map: Initialize ADVI with MAP and use MAP as starting point.
         * map: Use the MAP as starting point. This is discouraged.
-        * nuts: Run NUTS and estimate posterior mean and mass matrix from the trace.
         * adapt_full: Adapt a dense mass matrix using the sample covariances
     step: function or iterable of functions
         A step function or collection of functions. If there are variables without step methods,
         step methods for those variables will be assigned automatically.  By default the NUTS step
         method will be used, if appropriate to the model; this is a good default for beginning
         users.
     n_init: int
-        Number of iterations of initializer. Only works for 'nuts' and 'ADVI'.
-        If 'ADVI', number of iterations, if 'nuts', number of draws.
+        Number of iterations of initializer. Only works for 'ADVI' init methods.
     start: dict, or array of dict
         Starting point in parameter space (or partial point)
         Defaults to ``trace.point(-1))`` if there is a trace provided and model.test_point if not
@@ -1569,9 +1567,9 @@ def sample_posterior_predictive(
         nchain = 1
 
     if keep_size and samples is not None:
-        raise IncorrectArgumentsError("Should not specify both keep_size and samples argukments")
+        raise IncorrectArgumentsError("Should not specify both keep_size and samples arguments")
     if keep_size and size is not None:
-        raise IncorrectArgumentsError("Should not specify both keep_size and size argukments")
+        raise IncorrectArgumentsError("Should not specify both keep_size and size arguments")
 
     if samples is None:
         if isinstance(trace, MultiTrace):
@@ -1865,14 +1863,11 @@ def init_nuts(
         * advi: Run ADVI to estimate posterior mean and diagonal mass matrix.
         * advi_map: Initialize ADVI with MAP and use MAP as starting point.
         * map: Use the MAP as starting point. This is discouraged.
-        * nuts: Run NUTS and estimate posterior mean and mass matrix from
-          the trace.
         * adapt_full: Adapt a dense mass matrix using the sample covariances
     chains: int
         Number of jobs to start.
     n_init: int
-        Number of iterations of initializer
-        If 'ADVI', number of iterations, if 'nuts', number of draws.
+        Number of iterations of initializer. Only works for 'ADVI' init methods.
     model: Model (optional if in ``with`` context)
     progressbar: bool
         Whether or not to display a progressbar for advi sampling.
@@ -2001,13 +1996,6 @@ def init_nuts(
         cov = pm.find_hessian(point=start)
         start = [start] * chains
         potential = quadpotential.QuadPotentialFull(cov)
-    elif init == "nuts":
-        init_trace = pm.sample(
-            draws=n_init, step=pm.NUTS(), tune=n_init // 2, random_seed=random_seed
-        )
-        cov = np.atleast_1d(pm.trace_cov(init_trace))
-        start = list(np.random.choice(init_trace, chains))
-        potential = quadpotential.QuadPotentialFull(cov)
     elif init == "adapt_full":
         start = [model.test_point] * chains
         mean = np.mean([model.dict_to_array(vals) for vals in start], axis=0)

pymc3/tests/test_distributions_random.py

@@ -1176,6 +1176,31 @@ def sample_prior(
         with model:
             return pm.sample_prior_predictive(prior_samples)
 
+    @pytest.mark.parametrize(
+        ["prior_samples", "shape", "mu", "alpha"],
+        [
+            [10, (3,), (None, tuple()), (None, (3,))],
+            [10, (3,), (None, (3,)), (None, tuple())],
+            [10, (4, 3,), (None, (3,)), (None, (3,))],
+            [10, (4, 3,), (None, (3,)), (None, (4, 3))],
+        ],
+        ids=str,
+    )
+    def test_NegativeBinomial(
+        self,
+        prior_samples,
+        shape,
+        mu,
+        alpha,
+    ):
+        prior = self.sample_prior(
+            distribution=pm.NegativeBinomial,
+            shape=shape,
+            nested_rvs_info=dict(mu=mu, alpha=alpha),
+            prior_samples=prior_samples,
+        )
+        assert prior["target"].shape == (prior_samples,) + shape
+
     @pytest.mark.parametrize(
         ["prior_samples", "shape", "psi", "mu", "alpha"],
         [

pymc3/tests/test_sampling.py

@@ -87,7 +87,7 @@ def test_sample(self):
 
     def test_sample_init(self):
         with self.model:
-            for init in ("advi", "advi_map", "map", "nuts"):
+            for init in ("advi", "advi_map", "map"):
                 pm.sample(
                     init=init,
                     tune=0,
@@ -675,7 +675,6 @@ def test_sample_posterior_predictive_w(self):
         "advi+adapt_diag_grad",
         "map",
         "advi_map",
-        "nuts",
     ],
 )
 def test_exec_nuts_init(method):