Bring back univariate continuous logp methods

Also add parameter check to uniform logcdf

Author: ricardoV94

pymc/distributions/continuous.py

@@ -12,6 +12,8 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
+# Contains code from Aeppl, Copyright (c) 2021-2022, Aesara Developers.
+
 # coding: utf-8
 """
 A collection of common probability distributions for stochastic
@@ -26,7 +28,7 @@
 import aesara.tensor as at
 import numpy as np
 
-from aeppl.logprob import _logprob, logcdf
+from aeppl.logprob import _logprob, logcdf, logprob
 from aesara.graph.basic import Apply, Variable
 from aesara.graph.op import Op
 from aesara.raise_op import Assert
@@ -311,9 +313,22 @@ def moment(rv, size, lower, upper):
             moment = at.full(size, moment)
         return moment
 
+    def logp(value, lower, upper):
+        res = at.switch(
+            at.bitwise_and(at.ge(value, lower), at.le(value, upper)),
+            at.fill(value, -at.log(upper - lower)),
+            -np.inf,
+        )
+
+        return check_parameters(
+            res,
+            lower <= upper,
+            msg="lower <= upper",
+        )
+
     def logcdf(value, lower, upper):
-        return at.switch(
-            at.lt(value, lower) | at.lt(upper, lower),
+        res = at.switch(
+            at.lt(value, lower),
             -np.inf,
             at.switch(
                 at.lt(value, upper),
@@ -322,6 +337,12 @@ def logcdf(value, lower, upper):
             ),
         )
 
+        return check_parameters(
+            res,
+            lower <= upper,
+            msg="lower <= upper",
+        )
+
 
 @_default_transform.register(Uniform)
 def uniform_default_transform(op, rv):
@@ -495,6 +516,14 @@ def moment(rv, size, mu, sigma):
             mu = at.full(size, mu)
         return mu
 
+    def logp(value, mu, sigma):
+        res = -0.5 * at.pow((value - mu) / sigma, 2) - at.log(at.sqrt(2.0 * np.pi)) - at.log(sigma)
+        return check_parameters(
+            res,
+            sigma > 0,
+            msg="sigma > 0",
+        )
+
     def logcdf(value, mu, sigma):
         return check_parameters(
             normal_lcdf(mu, sigma, value),
@@ -780,6 +809,15 @@ def moment(rv, size, loc, sigma):
             moment = at.full(size, moment)
         return moment
 
+    def logp(value, loc, sigma):
+        res = -0.5 * at.pow((value - loc) / sigma, 2) + at.log(at.sqrt(2.0 / np.pi)) - at.log(sigma)
+        res = at.switch(at.ge(value, loc), res, -np.inf)
+        return check_parameters(
+            res,
+            sigma > 0,
+            msg="sigma > 0",
+        )
+
     def logcdf(value, loc, sigma):
         z = zvalue(value, mu=loc, sigma=sigma)
         logcdf = at.switch(
@@ -1079,6 +1117,20 @@ def get_alpha_beta(self, alpha=None, beta=None, mu=None, sigma=None):
 
         return alpha, beta
 
+    def logp(value, alpha, beta):
+        res = (
+            at.switch(at.eq(alpha, 1.0), 0.0, (alpha - 1.0) * at.log(value))
+            + at.switch(at.eq(beta, 1.0), 0.0, (beta - 1.0) * at.log1p(-value))
+            - (at.gammaln(alpha) + at.gammaln(beta) - at.gammaln(alpha + beta))
+        )
+        res = at.switch(at.bitwise_and(at.ge(value, 0.0), at.le(value, 1.0)), res, -np.inf)
+        return check_parameters(
+            res,
+            alpha > 0,
+            beta > 0,
+            msg="alpha > 0, beta > 0",
+        )
+
     def logcdf(value, alpha, beta):
         logcdf = at.switch(
             at.lt(value, 0),
@@ -1261,6 +1313,15 @@ def moment(rv, size, mu):
             mu = at.full(size, mu)
         return mu
 
+    def logp(value, mu):
+        res = -at.log(mu) - value / mu
+        res = at.switch(at.ge(value, 0.0), res, -np.inf)
+        return check_parameters(
+            res,
+            mu >= 0,
+            msg="mu >= 0",
+        )
+
     def logcdf(value, mu):
         lam = at.reciprocal(mu)
         res = at.switch(
@@ -1334,6 +1395,14 @@ def moment(rv, size, mu, b):
             mu = at.full(size, mu)
         return mu
 
+    def logp(value, mu, b):
+        res = -at.log(2 * b) - at.abs(value - mu) / b
+        return check_parameters(
+            res,
+            b > 0,
+            msg="b > 0",
+        )
+
     def logcdf(value, mu, b):
         y = (value - mu) / b
 
@@ -1524,6 +1593,20 @@ def moment(rv, size, mu, sigma):
             mean = at.full(size, mean)
         return mean
 
+    def logp(value, mu, sigma):
+        res = (
+            -0.5 * at.pow((at.log(value) - mu) / sigma, 2)
+            - 0.5 * at.log(2.0 * np.pi)
+            - at.log(sigma)
+            - at.log(value)
+        )
+        res = at.switch(at.gt(value, 0.0), res, -np.inf)
+        return check_parameters(
+            res,
+            sigma > 0,
+            msg="sigma > 0",
+        )
+
     def logcdf(value, mu, sigma):
         res = at.switch(
             at.le(value, 0),
@@ -1732,6 +1815,16 @@ def moment(rv, size, alpha, m):
             median = at.full(size, median)
         return median
 
+    def logp(value, alpha, m):
+        res = at.log(alpha) + logpow(m, alpha) - logpow(value, alpha + 1.0)
+        res = at.switch(at.ge(value, m), res, -np.inf)
+        return check_parameters(
+            res,
+            alpha > 0,
+            m > 0,
+            msg="alpha > 0, m > 0",
+        )
+
     def logcdf(value, alpha, m):
         arg = (m / value) ** alpha
 
@@ -1819,6 +1912,14 @@ def moment(rv, size, alpha, beta):
             alpha = at.full(size, alpha)
         return alpha
 
+    def logp(value, alpha, beta):
+        res = -at.log(np.pi) - at.log(beta) - at.log1p(at.pow((value - alpha) / beta, 2))
+        return check_parameters(
+            res,
+            beta > 0,
+            msg="beta > 0",
+        )
+
     def logcdf(value, alpha, beta):
         res = at.log(0.5 + at.arctan((value - alpha) / beta) / np.pi)
         return check_parameters(
@@ -1879,6 +1980,15 @@ def moment(rv, size, loc, beta):
             beta = at.full(size, beta)
         return beta
 
+    def logp(value, loc, beta):
+        res = at.log(2) + logprob(Cauchy.dist(loc, beta), value)
+        res = at.switch(at.ge(value, loc), res, -np.inf)
+        return check_parameters(
+            res,
+            beta > 0,
+            msg="beta > 0",
+        )
+
     def logcdf(value, loc, beta):
         res = at.switch(
             at.lt(value, loc),
@@ -1990,6 +2100,17 @@ def moment(rv, size, alpha, inv_beta):
             mean = at.full(size, mean)
         return mean
 
+    def logp(value, alpha, inv_beta):
+        beta = at.reciprocal(inv_beta)
+        res = -at.gammaln(alpha) + logpow(beta, alpha) - beta * value + logpow(value, alpha - 1)
+        res = at.switch(at.ge(value, 0.0), res, -np.inf)
+        return check_parameters(
+            res,
+            alpha > 0,
+            beta > 0,
+            msg="alpha > 0, beta > 0",
+        )
+
     def logcdf(value, alpha, inv_beta):
         beta = at.reciprocal(inv_beta)
         res = at.switch(
@@ -2091,6 +2212,16 @@ def _get_alpha_beta(cls, alpha, beta, mu, sigma):
 
         return alpha, beta
 
+    def logp(value, alpha, beta):
+        res = -at.gammaln(alpha) + logpow(beta, alpha) - beta / value + logpow(value, -alpha - 1)
+        res = at.switch(at.ge(value, 0.0), res, -np.inf)
+        return check_parameters(
+            res,
+            alpha > 0,
+            beta > 0,
+            msg="alpha > 0, beta > 0",
+        )
+
     def logcdf(value, alpha, beta):
         res = at.switch(
             at.lt(value, 0),
@@ -2158,6 +2289,9 @@ def moment(rv, size, nu):
             moment = at.full(size, moment)
         return moment
 
+    def logp(value, nu):
+        return logprob(Gamma.dist(alpha=nu / 2, beta=0.5), value)
+
     def logcdf(value, nu):
         return logcdf(Gamma.dist(alpha=nu / 2, beta=0.5), value)
 
@@ -2586,6 +2720,15 @@ def moment(rv, size, mu, kappa):
             mu = at.full(size, mu)
         return mu
 
+    def logp(value, mu, kappa):
+        res = kappa * at.cos(mu - value) - at.log(2 * np.pi) - at.log(at.i0(kappa))
+        res = at.switch(at.bitwise_and(at.ge(value, -np.pi), at.le(value, np.pi)), res, -np.inf)
+        return check_parameters(
+            res,
+            kappa > 0,
+            msg="kappa > 0",
+        )
+
 
 class SkewNormalRV(RandomVariable):
     name = "skewnormal"
@@ -2771,6 +2914,20 @@ def moment(rv, size, lower, c, upper):
             mean = at.full(size, mean)
         return mean
 
+    def logp(value, lower, c, upper):
+        res = at.switch(
+            at.lt(value, c),
+            at.log(2 * (value - lower) / ((upper - lower) * (c - lower))),
+            at.log(2 * (upper - value) / ((upper - lower) * (upper - c))),
+        )
+        res = at.switch(at.bitwise_and(at.le(lower, value), at.le(value, upper)), res, -np.inf)
+        return check_parameters(
+            res,
+            lower <= c,
+            c <= upper,
+            msg="lower <= c <= upper",
+        )
+
     def logcdf(value, lower, c, upper):
         res = at.switch(
             at.le(value, lower),
@@ -2863,6 +3020,15 @@ def moment(rv, size, mu, beta):
             mean = at.full(size, mean)
         return mean
 
+    def logp(value, mu, beta):
+        z = (value - mu) / beta
+        res = -z - at.exp(-z) - at.log(beta)
+        return check_parameters(
+            res,
+            beta > 0,
+            msg="beta > 0",
+        )
+
     def logcdf(value, mu, beta):
         res = -at.exp(-(value - mu) / beta)
 
@@ -3062,6 +3228,15 @@ def moment(rv, size, mu, s):
             mu = at.full(size, mu)
         return mu
 
+    def logp(value, mu, s):
+        z = (value - mu) / s
+        res = -z - at.log(s) - 2.0 * at.log1p(at.exp(-z))
+        return check_parameters(
+            res,
+            s > 0,
+            msg="s > 0",
+        )
+
     def logcdf(value, mu, s):
         res = -at.log1pexp(-(value - mu) / s)
 

pymc/tests/distributions/test_continuous.py

@@ -185,8 +185,10 @@ def test_uniform(self):
         # Custom logp / logcdf check for invalid parameters
         invalid_dist = pm.Uniform.dist(lower=1, upper=0)
         with aesara.config.change_flags(mode=Mode("py")):
-            assert logp(invalid_dist, np.array(0.5)).eval() == -np.inf
-            assert logcdf(invalid_dist, np.array(2.0)).eval() == -np.inf
+            with pytest.raises(ParameterValueError):
+                logp(invalid_dist, np.array(0.5)).eval()
+            with pytest.raises(ParameterValueError):
+                logcdf(invalid_dist, np.array(0.5)).eval()
 
     def test_triangular(self):
         check_logp(