diff --git a/pymc3/distributions/bound.py b/pymc3/distributions/bound.py index 5b99002d2c..c44684b061 100644 --- a/pymc3/distributions/bound.py +++ b/pymc3/distributions/bound.py @@ -68,17 +68,17 @@ def random(self, point=None, size=None): if self.lower is None and self.upper is None: return self._wrapped.random(point=point, size=size) elif self.lower is not None and self.upper is not None: - lower, upper = draw_values([self.lower, self.upper], point=point) + lower, upper = draw_values([self.lower, self.upper], point=point, size=size) return generate_samples(self._random, lower, upper, point, dist_shape=self.shape, size=size) elif self.lower is not None: - lower = draw_values([self.lower], point=point) + lower = draw_values([self.lower], point=point, size=size) return generate_samples(self._random, lower, np.inf, point, dist_shape=self.shape, size=size) else: - upper = draw_values([self.upper], point=point) + upper = draw_values([self.upper], point=point, size=size) return generate_samples(self._random, -np.inf, upper, point, dist_shape=self.shape, size=size) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 9272ad42d3..ca7eac1391 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -305,7 +305,7 @@ def __init__(self, mu=0, sd=None, tau=None, **kwargs): def random(self, point=None, size=None): mu, tau, _ = draw_values([self.mu, self.tau, self.sd], - point=point) + point=point, size=size) return generate_samples(stats.norm.rvs, loc=mu, scale=tau**-0.5, dist_shape=self.shape, size=size) @@ -406,7 +406,7 @@ def __init__(self, sd=None, tau=None, *args, **kwargs): assert_negative_support(sd, 'sd', 'HalfNormal') def random(self, point=None, size=None): - sd = draw_values([self.sd], point=point)[0] + sd = draw_values([self.sd], point=point, size=size)[0] return generate_samples(stats.halfnorm.rvs, loc=0., scale=sd, dist_shape=self.shape, size=size) @@ -547,7 +547,7 @@ def _random(self, mu, lam, alpha, size=None): def random(self, point=None, size=None): mu, lam, alpha = draw_values([self.mu, self.lam, self.alpha], - point=point) + point=point, size=size) return generate_samples(self._random, mu, lam, alpha, dist_shape=self.shape, @@ -672,7 +672,7 @@ def get_alpha_beta(self, alpha=None, beta=None, mu=None, sd=None): def random(self, point=None, size=None): alpha, beta = draw_values([self.alpha, self.beta], - point=point) + point=point, size=size) return generate_samples(stats.beta.rvs, alpha, beta, dist_shape=self.shape, size=size) @@ -2239,7 +2239,7 @@ def __init__(self, nu=None, sd=None, *args, **kwargs): def random(self, point=None, size=None, repeat=None): nu, sd = draw_values([self.nu, self.sd], - point=point) + point=point, size=size) return generate_samples(stats.rice.rvs, b=nu, scale=sd, loc=0, dist_shape=self.shape, size=size) diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index df4ac475ad..58ae48df2b 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -166,7 +166,7 @@ def _random(self, alpha, beta, n, size=None): def random(self, point=None, size=None): alpha, beta, n = \ - draw_values([self.alpha, self.beta, self.n], point=point) + draw_values([self.alpha, self.beta, self.n], point=point, size=size) return generate_samples(self._random, alpha=alpha, beta=beta, n=n, dist_shape=self.shape, size=size) @@ -247,7 +247,7 @@ def __init__(self, p=None, logit_p=None, *args, **kwargs): self.mode = tt.cast(tround(self.p), 'int8') def random(self, point=None, size=None): - p = draw_values([self.p], point=point)[0] + p = draw_values([self.p], point=point, size=size)[0] return generate_samples(stats.bernoulli.rvs, p, dist_shape=self.shape, size=size) @@ -343,7 +343,7 @@ def _random(self, q, beta, size=None): return np.ceil(np.power(np.log(1 - p) / np.log(q), 1. / beta)) - 1 def random(self, point=None, size=None): - q, beta = draw_values([self.q, self.beta], point=point) + q, beta = draw_values([self.q, self.beta], point=point, size=size) return generate_samples(self._random, q, beta, dist_shape=self.shape, @@ -410,7 +410,7 @@ def __init__(self, mu, *args, **kwargs): self.mode = tt.floor(mu).astype('int32') def random(self, point=None, size=None): - mu = draw_values([self.mu], point=point)[0] + mu = draw_values([self.mu], point=point, size=size)[0] return generate_samples(stats.poisson.rvs, mu, dist_shape=self.shape, size=size) @@ -490,7 +490,7 @@ def __init__(self, mu, alpha, *args, **kwargs): self.mode = tt.floor(mu).astype('int32') def random(self, point=None, size=None): - mu, alpha = draw_values([self.mu, self.alpha], point=point) + mu, alpha = draw_values([self.mu, self.alpha], point=point, size=size) g = generate_samples(stats.gamma.rvs, alpha, scale=mu / alpha, dist_shape=self.shape, size=size) @@ -564,7 +564,7 @@ def __init__(self, p, *args, **kwargs): self.mode = 1 def random(self, point=None, size=None): - p = draw_values([self.p], point=point)[0] + p = draw_values([self.p], point=point, size=size)[0] return generate_samples(np.random.geometric, p, dist_shape=self.shape, size=size) @@ -636,7 +636,7 @@ def _random(self, lower, upper, size=None): return samples def random(self, point=None, size=None): - lower, upper = draw_values([self.lower, self.upper], point=point) + lower, upper = draw_values([self.lower, self.upper], point=point, size=size) return generate_samples(self._random, lower, upper, dist_shape=self.shape, @@ -714,7 +714,7 @@ def random_choice(k, *args, **kwargs): else: return np.random.choice(k, *args, **kwargs) - p, k = draw_values([self.p, self.k], point=point) + p, k = draw_values([self.p, self.k], point=point, size=size) return generate_samples(partial(random_choice, np.arange(k)), p=p, broadcast_shape=p.shape[:-1] or (1,), @@ -764,7 +764,7 @@ def __init__(self, c, *args, **kwargs): self.mean = self.median = self.mode = self.c = c = tt.as_tensor_variable(c) def random(self, point=None, size=None): - c = draw_values([self.c], point=point)[0] + c = draw_values([self.c], point=point, size=size)[0] dtype = np.array(c).dtype def _random(c, dtype=dtype, size=None): @@ -845,7 +845,7 @@ def __init__(self, psi, theta, *args, **kwargs): self.mode = self.pois.mode def random(self, point=None, size=None): - theta, psi = draw_values([self.theta, self.psi], point=point) + theta, psi = draw_values([self.theta, self.psi], point=point, size=size) g = generate_samples(stats.poisson.rvs, theta, dist_shape=self.shape, size=size) @@ -938,7 +938,7 @@ def __init__(self, psi, n, p, *args, **kwargs): self.mode = self.bin.mode def random(self, point=None, size=None): - n, p, psi = draw_values([self.n, self.p, self.psi], point=point) + n, p, psi = draw_values([self.n, self.p, self.psi], point=point, size=size) g = generate_samples(stats.binom.rvs, n, p, dist_shape=self.shape, size=size) @@ -1056,7 +1056,7 @@ def __init__(self, psi, mu, alpha, *args, **kwargs): def random(self, point=None, size=None): mu, alpha, psi = draw_values( - [self.mu, self.alpha, self.psi], point=point) + [self.mu, self.alpha, self.psi], point=point, size=size) g = generate_samples(stats.gamma.rvs, alpha, scale=mu / alpha, dist_shape=self.shape, size=size) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index 997d0f5257..bdb7e266bd 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -177,8 +177,8 @@ def __init__(self, shape=(), dtype=None, defaults=('median', 'mean', 'mode'), class DensityDist(Distribution): """Distribution based on a given log density function. - - A distribution with the passed log density function is created. + + A distribution with the passed log density function is created. Requires a custom random function passed as kwarg `random` to enable sampling. @@ -200,7 +200,7 @@ def __init__(self, logp, shape=(), dtype=None, testval=0, random=None, *args, ** shape, dtype, testval, *args, **kwargs) self.logp = logp self.rand = random - + def random(self, *args, **kwargs): if self.rand is not None: return self.rand(*args, **kwargs) @@ -345,9 +345,7 @@ def _draw_value(param, point=None, givens=None, size=None): size : int, optional Number of samples """ - if isinstance(param, numbers.Number): - return param - elif isinstance(param, np.ndarray): + if isinstance(param, (numbers.Number, np.ndarray)): return param elif isinstance(param, tt.TensorConstant): return param.value @@ -357,7 +355,7 @@ def _draw_value(param, point=None, givens=None, size=None): if point and hasattr(param, 'model') and param.name in point: return point[param.name] elif hasattr(param, 'random') and param.random is not None: - return param.random(point=point, size=None) + return param.random(point=point, size=size) else: if givens: variables, values = list(zip(*givens)) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index 0a8f6d3001..d1786009a5 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -234,7 +234,7 @@ def random(self, point=None, size=None): size = [size] if self._cov_type == 'cov': - mu, cov = draw_values([self.mu, self.cov], point=point) + mu, cov = draw_values([self.mu, self.cov], point=point, size=size) if mu.shape[-1] != cov.shape[-1]: raise ValueError("Shapes for mu and cov don't match") @@ -246,7 +246,7 @@ def random(self, point=None, size=None): return np.nan * np.zeros(size) return dist.rvs(size) elif self._cov_type == 'chol': - mu, chol = draw_values([self.mu, self.chol_cov], point=point) + mu, chol = draw_values([self.mu, self.chol_cov], point=point, size=size) if mu.shape[-1] != chol[0].shape[-1]: raise ValueError("Shapes for mu and chol don't match") @@ -254,7 +254,7 @@ def random(self, point=None, size=None): standard_normal = np.random.standard_normal(size) return mu + np.dot(standard_normal, chol.T) else: - mu, tau = draw_values([self.mu, self.tau], point=point) + mu, tau = draw_values([self.mu, self.tau], point=point, size=size) if mu.shape[-1] != tau[0].shape[-1]: raise ValueError("Shapes for mu and tau don't match") @@ -338,15 +338,15 @@ def __init__(self, nu, Sigma=None, mu=None, cov=None, tau=None, chol=None, self.mean = self.median = self.mode = self.mu = self.mu def random(self, point=None, size=None): - nu, mu = draw_values([self.nu, self.mu], point=point) + nu, mu = draw_values([self.nu, self.mu], point=point, size=size) if self._cov_type == 'cov': - cov, = draw_values([self.cov], point=point) + cov, = draw_values([self.cov], point=point, size=size) dist = MvNormal.dist(mu=np.zeros_like(mu), cov=cov) elif self._cov_type == 'tau': - tau, = draw_values([self.tau], point=point) + tau, = draw_values([self.tau], point=point, size=size) dist = MvNormal.dist(mu=np.zeros_like(mu), tau=tau) else: - chol, = draw_values([self.chol_cov], point=point) + chol, = draw_values([self.chol_cov], point=point, size=size) dist = MvNormal.dist(mu=np.zeros_like(mu), chol=chol) samples = dist.random(point, size) @@ -422,7 +422,7 @@ def __init__(self, a, transform=transforms.stick_breaking, np.nan) def random(self, point=None, size=None): - a = draw_values([self.a], point=point)[0] + a = draw_values([self.a], point=point, size=size)[0] def _random(a, size=None): return stats.dirichlet.rvs(a, None if size == a.shape else size) @@ -545,7 +545,7 @@ def _random(self, n, p, size=None): return randnum.astype(original_dtype) def random(self, point=None, size=None): - n, p = draw_values([self.n, self.p], point=point) + n, p = draw_values([self.n, self.p], point=point, size=size) samples = generate_samples(self._random, n, p, dist_shape=self.shape, size=size) @@ -679,7 +679,7 @@ def __init__(self, nu, V, *args, **kwargs): np.nan) def random(self, point=None, size=None): - nu, V = draw_values([self.nu, self.V], point=point) + nu, V = draw_values([self.nu, self.V], point=point, size=size) size= 1 if size is None else size return generate_samples(stats.wishart.rvs, np.asscalar(nu), V, broadcast_shape=(size,)) @@ -1071,7 +1071,7 @@ def _random(self, n, eta, size=None): return C.transpose((1, 2, 0))[np.triu_indices(n, k=1)].T def random(self, point=None, size=None): - n, eta = draw_values([self.n, self.eta], point=point) + n, eta = draw_values([self.n, self.eta], point=point, size=size) size= 1 if size is None else size samples = generate_samples(self._random, n, eta, broadcast_shape=(size,)) @@ -1264,8 +1264,8 @@ def random(self, point=None, size=None): mu, colchol, rowchol = draw_values( [self.mu, self.colchol_cov, self.rowchol_cov], - point=point - ) + point=point, + size=size) standard_normal = np.random.standard_normal(size) return mu + np.matmul(rowchol, np.matmul(standard_normal, colchol.T)) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 9d450abda9..cf361c8d36 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -97,20 +97,6 @@ def test_random_sample_returns_nd_array(self): assert isinstance(mu, np.ndarray) assert isinstance(tau, np.ndarray) - def test_random_sample_returns_correctly(self): - # Based on what we discovered in #GH2909 - with pm.Model(): - a = pm.Uniform('a', lower=0, upper=1, shape=10) - b = pm.Binomial('b', n=1, p=a, shape=10) - array_of_uniform = a.random(size=10000).mean(axis=0) - array_of_binomial = b.random(size=10000).mean(axis=0) - npt.assert_allclose(array_of_uniform, [0.49886929, 0.49949713, 0.49946077, 0.49922606, 0.49927498, 0.50003914, - 0.49980687, 0.50180495, 0.500905, 0.50035121], rtol=1e-2, atol=0) - npt.assert_allclose(array_of_binomial, [0.7232, 0.131 , 0.9457, 0.8279, 0.2911, 0.8686, 0.57 , 0.9184, - 0.8177, 0.1625], rtol=1e-2, atol=0) - assert isinstance(array_of_binomial, np.ndarray) - assert isinstance(array_of_uniform, np.ndarray) - class BaseTestCases(object): class BaseTestCase(SeededTest): @@ -741,17 +727,17 @@ def test_lkj(self): for n in [2, 10, 50]: #pylint: disable=cell-var-from-loop shape = n*(n-1)//2 - + def ref_rand(size, eta): beta = eta - 1 + n/2 return (st.beta.rvs(size=(size, shape), a=beta, b=beta)-.5)*2 class TestedLKJCorr (pm.LKJCorr): - + def __init__(self, **kwargs): kwargs.pop('shape', None) super(TestedLKJCorr, self).__init__( - n=n, + n=n, **kwargs ) @@ -767,12 +753,12 @@ def ref_rand(size, w, mu, sd): pymc3_random(pm.NormalMixture, {'w': Simplex(2), 'mu': Domain([[.05, 2.5], [-5., 1.]], edges=(None, None)), - 'sd': Domain([[1, 1], [1.5, 2.]], edges=(None, None))}, + 'sd': Domain([[1, 1], [1.5, 2.]], edges=(None, None))}, size=1000, ref_rand=ref_rand) pymc3_random(pm.NormalMixture, {'w': Simplex(3), 'mu': Domain([[-5., 1., 2.5]], edges=(None, None)), - 'sd': Domain([[1.5, 2., 3.]], edges=(None, None))}, + 'sd': Domain([[1.5, 2., 3.]], edges=(None, None))}, size=1000, ref_rand=ref_rand)