Skip to content

Commit

Permalink
Fix exponential and gamma logp / random link
Browse files Browse the repository at this point in the history
  • Loading branch information
ricardoV94 authored and brandonwillard committed Mar 29, 2021
1 parent c9fa127 commit 862693f
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 30 deletions.
58 changes: 31 additions & 27 deletions pymc3/distributions/continuous.py
Original file line number Diff line number Diff line change
Expand Up @@ -1377,49 +1377,53 @@ class Exponential(PositiveContinuous):
@classmethod
def dist(cls, lam, *args, **kwargs):
lam = at.as_tensor_variable(floatX(lam))
# mean = 1.0 / lam
# median = mean * at.log(2)
# mode = at.zeros_like(lam)

# variance = lam ** -2

assert_negative_support(lam, "lam", "Exponential")
return super().dist([lam], **kwargs)

def logp(value, lam):
# Aesara exponential op is parametrized in terms of mu (1/lam)
return super().dist([at.inv(lam)], **kwargs)

def logp(value, mu):
"""
Calculate log-probability of Exponential distribution at specified value.
Parameters
----------
value: numeric
Value(s) for which log-probability is calculated. If the log probabilities for multiple
values are desired the values must be provided in a numpy array or aesara tensor
Value(s) for which log-probability is calculated. If the log
probabilities for multiple values are desired the values must be
provided in a numpy array or aesara tensor
Returns
-------
TensorVariable
"""
return bound(at.log(lam) - lam * value, value >= 0, lam > 0)
lam = at.inv(mu)
return bound(
at.log(lam) - lam * value,
value >= 0,
lam > 0,
)

def logcdf(value, lam):
def logcdf(value, mu):
r"""
Compute the log of cumulative distribution function for the Exponential distribution
at the specified value.
Parameters
----------
value: numeric or np.ndarray or aesara.tensor
Value(s) for which log CDF is calculated. If the log CDF for multiple
values are desired the values must be provided in a numpy array or aesara tensor.
Value(s) for which log CDF is calculated. If the log CDF for
multiple values are desired the values must be provided in a numpy
array or aesara tensor.
Returns
-------
TensorVariable
"""
a = lam * value
lam = at.inv(mu)
return bound(
log1mexp(a),
log1mexp(lam * value),
0 <= value,
0 <= lam,
)
Expand Down Expand Up @@ -2371,14 +2375,12 @@ def dist(cls, alpha=None, beta=None, mu=None, sigma=None, sd=None, no_assert=Fal
alpha, beta = cls.get_alpha_beta(alpha, beta, mu, sigma)
alpha = at.as_tensor_variable(floatX(alpha))
beta = at.as_tensor_variable(floatX(beta))
# mean = alpha / beta
# mode = at.maximum((alpha - 1) / beta, 0)
# variance = alpha / beta ** 2

if not no_assert:
assert_negative_support(alpha, "alpha", "Gamma")
assert_negative_support(beta, "beta", "Gamma")

# Aesara exponential op is parametrized in terms of 1/lam
return super().dist([alpha, at.inv(beta)], **kwargs)

@classmethod
Expand All @@ -2397,45 +2399,47 @@ def get_alpha_beta(cls, alpha=None, beta=None, mu=None, sigma=None):

return alpha, beta

def _distr_parameters_for_repr(self):
return ["alpha", "beta"]

def logp(value, alpha, beta):
def logp(value, alpha, inv_beta):
"""
Calculate log-probability of Gamma distribution at specified value.
Parameters
----------
value: numeric
Value(s) for which log-probability is calculated. If the log probabilities for multiple
values are desired the values must be provided in a numpy array or `TensorVariable`.
Value(s) for which log-probability is calculated. If the log
probabilities for multiple values are desired the values must be
provided in a numpy array or `TensorVariable`.
Returns
-------
TensorVariable
"""
beta = at.inv(inv_beta)
return bound(
-gammaln(alpha) + logpow(beta, alpha) - beta * value + logpow(value, alpha - 1),
value >= 0,
alpha > 0,
beta > 0,
)

def logcdf(value, alpha, beta):
def logcdf(value, alpha, inv_beta):
"""
Compute the log of the cumulative distribution function for Gamma distribution
at the specified value.
Parameters
----------
value: numeric or np.ndarray or `TensorVariable`
Value(s) for which log CDF is calculated. If the log CDF for multiple
values are desired the values must be provided in a numpy array or `TensorVariable`.
Value(s) for which log CDF is calculated. If the log CDF for
multiple values are desired the values must be provided in a numpy
array or `TensorVariable`.
Returns
-------
TensorVariable
"""
beta = at.inv(inv_beta)

# Avoid C-assertion when the gammainc function is called with invalid values (#4340)
safe_alpha = at.switch(at.lt(alpha, 0), 0, alpha)
safe_beta = at.switch(at.lt(beta, 0), 0, beta)
Expand Down
6 changes: 3 additions & 3 deletions pymc3/tests/test_distributions_random.py
Original file line number Diff line number Diff line change
Expand Up @@ -626,7 +626,7 @@ def ref_rand(size, alpha, beta):

pymc3_random(pm.Beta, {"alpha": Rplus, "beta": Rplus}, ref_rand=ref_rand)

@pytest.mark.skip(reason="This test is covered by Aesara")
# @pytest.mark.skip(reason="This test is covered by Aesara")
def test_exponential(self):
def ref_rand(size, lam):
return nr.exponential(scale=1.0 / lam, size=size)
Expand Down Expand Up @@ -680,14 +680,14 @@ def ref_rand(size, beta):

pymc3_random(pm.HalfCauchy, {"beta": Rplusbig}, ref_rand=ref_rand)

@pytest.mark.skip(reason="This test is covered by Aesara")
# @pytest.mark.skip(reason="This test is covered by Aesara")
def test_gamma_alpha_beta(self):
def ref_rand(size, alpha, beta):
return st.gamma.rvs(alpha, scale=1.0 / beta, size=size)

pymc3_random(pm.Gamma, {"alpha": Rplusbig, "beta": Rplusbig}, ref_rand=ref_rand)

@pytest.mark.skip(reason="This test is covered by Aesara")
# @pytest.mark.skip(reason="This test is covered by Aesara")
def test_gamma_mu_sigma(self):
def ref_rand(size, mu, sigma):
return st.gamma.rvs(mu ** 2 / sigma ** 2, scale=sigma ** 2 / mu, size=size)
Expand Down

0 comments on commit 862693f

Please sign in to comment.