From bfc384b06428c6805624a7d4bfd7104d48f25374 Mon Sep 17 00:00:00 2001 From: kc611 Date: Thu, 10 Jun 2021 08:10:41 -0400 Subject: [PATCH 1/5] Refactored continuous distributions to v4 (AsymmetricLaplace, HalfStudentT, ExGaussian, Interpolated) --- pymc3/distributions/continuous.py | 303 +++++++++++------------ pymc3/tests/test_distributions.py | 25 +- pymc3/tests/test_distributions_random.py | 13 +- 3 files changed, 161 insertions(+), 180 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 908a7b138a..31e6e9df77 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -1587,6 +1587,26 @@ def logcdf(value, mu, b): ) +class AsymmetricLaplaceRV(RandomVariable): + name = "asymmetriclaplace" + ndim_supp = 0 + ndims_params = [0, 0, 0] + dtype = "floatX" + _print_name = ("AsymmetricLaplace", "\\operatorname{AsymmetricLaplace}") + + @classmethod + def rng_fn(cls, rng, b, kappa, mu, size=None): + u = rng.uniform(size=size) + switch = kappa ** 2 / (1 + kappa ** 2) + non_positive_x = mu + kappa * np.log(u * (1 / switch)) / b + positive_x = mu - np.log((1 - u) * (1 + kappa ** 2)) / (kappa * b) + draws = non_positive_x * (u <= switch) + positive_x * (u > switch) + return draws + + +asymmetriclaplace = AsymmetricLaplaceRV() + + class AsymmetricLaplace(Continuous): r""" Asymmetric-Laplace log-likelihood. @@ -1622,49 +1642,23 @@ class AsymmetricLaplace(Continuous): -------- `Reference `_ """ + rv_op = asymmetriclaplace - def __init__(self, b, kappa, mu=0, *args, **kwargs): - self.b = at.as_tensor_variable(floatX(b)) - self.kappa = at.as_tensor_variable(floatX(kappa)) - self.mu = mu = at.as_tensor_variable(floatX(mu)) + @classmethod + def dist(cls, b, kappa, mu=0, *args, **kwargs): + b = at.as_tensor_variable(floatX(b)) + kappa = at.as_tensor_variable(floatX(kappa)) + mu = mu = at.as_tensor_variable(floatX(mu)) - self.mean = self.mu - (self.kappa - 1 / self.kappa) / b - self.variance = (1 + self.kappa ** 4) / (self.kappa ** 2 * self.b ** 2) + # mean = mu - (kappa - 1 / kappa) / b + # variance = (1 + kappa ** 4) / (kappa ** 2 * b ** 2) assert_negative_support(kappa, "kappa", "AsymmetricLaplace") assert_negative_support(b, "b", "AsymmetricLaplace") - super().__init__(*args, **kwargs) - - def _random(self, b, kappa, mu, size=None): - u = np.random.uniform(size=size) - switch = kappa ** 2 / (1 + kappa ** 2) - non_positive_x = mu + kappa * np.log(u * (1 / switch)) / b - positive_x = mu - np.log((1 - u) * (1 + kappa ** 2)) / (kappa * b) - draws = non_positive_x * (u <= switch) + positive_x * (u > switch) - return draws - - def random(self, point=None, size=None): - """ - Draw random samples from this distribution, using the inverse CDF method. - - 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 - """ - # b, kappa, mu = draw_values([self.b, self.kappa, self.mu], point=point, size=size) - # return generate_samples(self._random, b, kappa, mu, dist_shape=self.shape, size=size) + return super().dist([b, kappa, mu], *args, **kwargs) - def logp(self, value): + def logp(value, b, kappa, mu): """ Calculate log-probability of Asymmetric-Laplace distribution at specified value. @@ -1678,12 +1672,12 @@ def logp(self, value): ------- TensorVariable """ - value = value - self.mu + value = value - mu return bound( - at.log(self.b / (self.kappa + (self.kappa ** -1))) - + (-value * self.b * at.sgn(value) * (self.kappa ** at.sgn(value))), - 0 < self.b, - 0 < self.kappa, + at.log(b / (kappa + (kappa ** -1))) + + (-value * b * at.sgn(value) * (kappa ** at.sgn(value))), + 0 < b, + 0 < kappa, ) @@ -2759,6 +2753,21 @@ def logcdf(value, alpha, beta): ) +class HalfStudentTRV(RandomVariable): + name = "halfstudentt" + ndim_supp = 0 + ndims_params = [0, 0] + dtype = "floatX" + _print_name = ("HalfStudentT", "\\operatorname{HalfStudentT}") + + @classmethod + def rng_fn(cls, rng, nu, sigma, size=None): + return np.abs(stats.t.rvs(nu, sigma, size=size, random_state=rng)) + + +halfstudentt = HalfStudentTRV() + + class HalfStudentT(PositiveContinuous): r""" Half Student's T log-likelihood @@ -2816,46 +2825,29 @@ class HalfStudentT(PositiveContinuous): with pm.Model(): x = pm.HalfStudentT('x', lam=4, nu=10) """ + rv_op = halfstudentt + + @classmethod + def dist(cls, nu=1, sigma=None, lam=None, sd=None, *args, **kwargs): - def __init__(self, nu=1, sigma=None, lam=None, sd=None, *args, **kwargs): - super().__init__(*args, **kwargs) if sd is not None: sigma = sd - self.mode = at.as_tensor_variable(0) + nu = at.as_tensor_variable(floatX(nu)) lam, sigma = get_tau_sigma(lam, sigma) - self.median = at.as_tensor_variable(sigma) - self.sigma = self.sd = at.as_tensor_variable(sigma) - self.lam = at.as_tensor_variable(lam) - self.nu = nu = at.as_tensor_variable(floatX(nu)) + sigma = at.as_tensor_variable(sigma) - assert_negative_support(sigma, "sigma", "HalfStudentT") - assert_negative_support(lam, "lam", "HalfStudentT") - assert_negative_support(nu, "nu", "HalfStudentT") + # mode = at.as_tensor_variable(0) + # median = at.as_tensor_variable(sigma) + # sd = at.as_tensor_variable(sigma) - def random(self, point=None, size=None): - """ - Draw random values from HalfStudentT distribution. + assert_negative_support(nu, "nu", "HalfStudentT") + assert_negative_support(lam, "lam", "HalfStudentT") + assert_negative_support(sigma, "sigma", "HalfStudentT") - 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). + return super().dist([nu, sigma], *args, **kwargs) - Returns - ------- - array - """ - # nu, sigma = draw_values([self.nu, self.sigma], point=point, size=size) - # return np.abs( - # generate_samples(stats.t.rvs, nu, loc=0, scale=sigma, dist_shape=self.shape, size=size) - # ) - - def logp(self, value): + def logp(value, nu, sigma): """ Calculate log-probability of HalfStudentT distribution at specified value. @@ -2869,9 +2861,8 @@ def logp(self, value): ------- TensorVariable """ - nu = self.nu - sigma = self.sigma - lam = self.lam + + lam, sigma = get_tau_sigma(None, sigma) return bound( at.log(2) @@ -2889,6 +2880,21 @@ def _distr_parameters_for_repr(self): return ["nu", "lam"] +class ExGaussianRV(RandomVariable): + name = "exgaussian" + ndim_supp = 0 + ndims_params = [0, 0, 0] + dtype = "floatX" + _print_name = ("ExGaussian", "\\operatorname{ExGaussian}") + + @classmethod + def rng_fn(cls, rng, mu, sigma, nu, size=None): + return rng.normal(mu, sigma, size=size) + rng.exponential(scale=nu, size=size) + + +exgaussian = ExGaussianRV() + + class ExGaussian(Continuous): r""" Exponentially modified Gaussian log-likelihood. @@ -2954,49 +2960,28 @@ class ExGaussian(Continuous): Tutorials in Quantitative Methods for Psychology, Vol. 4, No. 1, pp 35-45. """ + rv_op = exgaussian - def __init__(self, mu=0.0, sigma=None, nu=None, sd=None, *args, **kwargs): - super().__init__(*args, **kwargs) + @classmethod + def dist(cls, mu=0.0, sigma=None, nu=None, sd=None, *args, **kwargs): if sd is not None: sigma = sd - self.mu = mu = at.as_tensor_variable(floatX(mu)) - self.sigma = self.sd = sigma = at.as_tensor_variable(floatX(sigma)) - self.nu = nu = at.as_tensor_variable(floatX(nu)) - self.mean = mu + nu - self.variance = (sigma ** 2) + (nu ** 2) + mu = at.as_tensor_variable(floatX(mu)) + sigma = at.as_tensor_variable(floatX(sigma)) + nu = at.as_tensor_variable(floatX(nu)) + + # sd = sigma + # mean = mu + nu + # variance = (sigma ** 2) + (nu ** 2) assert_negative_support(sigma, "sigma", "ExGaussian") assert_negative_support(nu, "nu", "ExGaussian") - def random(self, point=None, size=None): - """ - Draw random values from ExGaussian 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, sigma, nu = draw_values([self.mu, self.sigma, self.nu], point=point, size=size) - # - # def _random(mu, sigma, nu, size=None): - # return np.random.normal(mu, sigma, size=size) + np.random.exponential( - # scale=nu, size=size - # ) - # - # return generate_samples(_random, mu, sigma, nu, dist_shape=self.shape, size=size) + return super().dist([mu, sigma, nu], *args, **kwargs) - def logp(self, value): + def logp(value, mu, sigma, nu): """ Calculate log-probability of ExGaussian distribution at specified value. @@ -3010,9 +2995,6 @@ def logp(self, value): ------- TensorVariable """ - mu = self.mu - sigma = self.sigma - nu = self.nu # Alogithm is adapted from dexGAUS.R from gamlss return bound( @@ -3030,7 +3012,7 @@ def logp(self, value): 0 < nu, ) - def logcdf(self, value): + def logcdf(value, mu, sigma, nu): """ Compute the log of the cumulative distribution function for ExGaussian distribution at the specified value. @@ -3051,9 +3033,6 @@ def logcdf(self, value): ------- TensorVariable """ - mu = self.mu - sigma = self.sigma - nu = self.nu # Alogithm is adapted from pexGAUS.R from gamlss return bound( @@ -3822,6 +3801,33 @@ def logp(value, mu, sigma): ) +def _interpolated_argcdf(p, pdf, cdf, x): + index = np.searchsorted(cdf, p) - 1 + slope = (pdf[index + 1] - pdf[index]) / (x[index + 1] - x[index]) + + return x[index] + np.where( + np.abs(slope) <= 1e-8, + np.where(np.abs(pdf[index]) <= 1e-8, np.zeros(index.shape), (p - cdf[index]) / pdf[index]), + (-pdf[index] + np.sqrt(pdf[index] ** 2 + 2 * slope * (p - cdf[index]))) / slope, + ) + + +class InterpolatedRV(RandomVariable): + name = "interpolated" + ndim_supp = 0 + ndims_params = [1, 1, 1] + dtype = "floatX" + _print_name = ("Interpolated", "\\operatorname{Interpolated}") + + @classmethod + def rng_fn(cls, rng, x, pdf, cdf, size=None): + p = rng.uniform(size=size) + return _interpolated_argcdf(p, pdf, cdf, x) + + +interpolated = InterpolatedRV() + + class Interpolated(BoundedContinuous): r""" Univariate probability distribution defined as a linear interpolation @@ -3872,59 +3878,28 @@ class Interpolated(BoundedContinuous): Probability density function evaluated on lattice ``x_points`` """ - def __init__(self, x_points, pdf_points, *args, **kwargs): - self.lower = lower = at.as_tensor_variable(x_points[0]) - self.upper = upper = at.as_tensor_variable(x_points[-1]) + rv_op = interpolated - super().__init__(lower=lower, upper=upper, *args, **kwargs) + @classmethod + def dist(cls, x_points, pdf_points, *args, **kwargs): interp = InterpolatedUnivariateSpline(x_points, pdf_points, k=1, ext="zeros") - Z = interp.integral(x_points[0], x_points[-1]) - - self.Z = at.as_tensor_variable(Z) - self.interp_op = SplineWrapper(interp) - self.x_points = x_points - self.pdf_points = pdf_points / Z - self.cdf_points = interp.antiderivative()(x_points) / Z - - self.median = self._argcdf(0.5) - - def _argcdf(self, p): - pdf = self.pdf_points - cdf = self.cdf_points - x = self.x_points - index = np.searchsorted(cdf, p) - 1 - slope = (pdf[index + 1] - pdf[index]) / (x[index + 1] - x[index]) + Z = interp.integral(x_points[0], x_points[-1]) + cdf_points = interp.antiderivative()(x_points) / Z + pdf_points = pdf_points / Z - return x[index] + np.where( - np.abs(slope) <= 1e-8, - np.where( - np.abs(pdf[index]) <= 1e-8, np.zeros(index.shape), (p - cdf[index]) / pdf[index] - ), - (-pdf[index] + np.sqrt(pdf[index] ** 2 + 2 * slope * (p - cdf[index]))) / slope, - ) + x_points = at.as_tensor_variable(floatX(x_points)) + pdf_points = at.as_tensor_variable(floatX(pdf_points)) + cdf_points = at.as_tensor_variable(floatX(cdf_points)) - def _random(self, size=None): - return self._argcdf(np.random.uniform(size=size)) + # lower = at.as_tensor_variable(x_points[0]) + # upper = at.as_tensor_variable(x_points[-1]) + # median = _interpolated_argcdf(0.5, pdf_points, cdf_points, x_points) - def random(self, point=None, size=None): - """ - Draw random values from Interpolated distribution. + return super().dist([x_points, pdf_points, cdf_points], **kwargs) - Parameters - ---------- - size: int, optional - Desired size of random sample (returns one sample if not - specified). - - Returns - ------- - array - """ - # return generate_samples(self._random, dist_shape=self.shape, size=size) - - def logp(self, value): + def logp(value, x_points, pdf_points, cdf_points): """ Calculate log-probability of Interpolated distribution at specified value. @@ -3938,7 +3913,13 @@ def logp(self, value): ------- TensorVariable """ - return at.log(self.interp_op(value) / self.Z) + interp = InterpolatedUnivariateSpline(x_points.data, pdf_points.data, k=1, ext="zeros") + interp_op = SplineWrapper(interp) + + Z = interp.integral(x_points.data[0], x_points.data[-1]) + Z = at.as_tensor_variable(Z) + + return at.log(interp_op(value) / Z) def _distr_parameters_for_repr(self): return [] diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 30ee062291..fc4d9ab415 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -1309,7 +1309,6 @@ def test_laplace(self): lambda value, mu, b: sp.laplace.logcdf(value, mu, b), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_laplace_asymmetric(self): self.check_logp( AsymmetricLaplace, @@ -1518,7 +1517,6 @@ def test_weibull_logcdf(self): lambda value, alpha, beta: sp.exponweib.logcdf(value, 1, alpha, scale=beta), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_half_studentt(self): # this is only testing for nu=1 (halfcauchy) self.check_logp( @@ -2471,7 +2469,6 @@ def test_get_tau_sigma(self): (-1.0, 0.0, 0.1, 0.1, -51.022349), # Fails in previous pymc3 version ], ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_ex_gaussian(self, value, mu, sigma, nu, logp): """Log probabilities calculated using the dexGAUS function from the R package gamlss. See e.g., doi: 10.1111/j.1467-9876.2005.00510.x, or http://www.gamlss.org/.""" @@ -2486,7 +2483,7 @@ def test_ex_gaussian(self, value, mu, sigma, nu, logp): ) @pytest.mark.parametrize( - "value,mu,sigma,nu,logcdf", + "value,mu,sigma,nu,logcdf_val", [ (0.5, -50.000, 0.500, 0.500, 0.0000000), (1.0, -1.000, 0.001, 0.001, 0.0000000), @@ -2501,18 +2498,16 @@ def test_ex_gaussian(self, value, mu, sigma, nu, logp): (-0.72402009, 0.0, 0.1, 0.1, -31.26571842), # Previous 64-bit version failed here ], ) - @pytest.mark.xfail(reason="Distribution not refactored yet") - def test_ex_gaussian_cdf(self, value, mu, sigma, nu, logcdf): + def test_ex_gaussian_cdf(self, value, mu, sigma, nu, logcdf_val): """Log probabilities calculated using the pexGAUS function from the R package gamlss. See e.g., doi: 10.1111/j.1467-9876.2005.00510.x, or http://www.gamlss.org/.""" assert_almost_equal( - logcdf(ExGaussian.dist(mu=mu, sigma=sigma, nu=nu), value).tag.test_value, - logcdf, + logcdf(ExGaussian.dist(mu=mu, sigma=sigma, nu=nu), value).eval(), + logcdf_val, decimal=select_by_precision(float64=6, float32=2), - err_msg=str((value, mu, sigma, nu, logcdf)), + err_msg=str((value, mu, sigma, nu, logcdf_val)), ) - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_ex_gaussian_cdf_outside_edges(self): self.check_logcdf( ExGaussian, @@ -2615,7 +2610,6 @@ def test_moyal_logcdf(self): ) @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") - @pytest.mark.xfail(reason="Distribution not refactored yet") def test_interpolated(self): for mu in R.vals: for sigma in Rplus.vals: @@ -2623,11 +2617,16 @@ def test_interpolated(self): xmin = mu - 5 * sigma xmax = mu + 5 * sigma + from pymc3.distributions.continuous import interpolated + class TestedInterpolated(Interpolated): - def __init__(self, **kwargs): + rv_op = interpolated + + @classmethod + def dist(cls, **kwargs): x_points = np.linspace(xmin, xmax, 100000) pdf_points = sp.norm.pdf(x_points, loc=mu, scale=sigma) - super().__init__(x_points=x_points, pdf_points=pdf_points, **kwargs) + return super().dist(x_points=x_points, pdf_points=pdf_points, **kwargs) def ref_pdf(value): return np.where( diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index a4e2905d5b..4fa351866c 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -270,7 +270,6 @@ class TestWald(BaseTestCases.BaseTestCase): params = {"mu": 1.0, "lam": 1.0, "alpha": 0.0} -@pytest.mark.xfail(reason="This distribution has not been refactored for v4") class TestAsymmetricLaplace(BaseTestCases.BaseTestCase): distribution = pm.AsymmetricLaplace params = {"kappa": 1.0, "b": 1.0, "mu": 0.0} @@ -1257,7 +1256,6 @@ def ref_rand(size, mu, lam, alpha): ref_rand=ref_rand, ) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_laplace_asymmetric(self): def ref_rand(size, kappa, b, mu): u = np.random.uniform(size=size) @@ -1269,7 +1267,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_ex_gaussian(self): def ref_rand(size, mu, sigma, nu): return nr.normal(mu, sigma, size=size) + nr.exponential(scale=nu, size=size) @@ -1497,7 +1494,6 @@ def ref_rand(size, mu, sigma): pymc3_random(pm.Moyal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self): for mu in R.vals: @@ -1506,11 +1502,16 @@ def test_interpolated(self): def ref_rand(size): return st.norm.rvs(loc=mu, scale=sigma, size=size) + from pymc3.distributions.continuous import interpolated + class TestedInterpolated(pm.Interpolated): - def __init__(self, **kwargs): + rv_op = interpolated + + @classmethod + def dist(cls, **kwargs): x_points = np.linspace(mu - 5 * sigma, mu + 5 * sigma, 100) pdf_points = st.norm.pdf(x_points, loc=mu, scale=sigma) - super().__init__(x_points=x_points, pdf_points=pdf_points, **kwargs) + return super().dist(x_points=x_points, pdf_points=pdf_points, **kwargs) pymc3_random(TestedInterpolated, {}, ref_rand=ref_rand) From 7e34e8622b3cb41877d6b8fa97ba8e66d7a8e2a3 Mon Sep 17 00:00:00 2001 From: kc611 Date: Fri, 11 Jun 2021 01:47:53 -0400 Subject: [PATCH 2/5] Refactored tests for new distributions --- pymc3/tests/test_distributions_random.py | 130 ++++++++++++++--------- 1 file changed, 80 insertions(+), 50 deletions(-) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 4fa351866c..b2691e945e 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -32,7 +32,7 @@ import pymc3 as pm from pymc3.aesaraf import change_rv_size, floatX, intX -from pymc3.distributions.continuous import get_tau_sigma +from pymc3.distributions.continuous import get_tau_sigma, interpolated from pymc3.distributions.dist_math import clipped_beta_rvs from pymc3.distributions.multivariate import quaddist_matrix from pymc3.distributions.shape_utils import to_tuple @@ -270,17 +270,6 @@ class TestWald(BaseTestCases.BaseTestCase): params = {"mu": 1.0, "lam": 1.0, "alpha": 0.0} -class TestAsymmetricLaplace(BaseTestCases.BaseTestCase): - distribution = pm.AsymmetricLaplace - params = {"kappa": 1.0, "b": 1.0, "mu": 0.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 TestZeroInflatedNegativeBinomial(BaseTestCases.BaseTestCase): distribution = pm.ZeroInflatedNegativeBinomial @@ -464,6 +453,64 @@ class TestLaplace(BaseTestDistribution): ] +class TestAsymmetricLaplace(BaseTestDistribution): + def asymmetriclaplace_rng_fn(self, b, kappa, mu, size, uniform_rng_fct): + u = uniform_rng_fct(size=size) + switch = kappa ** 2 / (1 + kappa ** 2) + non_positive_x = mu + kappa * np.log(u * (1 / switch)) / b + positive_x = mu - np.log((1 - u) * (1 + kappa ** 2)) / (kappa * b) + draws = non_positive_x * (u <= switch) + positive_x * (u > switch) + return draws + + def seeded_asymmetriclaplace_rng_fn(self): + uniform_rng_fct = functools.partial( + getattr(np.random.RandomState, "uniform"), self.get_random_state() + ) + return functools.partial(self.asymmetriclaplace_rng_fn, uniform_rng_fct=uniform_rng_fct) + + pymc_dist = pm.AsymmetricLaplace + + pymc_dist_params = {"b": 1.0, "kappa": 1.0, "mu": 0.0} + expected_rv_op_params = {"b": 1.0, "kappa": 1.0, "mu": 0.0} + reference_dist_params = {"b": 1.0, "kappa": 1.0, "mu": 0.0} + reference_dist = seeded_asymmetriclaplace_rng_fn + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] + + +class TestExGaussian(BaseTestDistribution): + def exgaussian_rng_fn(self, mu, sigma, nu, size, normal_rng_fct, exponential_rng_fct): + return normal_rng_fct(mu, sigma, size=size) + exponential_rng_fct(scale=nu, size=size) + + def seeded_exgaussian_rng_fn(self): + normal_rng_fct = functools.partial( + getattr(np.random.RandomState, "normal"), self.get_random_state() + ) + exponential_rng_fct = functools.partial( + getattr(np.random.RandomState, "exponential"), self.get_random_state() + ) + return functools.partial( + self.exgaussian_rng_fn, + normal_rng_fct=normal_rng_fct, + exponential_rng_fct=exponential_rng_fct, + ) + + pymc_dist = pm.ExGaussian + + pymc_dist_params = {"mu": 1.0, "sigma": 1.0, "nu": 1.0} + expected_rv_op_params = {"mu": 1.0, "sigma": 1.0, "nu": 1.0} + reference_dist_params = {"mu": 1.0, "sigma": 1.0, "nu": 1.0} + reference_dist = seeded_exgaussian_rng_fn + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] + + class TestGumbel(BaseTestDistribution): pymc_dist = pm.Gumbel pymc_dist_params = {"mu": 1.5, "beta": 3.0} @@ -1195,6 +1242,27 @@ class TestOrderedProbit(BaseTestDistribution): ] +class TestInterpolated(SeededTest): + @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") + def test_interpolated(self): + for mu in R.vals: + for sigma in Rplus.vals: + # pylint: disable=cell-var-from-loop + def ref_rand(size): + return st.norm.rvs(loc=mu, scale=sigma, size=size) + + class TestedInterpolated(pm.Interpolated): + rv_op = interpolated + + @classmethod + def dist(cls, **kwargs): + x_points = np.linspace(mu - 5 * sigma, mu + 5 * sigma, 100) + pdf_points = st.norm.pdf(x_points, loc=mu, scale=sigma) + return super().dist(x_points=x_points, pdf_points=pdf_points, **kwargs) + + pymc3_random(TestedInterpolated, {}, ref_rand=ref_rand) + + class TestScalarParameterSamples(SeededTest): @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_bounded(self): @@ -1256,23 +1324,6 @@ def ref_rand(size, mu, lam, alpha): ref_rand=ref_rand, ) - def test_laplace_asymmetric(self): - def ref_rand(size, kappa, b, mu): - u = np.random.uniform(size=size) - switch = kappa ** 2 / (1 + kappa ** 2) - non_positive_x = mu + kappa * np.log(u * (1 / switch)) / b - positive_x = mu - np.log((1 - u) * (1 + kappa ** 2)) / (kappa * b) - draws = non_positive_x * (u <= switch) + positive_x * (u > switch) - return draws - - pymc3_random(pm.AsymmetricLaplace, {"b": Rplus, "kappa": Rplus, "mu": R}, ref_rand=ref_rand) - - def test_ex_gaussian(self): - def ref_rand(size, mu, sigma, nu): - return nr.normal(mu, sigma, size=size) + nr.exponential(scale=nu, size=size) - - 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_matrix_normal(self): def ref_rand(size, mu, rowcov, colcov): @@ -1494,27 +1545,6 @@ def ref_rand(size, mu, sigma): pymc3_random(pm.Moyal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand) - @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") - def test_interpolated(self): - for mu in R.vals: - for sigma in Rplus.vals: - # pylint: disable=cell-var-from-loop - def ref_rand(size): - return st.norm.rvs(loc=mu, scale=sigma, size=size) - - from pymc3.distributions.continuous import interpolated - - class TestedInterpolated(pm.Interpolated): - rv_op = interpolated - - @classmethod - def dist(cls, **kwargs): - x_points = np.linspace(mu - 5 * sigma, mu + 5 * sigma, 100) - pdf_points = st.norm.pdf(x_points, loc=mu, scale=sigma) - return super().dist(x_points=x_points, pdf_points=pdf_points, **kwargs) - - pymc3_random(TestedInterpolated, {}, ref_rand=ref_rand) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") @pytest.mark.skip( "Wishart random sampling not implemented.\n" From 53e4e0ef1bbcaf4ab50af58ca3e01e7aec0b9835 Mon Sep 17 00:00:00 2001 From: kc611 Date: Fri, 11 Jun 2021 04:37:38 -0400 Subject: [PATCH 3/5] Added RV size testing for Interpolated --- pymc3/tests/test_distributions_random.py | 29 ++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index b2691e945e..6af81782cc 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -32,7 +32,11 @@ import pymc3 as pm from pymc3.aesaraf import change_rv_size, floatX, intX -from pymc3.distributions.continuous import get_tau_sigma, interpolated +from pymc3.distributions.continuous import ( + _interpolated_argcdf, + get_tau_sigma, + interpolated, +) from pymc3.distributions.dist_math import clipped_beta_rvs from pymc3.distributions.multivariate import quaddist_matrix from pymc3.distributions.shape_utils import to_tuple @@ -1242,7 +1246,28 @@ class TestOrderedProbit(BaseTestDistribution): ] -class TestInterpolated(SeededTest): +class TestInterpolated(BaseTestDistribution): + def interpolated_rng_fn(self, size, mu, sigma, rng): + return st.norm.rvs(loc=mu, scale=sigma, size=size) + + pymc_dist = pm.Interpolated + + # Dummy values for RV size testing + mu = sigma = 1 + x_points = pdf_points = np.linspace(1, 100, 100) + + pymc_dist_params = {"x_points": x_points, "pdf_points": pdf_points} + reference_dist_params = {"mu": mu, "sigma": sigma} + + reference_dist = lambda self: functools.partial( + self.interpolated_rng_fn, rng=self.get_random_state() + ) + tests_to_run = [ + "check_rv_size", + ] + + +class TestInterpolatedSeeded(SeededTest): @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self): for mu in R.vals: From c421634627cb2508a4a128f6b165121d0b861cbb Mon Sep 17 00:00:00 2001 From: kc611 Date: Fri, 11 Jun 2021 04:41:43 -0400 Subject: [PATCH 4/5] Removed unused import --- pymc3/tests/test_distributions_random.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 6af81782cc..b907a7f517 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -32,11 +32,7 @@ import pymc3 as pm from pymc3.aesaraf import change_rv_size, floatX, intX -from pymc3.distributions.continuous import ( - _interpolated_argcdf, - get_tau_sigma, - interpolated, -) +from pymc3.distributions.continuous import get_tau_sigma, interpolated from pymc3.distributions.dist_math import clipped_beta_rvs from pymc3.distributions.multivariate import quaddist_matrix from pymc3.distributions.shape_utils import to_tuple From e684ede94d9b2a26eb2a4130111a1b9539f580ae Mon Sep 17 00:00:00 2001 From: Ricardo Date: Sat, 12 Jun 2021 19:21:35 +0200 Subject: [PATCH 5/5] Convert to constant and add comments --- pymc3/distributions/continuous.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 31e6e9df77..47920f3518 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -3873,9 +3873,10 @@ class Interpolated(BoundedContinuous): Parameters ---------- x_points: array-like - A monotonically growing list of values + A monotonically growing list of values. Must be non-symbolic pdf_points: array-like - Probability density function evaluated on lattice ``x_points`` + Probability density function evaluated on lattice ``x_points``. Must + be non-symbolic """ rv_op = interpolated @@ -3889,9 +3890,9 @@ def dist(cls, x_points, pdf_points, *args, **kwargs): cdf_points = interp.antiderivative()(x_points) / Z pdf_points = pdf_points / Z - x_points = at.as_tensor_variable(floatX(x_points)) - pdf_points = at.as_tensor_variable(floatX(pdf_points)) - cdf_points = at.as_tensor_variable(floatX(cdf_points)) + x_points = at.constant(floatX(x_points)) + pdf_points = at.constant(floatX(pdf_points)) + cdf_points = at.constant(floatX(cdf_points)) # lower = at.as_tensor_variable(x_points[0]) # upper = at.as_tensor_variable(x_points[-1]) @@ -3913,11 +3914,14 @@ def logp(value, x_points, pdf_points, cdf_points): ------- TensorVariable """ + # x_points and pdf_points are expected to be non-symbolic arrays wrapped + # within a tensor.constant. We use the .data method to retrieve them interp = InterpolatedUnivariateSpline(x_points.data, pdf_points.data, k=1, ext="zeros") - interp_op = SplineWrapper(interp) - Z = interp.integral(x_points.data[0], x_points.data[-1]) - Z = at.as_tensor_variable(Z) + + # interp and Z are converted to symbolic variables here + interp_op = SplineWrapper(interp) + Z = at.constant(Z) return at.log(interp_op(value) / Z)