diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 88cff6779f..c60898ab16 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -1669,14 +1669,14 @@ class AsymmetricLaplace(Continuous): The pdf of this distribution is .. math:: - {f(x|\\b,\kappa) = - \left({\frac{\\b}{\kappa + 1/\kappa}}\right)\,e^{-(x)\\b\,s\kappa ^{s}}} + {f(x|\\b,\kappa,\mu) = + \left({\frac{\\b}{\kappa + 1/\kappa}}\right)\,e^{-(x-\mu)\\b\,s\kappa ^{s}}} where .. math:: - s = sgn(x) + s = sgn(x-\mu) ======== ======================== Support :math:`x \in \mathbb{R}` @@ -1706,13 +1706,16 @@ def __init__(self, b, kappa, mu=0, *args, **kwargs): self.mean = self.mu - (self.kappa - 1 / self.kappa) / b self.variance = (1 + self.kappa ** 4) / (self.kappa ** 2 * self.b ** 2) + assert_negative_support(kappa, "kappa", "AsymmetricLaplace") + assert_negative_support(b, "b", "AsymmetricLaplace") + super().__init__(*args, **kwargs) - def _random(self, b, kappa, size=None): + def _random(self, b, kappa, mu, size=None): u = np.random.uniform(size=size) switch = kappa ** 2 / (1 + kappa ** 2) - non_positive_x = kappa * np.log(u * (1 / switch)) / b - positive_x = -np.log((1 - u) * (1 + kappa ** 2)) / (kappa * b) + 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 @@ -1733,8 +1736,8 @@ def random(self, point=None, size=None): ------- array """ - b, kappa = draw_values([self.b, self.kappa], point=point, size=size) - return generate_samples(self._random, b, kappa, dist_shape=self.shape, size=size) + 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) def logp(self, value): """ diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 7f7e1195e9..c2548423f6 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -220,8 +220,9 @@ def build_model(distfam, valuedomain, vardomains, extra_args=None): return m -def laplace_asymmetric_logpdf(value, kappa, b): +def laplace_asymmetric_logpdf(value, kappa, b, mu): kapinv = 1 / kappa + value = value - mu lPx = value * b * np.where(value >= 0, -kappa, kapinv) lPx += np.log(b / (kappa + kapinv)) return lPx @@ -999,7 +1000,7 @@ def test_laplace_asymmetric(self): self.pymc3_matches_scipy( AsymmetricLaplace, R, - {"b": Rplus, "kappa": Rplus}, + {"b": Rplus, "kappa": Rplus, "mu": R}, laplace_asymmetric_logpdf, ) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 19b7901f8e..6cc4fe8042 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -377,7 +377,7 @@ class TestLaplace(BaseTestCases.BaseTestCase): class TestAsymmetricLaplace(BaseTestCases.BaseTestCase): distribution = pm.AsymmetricLaplace - params = {"kappa": 1.0, "b": 1.0} + params = {"kappa": 1.0, "b": 1.0, "mu": 0.0} class TestLognormal(BaseTestCases.BaseTestCase): @@ -632,17 +632,15 @@ def ref_rand(size, mu, b): pymc3_random(pm.Laplace, {"mu": R, "b": Rplus}, ref_rand=ref_rand) def test_laplace_asymmetric(self): - def ref_rand(size, kappa, b): + def ref_rand(size, kappa, b, mu): u = np.random.uniform(size=size) - x = -np.log((1 - u) * (1 + kappa ** 2)) / (kappa * b) * ( - u > ((kappa ** 2) / (1 + kappa ** 2)) - ) + kappa * np.log(u * (1 + kappa ** 2) / (kappa ** 2)) / b * ( - u < ((kappa ** 2) / (1 + kappa ** 2)) - ) - - return x + 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}, ref_rand=ref_rand) + pymc3_random(pm.AsymmetricLaplace, {"b": Rplus, "kappa": Rplus, "mu": R}, ref_rand=ref_rand) def test_lognormal(self): def ref_rand(size, mu, tau):