diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 3ffa68df93..9e03eb6d61 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -26,6 +26,7 @@ from aesara.assert_op import Assert from aesara.tensor.random.basic import ( BetaRV, + WeibullRV, cauchy, exponential, gamma, @@ -33,9 +34,13 @@ halfcauchy, halfnormal, invgamma, + logistic, + lognormal, normal, pareto, + triangular, uniform, + vonmises, ) from aesara.tensor.random.op import RandomVariable from aesara.tensor.var import TensorVariable @@ -61,6 +66,7 @@ from pymc3.distributions.distribution import Continuous from pymc3.distributions.special import log_i0 from pymc3.math import invlogit, log1mexp, log1pexp, logdiffexp, logit +from pymc3.util import UNSET __all__ = [ "Uniform", @@ -106,8 +112,8 @@ class UnitContinuous(Continuous): """Base class for continuous distributions on [0,1]""" -class BoundedContinuous(Continuous): - """Base class for bounded continuous distributions""" +class CircularContinuous(Continuous): + """Base class for circular continuous distributions""" @logp_transform.register(PositiveContinuous) @@ -120,15 +126,39 @@ def unit_cont_transform(op): return transforms.logodds -@logp_transform.register(BoundedContinuous) -def bounded_cont_transform(op): - def transform_params(rv_var): - _, _, _, lower, upper = rv_var.owner.inputs - lower = at.as_tensor_variable(lower) if lower is not None else None - upper = at.as_tensor_variable(upper) if upper is not None else None - return lower, upper +@logp_transform.register(CircularContinuous) +def circ_cont_transform(op): + return transforms.circular + + +class BoundedContinuous(Continuous): + """Base class for bounded continuous distributions""" + + # Indices of the arguments that define the lower and upper bounds of the distribution + bound_args_indices = None + + def __new__(cls, *args, **kwargs): + transform = kwargs.get("transform", UNSET) + if transform is UNSET: + kwargs["transform"] = cls.default_transform() + return super().__new__(cls, *args, **kwargs) + + @classmethod + def default_transform(cls): + if cls.bound_args_indices is None: + raise ValueError( + f"Must specify bound_args_indices for {cls.__name__} bounded distribution" + ) + + def transform_params(rv_var): + _, _, _, *args = rv_var.owner.inputs + lower = args[cls.bound_args_indices[0]] + upper = args[cls.bound_args_indices[1]] + lower = at.as_tensor_variable(lower) if lower is not None else None + upper = at.as_tensor_variable(upper) if upper is not None else None + return lower, upper - return transforms.interval(transform_params) + return transforms.interval(transform_params) def assert_negative_support(var, label, distname, value=-1e-6): @@ -217,6 +247,7 @@ class Uniform(BoundedContinuous): Upper limit. """ rv_op = uniform + bound_args_indices = (0, 1) # Lower, Upper @classmethod def dist(cls, lower=0, upper=1, **kwargs): @@ -1716,50 +1747,24 @@ class Lognormal(PositiveContinuous): x = pm.Lognormal('x', mu=2, tau=1/100) """ - def __init__(self, mu=0, sigma=None, tau=None, sd=None, *args, **kwargs): - super().__init__(*args, **kwargs) + rv_op = lognormal + + @classmethod + def dist(cls, mu=0, sigma=None, tau=None, sd=None, *args, **kwargs): if sd is not None: sigma = sd tau, sigma = get_tau_sigma(tau=tau, sigma=sigma) - self.mu = mu = at.as_tensor_variable(floatX(mu)) - self.tau = tau = at.as_tensor_variable(tau) - self.sigma = self.sd = sigma = at.as_tensor_variable(sigma) - - self.mean = at.exp(self.mu + 1.0 / (2 * self.tau)) - self.median = at.exp(self.mu) - self.mode = at.exp(self.mu - 1.0 / self.tau) - self.variance = (at.exp(1.0 / self.tau) - 1) * at.exp(2 * self.mu + 1.0 / self.tau) - - assert_negative_support(tau, "tau", "Lognormal") - assert_negative_support(sigma, "sigma", "Lognormal") - - def _random(self, mu, tau, size=None): - samples = np.random.normal(size=size) - return np.exp(mu + (tau ** -0.5) * samples) - - def random(self, point=None, size=None): - """ - Draw random values from Lognormal distribution. + mu = at.as_tensor_variable(floatX(mu)) + sigma = at.as_tensor_variable(floatX(sigma)) - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). + assert_negative_support(tau, "tau", "LogNormal") + assert_negative_support(sigma, "sigma", "LogNormal") - Returns - ------- - array - """ - # mu, tau = draw_values([self.mu, self.tau], point=point, size=size) - # return generate_samples(self._random, mu, tau, dist_shape=self.shape, size=size) + return super().dist([mu, sigma], *args, **kwargs) - def logp(self, value): + def logp(value, mu, sigma): """ Calculate log-probability of Lognormal distribution at specified value. @@ -1773,8 +1778,7 @@ def logp(self, value): ------- TensorVariable """ - mu = self.mu - tau = self.tau + tau, sigma = get_tau_sigma(tau=None, sigma=sigma) return bound( -0.5 * tau * (at.log(value) - mu) ** 2 + 0.5 * at.log(tau / (2.0 * np.pi)) @@ -1782,10 +1786,7 @@ def logp(self, value): tau > 0, ) - def _distr_parameters_for_repr(self): - return ["mu", "tau"] - - def logcdf(self, value): + def logcdf(value, mu, sigma): """ Compute the log of the cumulative distribution function for Lognormal distribution at the specified value. @@ -1800,14 +1801,11 @@ def logcdf(self, value): ------- TensorVariable """ - mu = self.mu - sigma = self.sigma - tau = self.tau return bound( normal_lcdf(mu, sigma, at.log(value)), 0 < value, - 0 < tau, + 0 < sigma, ) @@ -2646,6 +2644,18 @@ def __init__(self, nu, *args, **kwargs): super().__init__(alpha=nu / 2.0, beta=0.5, *args, **kwargs) +# TODO: Remove this once logpt for multiplication is working! +class WeibullBetaRV(WeibullRV): + ndims_params = [0, 0] + + @classmethod + def rng_fn(cls, rng, alpha, beta, size): + return beta * rng.weibull(alpha, size=size) + + +weibull_beta = WeibullBetaRV() + + class Weibull(PositiveContinuous): r""" Weibull log-likelihood. @@ -2691,45 +2701,19 @@ class Weibull(PositiveContinuous): Scale parameter (beta > 0). """ - def __init__(self, alpha, beta, *args, **kwargs): - super().__init__(*args, **kwargs) - self.alpha = alpha = at.as_tensor_variable(floatX(alpha)) - self.beta = beta = at.as_tensor_variable(floatX(beta)) - self.mean = beta * at.exp(gammaln(1 + 1.0 / alpha)) - self.median = beta * at.exp(gammaln(at.log(2))) ** (1.0 / alpha) - self.variance = beta ** 2 * at.exp(gammaln(1 + 2.0 / alpha)) - self.mean ** 2 - self.mode = at.switch( - alpha >= 1, beta * ((alpha - 1) / alpha) ** (1 / alpha), 0 - ) # Reference: https://en.wikipedia.org/wiki/Weibull_distribution + rv_op = weibull_beta + + @classmethod + def dist(cls, alpha, beta, *args, **kwargs): + alpha = at.as_tensor_variable(floatX(alpha)) + beta = at.as_tensor_variable(floatX(beta)) assert_negative_support(alpha, "alpha", "Weibull") assert_negative_support(beta, "beta", "Weibull") - def random(self, point=None, size=None): - """ - Draw random values from Weibull distribution. + return super().dist([alpha, beta], *args, **kwargs) - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). - - Returns - ------- - array - """ - # alpha, beta = draw_values([self.alpha, self.beta], point=point, size=size) - # - # def _random(a, b, size=None): - # return b * (-np.log(np.random.uniform(size=size))) ** (1 / a) - # - # return generate_samples(_random, alpha, beta, dist_shape=self.shape, size=size) - - def logp(self, value): + def logp(value, alpha, beta): """ Calculate log-probability of Weibull distribution at specified value. @@ -2743,8 +2727,6 @@ def logp(self, value): ------- TensorVariable """ - alpha = self.alpha - beta = self.beta return bound( at.log(alpha) - at.log(beta) @@ -2755,7 +2737,7 @@ def logp(self, value): beta > 0, ) - def logcdf(self, value): + def logcdf(value, alpha, beta): r""" Compute the log of the cumulative distribution function for Weibull distribution at the specified value. @@ -2770,8 +2752,6 @@ def logcdf(self, value): ------- TensorVariable """ - alpha = self.alpha - beta = self.beta a = (value / beta) ** alpha return bound( log1mexp(a), @@ -3099,7 +3079,7 @@ def _distr_parameters_for_repr(self): return ["mu", "sigma", "nu"] -class VonMises(Continuous): +class VonMises(CircularContinuous): r""" Univariate VonMises log-likelihood. @@ -3144,38 +3124,16 @@ class VonMises(Continuous): Concentration (\frac{1}{kappa} is analogous to \sigma^2). """ - def __init__(self, mu=0.0, kappa=None, transform="circular", *args, **kwargs): - if transform == "circular": - transform = transforms.Circular() - super().__init__(transform=transform, *args, **kwargs) - self.mean = self.median = self.mode = self.mu = mu = at.as_tensor_variable(floatX(mu)) - self.kappa = kappa = at.as_tensor_variable(floatX(kappa)) + rv_op = vonmises + @classmethod + def dist(cls, mu=0.0, kappa=None, *args, **kwargs): + mu = at.as_tensor_variable(floatX(mu)) + kappa = at.as_tensor_variable(floatX(kappa)) assert_negative_support(kappa, "kappa", "VonMises") + return super().dist([mu, kappa], *args, **kwargs) - def random(self, point=None, size=None): - """ - Draw random values from VonMises distribution. - - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). - - Returns - ------- - array - """ - # mu, kappa = draw_values([self.mu, self.kappa], point=point, size=size) - # return generate_samples( - # stats.vonmises.rvs, loc=mu, kappa=kappa, dist_shape=self.shape, size=size - # ) - - def logp(self, value): + def logp(value, mu, kappa): """ Calculate log-probability of VonMises distribution at specified value. @@ -3189,8 +3147,6 @@ def logp(self, value): ------- TensorVariable """ - mu = self.mu - kappa = self.kappa return bound( kappa * at.cos(mu - value) - (at.log(2 * np.pi) + log_i0(kappa)), kappa > 0, @@ -3198,9 +3154,6 @@ def logp(self, value): value <= np.pi, ) - def _distr_parameters_for_repr(self): - return ["mu", "kappa"] - class SkewNormal(Continuous): r""" @@ -3388,45 +3341,18 @@ class Triangular(BoundedContinuous): Upper limit. """ - def __init__(self, lower=0, upper=1, c=0.5, *args, **kwargs): - self.median = self.mean = self.c = c = at.as_tensor_variable(floatX(c)) - self.lower = lower = at.as_tensor_variable(floatX(lower)) - self.upper = upper = at.as_tensor_variable(floatX(upper)) - - super().__init__(lower=lower, upper=upper, *args, **kwargs) - - def random(self, point=None, size=None): - """ - Draw random values from Triangular distribution. + rv_op = triangular + bound_args_indices = (0, 2) # lower, upper - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). - - Returns - ------- - array - """ - # c, lower, upper = draw_values([self.c, self.lower, self.upper], point=point, size=size) - # return generate_samples( - # self._random, c=c, lower=lower, upper=upper, size=size, dist_shape=self.shape - # ) + @classmethod + def dist(cls, lower=0, upper=1, c=0.5, *args, **kwargs): + lower = at.as_tensor_variable(floatX(lower)) + upper = at.as_tensor_variable(floatX(upper)) + c = at.as_tensor_variable(floatX(c)) - def _random(self, c, lower, upper, size): - """Wrapper around stats.triang.rvs that converts Triangular's - parametrization to scipy.triang. All parameter arrays should have - been broadcasted properly by generate_samples at this point and size is - the scipy.rvs representation. - """ - scale = upper - lower - return stats.triang.rvs(c=(c - lower) / scale, loc=lower, scale=scale, size=size) + return super().dist([lower, c, upper], *args, **kwargs) - def logp(self, value): + def logp(value, lower, c, upper): """ Calculate log-probability of Triangular distribution at specified value. @@ -3440,9 +3366,6 @@ def logp(self, value): ------- TensorVariable """ - c = self.c - lower = self.lower - upper = self.upper return bound( at.switch( at.lt(value, c), @@ -3451,9 +3374,11 @@ def logp(self, value): ), lower <= value, value <= upper, + lower <= c, + c <= upper, ) - def logcdf(self, value): + def logcdf(value, lower, c, upper): """ Compute the log of the cumulative distribution function for Triangular distribution at the specified value. @@ -3468,9 +3393,6 @@ def logcdf(self, value): ------- TensorVariable """ - c = self.c - lower = self.lower - upper = self.upper return bound( at.switch( at.le(value, lower), @@ -3485,7 +3407,8 @@ def logcdf(self, value): ), ), ), - lower <= upper, + lower <= c, + c <= upper, ) @@ -3810,39 +3733,15 @@ class Logistic(Continuous): Scale (s > 0). """ - def __init__(self, mu=0.0, s=1.0, *args, **kwargs): - super().__init__(*args, **kwargs) - - self.mu = at.as_tensor_variable(floatX(mu)) - self.s = at.as_tensor_variable(floatX(s)) - - self.mean = self.mode = mu - self.variance = s ** 2 * np.pi ** 2 / 3.0 + rv_op = logistic - def random(self, point=None, size=None): - """ - Draw random values from Logistic distribution. - - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). - - Returns - ------- - array - """ - # mu, s = draw_values([self.mu, self.s], point=point, size=size) - # - # return generate_samples( - # stats.logistic.rvs, loc=mu, scale=s, dist_shape=self.shape, size=size - # ) + @classmethod + def dist(cls, mu=0.0, s=1.0, *args, **kwargs): + mu = at.as_tensor_variable(floatX(mu)) + s = at.as_tensor_variable(floatX(s)) + return super().dist([mu, s], *args, **kwargs) - def logp(self, value): + def logp(value, mu, s): """ Calculate log-probability of Logistic distribution at specified value. @@ -3856,15 +3755,13 @@ def logp(self, value): ------- TensorVariable """ - mu = self.mu - s = self.s return bound( -(value - mu) / s - at.log(s) - 2 * at.log1p(at.exp(-(value - mu) / s)), s > 0, ) - def logcdf(self, value): + def logcdf(value, mu, s): r""" Compute the log of the cumulative distribution function for Logistic distribution at the specified value. @@ -3879,8 +3776,7 @@ def logcdf(self, value): ------- TensorVariable """ - mu = self.mu - s = self.s + return bound( -log1pexp(-(value - mu) / s), 0 < s, diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index 633539d178..03672ac41e 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -19,8 +19,11 @@ from aesara.tensor.random.basic import ( RandomVariable, bernoulli, + betabinom, binomial, categorical, + geometric, + hypergeometric, nbinom, poisson, ) @@ -39,7 +42,7 @@ normal_lcdf, ) from pymc3.distributions.distribution import Discrete -from pymc3.math import log1mexp, logaddexp, logsumexp, sigmoid, tround +from pymc3.math import log1mexp, logaddexp, logsumexp, sigmoid __all__ = [ "Binomial", @@ -225,58 +228,16 @@ def BetaBinom(a, b, n, x): beta > 0. """ - def __init__(self, alpha, beta, n, *args, **kwargs): - super().__init__(*args, **kwargs) - self.alpha = alpha = at.as_tensor_variable(floatX(alpha)) - self.beta = beta = at.as_tensor_variable(floatX(beta)) - self.n = n = at.as_tensor_variable(intX(n)) - self.mode = at.cast(tround(alpha / (alpha + beta)), "int8") - - def _random(self, alpha, beta, n, size=None): - size = size or () - p = stats.beta.rvs(a=alpha, b=beta, size=size).flatten() - # Sometimes scipy.beta returns nan. Ugh. - while np.any(np.isnan(p)): - i = np.isnan(p) - p[i] = stats.beta.rvs(a=alpha, b=beta, size=np.sum(i)) - # Sigh... - _n, _p, _size = np.atleast_1d(n).flatten(), p.flatten(), p.shape[0] - - quotient, remainder = divmod(_p.shape[0], _n.shape[0]) - if remainder != 0: - raise TypeError( - "n has a bad size! Was cast to {}, must evenly divide {}".format( - _n.shape[0], _p.shape[0] - ) - ) - if quotient != 1: - _n = np.tile(_n, quotient) - samples = np.reshape(stats.binom.rvs(n=_n, p=_p, size=_size), size) - return samples - - def random(self, point=None, size=None): - r""" - Draw random values from BetaBinomial distribution. - - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). + rv_op = betabinom - Returns - ------- - array - """ - # alpha, beta, n = 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 - # ) + @classmethod + def dist(cls, alpha, beta, n, *args, **kwargs): + alpha = at.as_tensor_variable(floatX(alpha)) + beta = at.as_tensor_variable(floatX(beta)) + n = at.as_tensor_variable(intX(n)) + return super().dist([n, alpha, beta], **kwargs) - def logp(self, value): + def logp(value, n, alpha, beta): r""" Calculate log-probability of BetaBinomial distribution at specified value. @@ -290,9 +251,6 @@ def logp(self, value): ------- TensorVariable """ - alpha = self.alpha - beta = self.beta - n = self.n return bound( binomln(n, value) + betaln(value + alpha, n - value + beta) - betaln(alpha, beta), value >= 0, @@ -301,7 +259,7 @@ def logp(self, value): beta > 0, ) - def logcdf(self, value): + def logcdf(value, n, alpha, beta): """ Compute the log of the cumulative distribution function for BetaBinomial distribution at the specified value. @@ -321,15 +279,15 @@ def logcdf(self, value): f"BetaBinomial.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object." ) - alpha = self.alpha - beta = self.beta - n = self.n safe_lower = at.switch(at.lt(value, 0), value, 0) return bound( at.switch( at.lt(value, n), - logsumexp(self.logp(at.arange(safe_lower, value + 1)), keepdims=False), + logsumexp( + BetaBinomial.logp(at.arange(safe_lower, value + 1), n, alpha, beta), + keepdims=False, + ), 0, ), 0 <= value, @@ -833,32 +791,14 @@ class Geometric(Discrete): Probability of success on an individual trial (0 < p <= 1). """ - def __init__(self, p, *args, **kwargs): - super().__init__(*args, **kwargs) - self.p = p = at.as_tensor_variable(floatX(p)) - self.mode = 1 - - def random(self, point=None, size=None): - r""" - Draw random values from Geometric distribution. - - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). + rv_op = geometric - Returns - ------- - array - """ - # p = draw_values([self.p], point=point, size=size)[0] - # return generate_samples(np.random.geometric, p, dist_shape=self.shape, size=size) + @classmethod + def dist(cls, p, *args, **kwargs): + p = at.as_tensor_variable(floatX(p)) + return super().dist([p], *args, **kwargs) - def logp(self, value): + def logp(value, p): r""" Calculate log-probability of Geometric distribution at specified value. @@ -872,10 +812,14 @@ def logp(self, value): ------- TensorVariable """ - p = self.p - return bound(at.log(p) + logpow(1 - p, value - 1), 0 <= p, p <= 1, value >= 1) + return bound( + at.log(p) + logpow(1 - p, value - 1), + 0 <= p, + p <= 1, + value >= 1, + ) - def logcdf(self, value): + def logcdf(value, p): """ Compute the log of the cumulative distribution function for Geometric distribution at the specified value. @@ -890,7 +834,6 @@ def logcdf(self, value): ------- TensorVariable """ - p = self.p return bound( log1mexp(-at.log1p(-p) * value), @@ -947,43 +890,16 @@ class HyperGeometric(Discrete): Number of samples drawn from the population """ - def __init__(self, N, k, n, *args, **kwargs): - super().__init__(*args, **kwargs) - self.N = intX(N) - self.k = intX(k) - self.n = intX(n) - self.mode = intX(at.floor((n + 1) * (k + 1) / (N + 2))) - - def random(self, point=None, size=None): - r""" - Draw random values from HyperGeometric distribution. - - Parameters - ---------- - point : dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size : int, optional - Desired size of random sample (returns one sample if not - specified). - - Returns - ------- - array - """ - - # N, k, n = draw_values([self.N, self.k, self.n], point=point, size=size) - # return generate_samples(self._random, N, k, n, dist_shape=self.shape, size=size) + rv_op = hypergeometric - def _random(self, M, n, N, size=None): - r"""Wrapper around scipy stat's hypergeom.rvs""" - try: - samples = stats.hypergeom.rvs(M=M, n=n, N=N, size=size) - return samples - except ValueError: - raise ValueError("Domain error in arguments") + @classmethod + def dist(cls, N, k, n, *args, **kwargs): + good = at.as_tensor_variable(intX(k)) + bad = at.as_tensor_variable(intX(N - k)) + n = at.as_tensor_variable(intX(n)) + return super().dist([good, bad, n], *args, **kwargs) - def logp(self, value): + def logp(value, good, bad, n): r""" Calculate log-probability of HyperGeometric distribution at specified value. @@ -997,11 +913,8 @@ def logp(self, value): ------- TensorVariable """ - N = self.N - k = self.k - n = self.n - tot, good = N, k - bad = tot - good + + tot = good + bad result = ( betaln(good + 1, 1) + betaln(bad + 1, 1) @@ -1011,11 +924,11 @@ def logp(self, value): - betaln(tot + 1, 1) ) # value in [max(0, n - N + k), min(k, n)] - lower = at.switch(at.gt(n - N + k, 0), n - N + k, 0) - upper = at.switch(at.lt(k, n), k, n) + lower = at.switch(at.gt(n - tot + good, 0), n - tot + good, 0) + upper = at.switch(at.lt(good, n), good, n) return bound(result, lower <= value, value <= upper) - def logcdf(self, value): + def logcdf(value, good, bad, n): """ Compute the log of the cumulative distribution function for HyperGeometric distribution at the specified value. @@ -1035,23 +948,24 @@ def logcdf(self, value): f"HyperGeometric.logcdf expects a scalar value but received a {np.ndim(value)}-dimensional object." ) + N = good + bad # TODO: Use lower upper in locgdf for smarter logsumexp? - N = self.N - n = self.n - k = self.k safe_lower = at.switch(at.lt(value, 0), value, 0) return bound( at.switch( at.lt(value, n), - logsumexp(self.logp(at.arange(safe_lower, value + 1)), keepdims=False), + logsumexp( + HyperGeometric.logp(at.arange(safe_lower, value + 1), good, bad, n), + keepdims=False, + ), 0, ), 0 <= value, 0 < N, - 0 <= k, + 0 <= good, 0 <= n, - k <= N, + good <= N, n <= N, ) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index 59d2e185a1..a6f9922371 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -109,13 +109,6 @@ def logcdf(op, var, rvs_to_values, *dist_params, **kwargs): value_var = rvs_to_values.get(var, var) return class_logcdf(value_var, *dist_params, **kwargs) - # class_transform = clsdict.get("transform") - # if class_transform: - # - # @logp_transform.register(rv_type) - # def transform(op, *args, **kwargs): - # return class_transform(*args, **kwargs) - # Register the Aesara `RandomVariable` type as a subclass of this # `Distribution` type. new_cls.register(rv_type) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index c6869ec110..397d2ba6af 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -883,7 +883,6 @@ def test_uniform(self): assert logpt(invalid_dist, np.array(0.5)).eval() == -np.inf assert logcdf(invalid_dist, np.array(2.0)).eval() == -np.inf - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_triangular(self): self.check_logp( Triangular, @@ -898,12 +897,20 @@ def test_triangular(self): lambda value, c, lower, upper: sp.triang.logcdf(value, c - lower, lower, upper - lower), skip_paramdomain_outside_edge_test=True, ) - # Custom logp check for invalid value - valid_dist = Triangular.dist(lower=0, upper=1, c=2.0) - assert np.all(logpt(valid_dist, np.array([1.9, 2.0, 2.1])).tag.test_value == -np.inf) + + # Custom logp/logcdf check for values outside of domain + valid_dist = Triangular.dist(lower=0, upper=1, c=0.9, size=2) + with aesara.config.change_flags(mode=Mode("py")): + assert np.all(logpt(valid_dist, np.array([-1, 2])).eval() == -np.inf) + assert np.all(logcdf(valid_dist, np.array([-1, 2])).eval() == [-np.inf, 0]) # Custom logp / logcdf check for invalid parameters - invalid_dist = Triangular.dist(lower=1, upper=0, c=2.0) + invalid_dist = Triangular.dist(lower=1, upper=0, c=0.1) + with aesara.config.change_flags(mode=Mode("py")): + assert logpt(invalid_dist, 0.5).eval() == -np.inf + assert logcdf(invalid_dist, 2).eval() == -np.inf + + invalid_dist = Triangular.dist(lower=0, upper=1, c=2.0) with aesara.config.change_flags(mode=Mode("py")): assert logpt(invalid_dist, 0.5).eval() == -np.inf assert logcdf(invalid_dist, 2).eval() == -np.inf @@ -1121,7 +1128,6 @@ def test_exponential(self): lambda value, lam: sp.expon.logcdf(value, 0, 1 / lam), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_geometric(self): self.check_logp( Geometric, @@ -1141,7 +1147,6 @@ def test_geometric(self): {"p": Unit}, ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_hypergeometric(self): def modified_scipy_hypergeom_logpmf(value, N, k, n): # Convert nan to -np.inf @@ -1263,7 +1268,6 @@ def test_laplace_asymmetric(self): laplace_asymmetric_logpdf, ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_lognormal(self): self.check_logp( Lognormal, @@ -1271,12 +1275,26 @@ def test_lognormal(self): {"mu": R, "tau": Rplusbig}, lambda value, mu, tau: floatX(sp.lognorm.logpdf(value, tau ** -0.5, 0, np.exp(mu))), ) + self.check_logp( + Lognormal, + Rplus, + {"mu": R, "sigma": Rplusbig}, + lambda value, mu, sigma: floatX(sp.lognorm.logpdf(value, sigma, 0, np.exp(mu))), + n_samples=5, # Just testing alternative parametrization + ) self.check_logcdf( Lognormal, Rplus, {"mu": R, "tau": Rplusbig}, lambda value, mu, tau: sp.lognorm.logcdf(value, tau ** -0.5, 0, np.exp(mu)), ) + self.check_logcdf( + Lognormal, + Rplus, + {"mu": R, "sigma": Rplusbig}, + lambda value, mu, sigma: sp.lognorm.logcdf(value, sigma, 0, np.exp(mu)), + n_samples=5, # Just testing alternative parametrization + ) @pytest.mark.xfail(reason="Distribution not refactored yet") def test_t(self): @@ -1413,7 +1431,10 @@ def test_pareto(self): lambda value, alpha, m: sp.pareto.logcdf(value, alpha, scale=m), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") + @pytest.mark.xfail( + condition=(aesara.config.floatX == "float32"), + reason="Fails on float32 due to numerical issues", + ) def test_weibull_logp(self): self.check_logp( Weibull, @@ -1422,7 +1443,6 @@ def test_weibull_logp(self): lambda value, alpha, beta: sp.exponweib.logpdf(value, 1, alpha, scale=beta), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") @pytest.mark.xfail( condition=(aesara.config.floatX == "float32"), reason="Fails on float32 due to inf issues", @@ -1476,8 +1496,7 @@ def test_binomial(self): n_samples=10, ) - # Too lazy to propagate decimal parameter through the whole chain of deps - @pytest.mark.xfail(reason="Distribution not refactored yet") + @pytest.mark.xfail(reason="checkd tests has not been refactored") @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") def test_beta_binomial_distribution(self): self.checkd( @@ -1486,7 +1505,6 @@ def test_beta_binomial_distribution(self): {"alpha": Rplus, "beta": Rplus, "n": NatSmall}, ) - @pytest.mark.xfail(reason="Distribution not refactored yet") @pytest.mark.skipif( condition=(SCIPY_VERSION < parse("1.4.0")), reason="betabinom is new in Scipy 1.4.0" ) @@ -1498,7 +1516,6 @@ def test_beta_binomial_logp(self): lambda value, alpha, beta, n: sp.betabinom.logpmf(value, a=alpha, b=beta, n=n), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") @pytest.mark.skipif( condition=(SCIPY_VERSION < parse("1.4.0")), reason="betabinom is new in Scipy 1.4.0" @@ -1511,7 +1528,6 @@ def test_beta_binomial_logcdf(self): lambda value, alpha, beta, n: sp.betabinom.logcdf(value, a=alpha, b=beta, n=n), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_beta_binomial_selfconsistency(self): self.check_selfconsistency_discrete_logcdf( BetaBinomial, @@ -2381,7 +2397,6 @@ def test_ex_gaussian_cdf_outside_edges(self): ) @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_vonmises(self): self.check_logp( VonMises, @@ -2402,7 +2417,6 @@ def gumbellcdf(value, mu, beta): self.check_logcdf(Gumbel, R, {"mu": R, "beta": Rplusbig}, gumbellcdf) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_logistic(self): self.check_logp( Logistic, diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 03190f3bcd..1a7aa185a4 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -13,7 +13,6 @@ # limitations under the License. import functools import itertools -import sys from contextlib import ExitStack as does_not_raise from typing import Callable, List, Optional @@ -26,6 +25,8 @@ import scipy.stats as st from numpy.testing import assert_almost_equal, assert_array_almost_equal +from packaging.version import parse +from scipy import __version__ as scipy_version from scipy.special import expit import pymc3 as pm @@ -36,7 +37,7 @@ from pymc3.distributions.multivariate import quaddist_matrix from pymc3.distributions.shape_utils import to_tuple from pymc3.exceptions import ShapeError -from pymc3.tests.helpers import SeededTest +from pymc3.tests.helpers import SeededTest, select_by_precision from pymc3.tests.test_distributions import ( Domain, I, @@ -49,15 +50,14 @@ RealMatrix, Rplus, Rplusbig, - Rplusdunif, - Runif, Simplex, - Unit, Vector, build_model, product, ) +SCIPY_VERSION = parse(scipy_version) + def pymc3_random( dist, @@ -267,12 +267,6 @@ class TestSkewNormal(BaseTestCases.BaseTestCase): params = {"mu": 0.0, "sigma": 1.0, "alpha": 5.0} -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestTriangular(BaseTestCases.BaseTestCase): - distribution = pm.Triangular - params = {"c": 0.5, "lower": 0.0, "upper": 1.0} - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") class TestWald(BaseTestCases.BaseTestCase): distribution = pm.Wald @@ -297,12 +291,6 @@ class TestAsymmetricLaplace(BaseTestCases.BaseTestCase): params = {"kappa": 1.0, "b": 1.0, "mu": 0.0} -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestLognormal(BaseTestCases.BaseTestCase): - distribution = pm.Lognormal - params = {"mu": 1.0, "tau": 1.0} - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") class TestStudentT(BaseTestCases.BaseTestCase): distribution = pm.StudentT @@ -315,42 +303,18 @@ class TestChiSquared(BaseTestCases.BaseTestCase): params = {"nu": 2.0} -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestWeibull(BaseTestCases.BaseTestCase): - distribution = pm.Weibull - params = {"alpha": 1.0, "beta": 1.0} - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") class TestExGaussian(BaseTestCases.BaseTestCase): distribution = pm.ExGaussian params = {"mu": 0.0, "sigma": 1.0, "nu": 1.0} -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestVonMises(BaseTestCases.BaseTestCase): - distribution = pm.VonMises - params = {"mu": 0.0, "kappa": 1.0} - - -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestLogistic(BaseTestCases.BaseTestCase): - distribution = pm.Logistic - params = {"mu": 0.0, "s": 1.0} - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") class TestLogitNormal(BaseTestCases.BaseTestCase): distribution = pm.LogitNormal params = {"mu": 0.0, "sigma": 1.0} -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestBetaBinomial(BaseTestCases.BaseTestCase): - distribution = pm.BetaBinomial - params = {"n": 5, "alpha": 1.0, "beta": 1.0} - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") class TestConstant(BaseTestCases.BaseTestCase): distribution = pm.Constant @@ -381,18 +345,6 @@ class TestDiscreteUniform(BaseTestCases.BaseTestCase): params = {"lower": 0.0, "upper": 10.0} -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestGeometric(BaseTestCases.BaseTestCase): - distribution = pm.Geometric - params = {"p": 0.5} - - -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") -class TestHyperGeometric(BaseTestCases.BaseTestCase): - distribution = pm.HyperGeometric - params = {"N": 50, "k": 25, "n": 10} - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") class TestMoyal(BaseTestCases.BaseTestCase): distribution = pm.Moyal @@ -407,7 +359,7 @@ class BaseTestDistribution(SeededTest): expected_rv_op_params = dict() tests_to_run = [] size = 15 - decimal = 6 + decimal = select_by_precision(float64=6, float32=3) sizes_to_check: Optional[List] = None sizes_expected: Optional[List] = None @@ -454,15 +406,15 @@ def check_rv_size(self): sizes_expected = self.sizes_expected or [(), (), (1,), (1,), (5,), (4, 5), (2, 4, 2)] for size, expected in zip(sizes_to_check, sizes_expected): actual = change_rv_size(self.pymc_rv, size).eval().shape - assert actual == expected + assert actual == expected, f"size={size}, expected={expected}, actual={actual}" # test negative sizes raise for size in [-2, (3, -2)]: with pytest.raises(ValueError): change_rv_size(self.pymc_rv, size).eval() - # test multi-parameters sampling for univariate distributions - if self.pymc_dist.rv_op.ndim_supp == 0: + # test multi-parameters sampling for univariate distributions (with univariate inputs) + if self.pymc_dist.rv_op.ndim_supp == 0 and sum(self.pymc_dist.rv_op.ndims_params) == 0: params = { k: p * np.ones(self.repeated_params_shape) for k, p in self.pymc_dist_params.items() } @@ -512,8 +464,8 @@ def seeded_discrete_weibul_rng_fn(self): reference_dist = seeded_discrete_weibul_rng_fn tests_to_run = [ "check_pymc_params_match_rv_op", - "check_rv_size", "check_pymc_draws_match_reference", + "check_rv_size", ] @@ -522,11 +474,9 @@ class TestGumbel(BaseTestDistribution): pymc_dist_params = {"mu": 1.5, "beta": 3.0} expected_rv_op_params = {"mu": 1.5, "beta": 3.0} reference_dist_params = {"loc": 1.5, "scale": 3.0} - size = 15 reference_dist = seeded_scipy_distribution_builder("gumbel_r") tests_to_run = [ "check_pymc_params_match_rv_op", - "check_rv_size", "check_pymc_draws_match_reference", ] @@ -540,16 +490,16 @@ class TestNormal(BaseTestDistribution): reference_dist = seeded_numpy_distribution_builder("normal") tests_to_run = [ "check_pymc_params_match_rv_op", - "check_rv_size", "check_pymc_draws_match_reference", + "check_rv_size", ] class TestNormalTau(BaseTestDistribution): pymc_dist = pm.Normal tau, sigma = get_tau_sigma(tau=25.0) - pymc_dist_params = {"mu": 1.0, "sigma": sigma} - expected_rv_op_params = {"mu": 1.0, "sigma": 0.2} + pymc_dist_params = {"mu": 1.0, "tau": tau} + expected_rv_op_params = {"mu": 1.0, "sigma": sigma} tests_to_run = ["check_pymc_params_match_rv_op"] @@ -571,14 +521,19 @@ class TestHalfNormal(BaseTestDistribution): pymc_dist = pm.HalfNormal pymc_dist_params = {"sigma": 10.0} expected_rv_op_params = {"mean": 0, "sigma": 10.0} - tests_to_run = ["check_pymc_params_match_rv_op"] + reference_dist_params = {"loc": 0, "scale": 10.0} + reference_dist = seeded_scipy_distribution_builder("halfnorm") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] class TestHalfNormalTau(BaseTestDistribution): pymc_dist = pm.Normal tau, sigma = get_tau_sigma(tau=25.0) - pymc_dist_params = {"sigma": sigma} - expected_rv_op_params = {"mu": 0.0, "sigma": 0.2} + pymc_dist_params = {"tau": tau} + expected_rv_op_params = {"mu": 0.0, "sigma": sigma} tests_to_run = ["check_pymc_params_match_rv_op"] @@ -600,8 +555,8 @@ class TestBeta(BaseTestDistribution): ) tests_to_run = [ "check_pymc_params_match_rv_op", - "check_rv_size", "check_pymc_draws_match_reference", + "check_rv_size", ] @@ -618,29 +573,49 @@ class TestBetaMuSigma(BaseTestDistribution): class TestExponential(BaseTestDistribution): pymc_dist = pm.Exponential pymc_dist_params = {"lam": 10.0} - expected_rv_op_params = {"lam": 1.0 / pymc_dist_params["lam"]} - tests_to_run = ["check_pymc_params_match_rv_op"] + expected_rv_op_params = {"mu": 1.0 / pymc_dist_params["lam"]} + reference_dist_params = {"scale": 1.0 / pymc_dist_params["lam"]} + reference_dist = seeded_numpy_distribution_builder("exponential") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] class TestCauchy(BaseTestDistribution): pymc_dist = pm.Cauchy pymc_dist_params = {"alpha": 2.0, "beta": 5.0} expected_rv_op_params = {"alpha": 2.0, "beta": 5.0} - tests_to_run = ["check_pymc_params_match_rv_op"] + reference_dist_params = {"loc": 2.0, "scale": 5.0} + reference_dist = seeded_scipy_distribution_builder("cauchy") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] -class TestHalfCauchyn(BaseTestDistribution): +class TestHalfCauchy(BaseTestDistribution): pymc_dist = pm.HalfCauchy pymc_dist_params = {"beta": 5.0} expected_rv_op_params = {"alpha": 0.0, "beta": 5.0} - tests_to_run = ["check_pymc_params_match_rv_op"] + reference_dist_params = {"loc": 0.0, "scale": 5.0} + reference_dist = seeded_scipy_distribution_builder("halfcauchy") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] class TestGamma(BaseTestDistribution): pymc_dist = pm.Gamma pymc_dist_params = {"alpha": 2.0, "beta": 5.0} expected_rv_op_params = {"alpha": 2.0, "beta": 1 / 5.0} - tests_to_run = ["check_pymc_params_match_rv_op"] + reference_dist_params = {"shape": 2.0, "scale": 1 / 5.0} + reference_dist = seeded_numpy_distribution_builder("gamma") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] class TestGammaMuSigma(BaseTestDistribution): @@ -657,7 +632,12 @@ class TestInverseGamma(BaseTestDistribution): pymc_dist = pm.InverseGamma pymc_dist_params = {"alpha": 2.0, "beta": 5.0} expected_rv_op_params = {"alpha": 2.0, "beta": 5.0} - tests_to_run = ["check_pymc_params_match_rv_op"] + reference_dist_params = {"a": 2.0, "scale": 5.0} + reference_dist = seeded_scipy_distribution_builder("invgamma") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] class TestInverseGammaMuSigma(BaseTestDistribution): @@ -704,7 +684,12 @@ class TestBernoulli(BaseTestDistribution): pymc_dist = pm.Bernoulli pymc_dist_params = {"p": 0.33} expected_rv_op_params = {"p": 0.33} - tests_to_run = ["check_pymc_params_match_rv_op"] + reference_dist_params = {"p": 0.33} + reference_dist = seeded_scipy_distribution_builder("bernoulli") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] @pytest.mark.skip("Still not implemented") @@ -734,7 +719,16 @@ class TestMvNormal(BaseTestDistribution): } sizes_to_check = [None, (1), (2, 3)] sizes_expected = [(2,), (1, 2), (2, 3, 2)] - tests_to_run = ["check_pymc_params_match_rv_op", "check_rv_size"] + reference_dist_params = { + "mean": np.array([1.0, 2.0]), + "cov": np.array([[2.0, 0.0], [0.0, 3.5]]), + } + reference_dist = seeded_numpy_distribution_builder("multivariate_normal") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] class TestMvNormalChol(BaseTestDistribution): @@ -767,7 +761,15 @@ class TestDirichlet(BaseTestDistribution): pymc_dist = pm.Dirichlet pymc_dist_params = {"a": np.array([1.0, 2.0])} expected_rv_op_params = {"a": np.array([1.0, 2.0])} - tests_to_run = ["check_pymc_params_match_rv_op"] + sizes_to_check = [None, (1), (4,), (3, 4)] + sizes_expected = [(2,), (1, 2), (4, 2), (3, 4, 2)] + reference_dist_params = {"alpha": np.array([1.0, 2.0])} + reference_dist = seeded_numpy_distribution_builder("dirichlet") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] class TestMultinomial(BaseTestDistribution): @@ -776,16 +778,135 @@ class TestMultinomial(BaseTestDistribution): expected_rv_op_params = {"n": 85, "p": np.array([0.28, 0.62, 0.10])} sizes_to_check = [None, (1), (4,), (3, 2)] sizes_expected = [(3,), (1, 3), (4, 3), (3, 2, 3)] - tests_to_run = ["check_pymc_params_match_rv_op", "check_rv_size"] + reference_dist_params = {"n": 85, "pvals": np.array([0.28, 0.62, 0.10])} + reference_dist = seeded_numpy_distribution_builder("multinomial") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] class TestCategorical(BaseTestDistribution): pymc_dist = pm.Categorical pymc_dist_params = {"p": np.array([0.28, 0.62, 0.10])} expected_rv_op_params = {"p": np.array([0.28, 0.62, 0.10])} + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_rv_size", + ] + + +class TestGeometric(BaseTestDistribution): + pymc_dist = pm.Geometric + pymc_dist_params = {"p": 0.9} + expected_rv_op_params = {"p": 0.9} tests_to_run = ["check_pymc_params_match_rv_op"] +class TestHyperGeometric(BaseTestDistribution): + pymc_dist = pm.HyperGeometric + pymc_dist_params = {"N": 20, "k": 12, "n": 5} + expected_rv_op_params = { + "ngood": pymc_dist_params["k"], + "nbad": pymc_dist_params["N"] - pymc_dist_params["k"], + "nsample": pymc_dist_params["n"], + } + reference_dist_params = expected_rv_op_params + reference_dist = seeded_numpy_distribution_builder("hypergeometric") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] + + +class TestLogistic(BaseTestDistribution): + pymc_dist = pm.Logistic + pymc_dist_params = {"mu": 1.0, "s": 2.0} + expected_rv_op_params = {"mu": 1.0, "s": 2.0} + tests_to_run = ["check_pymc_params_match_rv_op"] + + +class TestLognormal(BaseTestDistribution): + pymc_dist = pm.Lognormal + pymc_dist_params = {"mu": 1.0, "sigma": 5.0} + expected_rv_op_params = {"mu": 1.0, "sigma": 5.0} + tests_to_run = ["check_pymc_params_match_rv_op"] + + +class TestLognormalTau(BaseTestDistribution): + pymc_dist = pm.Lognormal + tau, sigma = get_tau_sigma(tau=25.0) + pymc_dist_params = {"mu": 1.0, "tau": 25.0} + expected_rv_op_params = {"mu": 1.0, "sigma": sigma} + tests_to_run = ["check_pymc_params_match_rv_op"] + + +class TestLognormalSd(BaseTestDistribution): + pymc_dist = pm.Lognormal + pymc_dist_params = {"mu": 1.0, "sd": 5.0} + expected_rv_op_params = {"mu": 1.0, "sigma": 5.0} + tests_to_run = ["check_pymc_params_match_rv_op"] + + +class TestTriangular(BaseTestDistribution): + pymc_dist = pm.Triangular + pymc_dist_params = {"lower": 0, "upper": 1, "c": 0.5} + expected_rv_op_params = {"lower": 0, "c": 0.5, "upper": 1} + reference_dist_params = {"left": 0, "mode": 0.5, "right": 1} + reference_dist = seeded_numpy_distribution_builder("triangular") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + ] + + +class TestVonMises(BaseTestDistribution): + pymc_dist = pm.VonMises + pymc_dist_params = {"mu": -2.1, "kappa": 5} + expected_rv_op_params = {"mu": -2.1, "kappa": 5} + tests_to_run = ["check_pymc_params_match_rv_op"] + + +class TestWeibull(BaseTestDistribution): + def weibull_rng_fn(self, size, alpha, beta, std_weibull_rng_fct): + return beta * std_weibull_rng_fct(alpha, size=size) + + def seeded_weibul_rng_fn(self): + std_weibull_rng_fct = functools.partial( + getattr(np.random.RandomState, "weibull"), self.get_random_state() + ) + return functools.partial(self.weibull_rng_fn, std_weibull_rng_fct=std_weibull_rng_fct) + + pymc_dist = pm.Weibull + pymc_dist_params = {"alpha": 1.0, "beta": 2.0} + expected_rv_op_params = {"alpha": 1.0, "beta": 2.0} + reference_dist_params = {"alpha": 1.0, "beta": 2.0} + reference_dist = seeded_weibul_rng_fn + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] + + +@pytest.mark.skipif( + condition=(SCIPY_VERSION < parse("1.4.0")), + reason="betabinom is new in Scipy 1.4.0", +) +class TestBetaBinomial(BaseTestDistribution): + pymc_dist = pm.BetaBinomial + pymc_dist_params = {"alpha": 2.0, "beta": 1.0, "n": 5} + expected_rv_op_params = {"n": 5, "alpha": 2.0, "beta": 1.0} + reference_dist_params = {"n": 5, "a": 2.0, "b": 1.0} + reference_dist = seeded_scipy_distribution_builder("betabinom") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] + + class TestScalarParameterSamples(SeededTest): @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_bounded(self): @@ -867,13 +988,6 @@ def ref_rand(size, kappa, b, mu): pymc3_random(pm.AsymmetricLaplace, {"b": Rplus, "kappa": Rplus, "mu": R}, ref_rand=ref_rand) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - def test_lognormal(self): - def ref_rand(size, mu, tau): - return np.exp(mu + (tau ** -0.5) * st.norm.rvs(loc=0.0, scale=1.0, size=size)) - - pymc3_random(pm.Lognormal, {"mu": R, "tau": Rplusbig}, ref_rand=ref_rand) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_student_t(self): def ref_rand(size, nu, mu, lam): @@ -888,24 +1002,6 @@ def ref_rand(size, mu, sigma, nu): pymc3_random(pm.ExGaussian, {"mu": R, "sigma": Rplus, "nu": Rplus}, ref_rand=ref_rand) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - def test_vonmises(self): - def ref_rand(size, mu, kappa): - return st.vonmises.rvs(size=size, loc=mu, kappa=kappa) - - pymc3_random(pm.VonMises, {"mu": R, "kappa": Rplus}, ref_rand=ref_rand) - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - def test_triangular(self): - def ref_rand(size, lower, upper, c): - scale = upper - lower - c_ = (c - lower) / scale - return st.triang.rvs(size=size, loc=lower, scale=scale, c=c_) - - pymc3_random( - pm.Triangular, {"lower": Runif, "upper": Runif + 3, "c": Runif + 1}, ref_rand=ref_rand - ) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_flat(self): with pm.Model(): @@ -920,40 +1016,6 @@ def test_half_flat(self): with pytest.raises(ValueError): f.random(1) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - @pytest.mark.xfail( - sys.platform.startswith("win"), - reason="Known issue: https://github.com/pymc-devs/pymc3/pull/4269", - ) - def test_beta_binomial(self): - pymc3_random_discrete( - pm.BetaBinomial, {"n": Nat, "alpha": Rplus, "beta": Rplus}, ref_rand=self._beta_bin - ) - - def _beta_bin(self, n, alpha, beta, size=None): - return st.binom.rvs(n, st.beta.rvs(a=alpha, b=beta, size=size)) - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - def test_geometric(self): - pymc3_random_discrete(pm.Geometric, {"p": Unit}, size=500, fails=50, ref_rand=nr.geometric) - - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - def test_hypergeometric(self): - def ref_rand(size, N, k, n): - return st.hypergeom.rvs(M=N, n=k, N=n, size=size) - - pymc3_random_discrete( - pm.HyperGeometric, - { - "N": Domain([10, 11, 12, 13], "int64"), - "k": Domain([4, 5, 6, 7], "int64"), - "n": Domain([6, 7, 8, 9], "int64"), - }, - size=500, - fails=50, - ref_rand=ref_rand, - ) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_discrete_uniform(self): def ref_rand(size, lower, upper): @@ -963,16 +1025,6 @@ def ref_rand(size, lower, upper): pm.DiscreteUniform, {"lower": -NatSmall, "upper": NatSmall}, ref_rand=ref_rand ) - def test_discrete_weibull(self): - def ref_rand(size, q, beta): - u = np.random.uniform(size=size) - - return np.ceil(np.power(np.log(1 - u) / np.log(q), 1.0 / beta)) - 1 - - pymc3_random_discrete( - pm.DiscreteWeibull, {"q": Unit, "beta": Rplusdunif}, ref_rand=ref_rand - ) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_constant_dist(self): def ref_rand(size, c): @@ -1189,13 +1241,6 @@ def test_dirichlet_multinomial_dist_ShapeError(self, n, a, shape, expectation): with expectation: m.random() - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - def test_logistic(self): - def ref_rand(size, mu, s): - return st.logistic.rvs(loc=mu, scale=s, size=size) - - pymc3_random(pm.Logistic, {"mu": R, "s": Rplus}, ref_rand=ref_rand) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_logitnormal(self): def ref_rand(size, mu, sigma): diff --git a/pymc3/tests/test_transforms.py b/pymc3/tests/test_transforms.py index e040e6e244..c2a8f0eb00 100644 --- a/pymc3/tests/test_transforms.py +++ b/pymc3/tests/test_transforms.py @@ -368,10 +368,30 @@ def transform_params(rv_var): ) self.check_transform_elementwise_logp(model) + @pytest.mark.parametrize( + "lower, c, upper, size", + [ + (0.0, 1.0, 2.0, 2), + (-10, 0, 200, (2, 3)), + (np.zeros(3), np.ones(3), np.ones(3), (4, 3)), + ], + ) + def test_triangular(self, lower, c, upper, size): + def transform_params(rv_var): + _, _, _, lower, _, upper = rv_var.owner.inputs + lower = at.as_tensor_variable(lower) if lower is not None else None + upper = at.as_tensor_variable(upper) if upper is not None else None + return lower, upper + + interval = tr.Interval(transform_params) + model = self.build_model( + pm.Triangular, {"lower": lower, "c": c, "upper": upper}, size=size, transform=interval + ) + self.check_transform_elementwise_logp(model) + @pytest.mark.parametrize( "mu,kappa,shape", [(0.0, 1.0, 2), (-0.5, 5.5, (2, 3)), (np.zeros(3), np.ones(3), (4, 3))] ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_vonmises(self, mu, kappa, shape): model = self.build_model( pm.VonMises, {"mu": mu, "kappa": kappa}, shape=shape, transform=tr.circular @@ -470,7 +490,6 @@ def transform_params(rv_var): @pytest.mark.parametrize( "mu,kappa,shape", [(0.0, 1.0, (2,)), (np.zeros(3), np.ones(3), (4, 3))] ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_vonmises_ordered(self, mu, kappa, shape): testval = np.sort(np.abs(np.random.rand(*shape))) model = self.build_model( @@ -514,3 +533,12 @@ def test_mvnormal_ordered(self, mu, cov, size, shape): pm.MvNormal, {"mu": mu, "cov": cov}, size=size, testval=testval, transform=tr.ordered ) self.check_vectortransform_elementwise_logp(model, vect_opt=1) + + +def test_triangular_transform(): + with pm.Model() as m: + x = pm.Triangular("x", lower=0, c=1, upper=2) + + transform = x.tag.value_var.tag.transform + assert np.isclose(transform.backward(x, -np.inf).eval(), 0) + assert np.isclose(transform.backward(x, np.inf).eval(), 2)