add distributions from student t family

distributions/cauchy.py

@@ -0,0 +1,78 @@
+import pytensor.tensor as pt
+
+from .helper import ppf_bounds_cont
+
+
+def mean(alpha, beta):
+    alpha_b, _ = pt.broadcast_arrays(alpha, beta)
+    return pt.full_like(alpha_b, pt.nan)
+
+
+def mode(alpha, beta):
+    return pt.broadcast_arrays(alpha, beta)[0]
+
+
+def median(alpha, beta):
+    return pt.broadcast_arrays(alpha, beta)[0]
+
+
+def var(alpha, beta):
+    alpha_b, _ = pt.broadcast_arrays(alpha, beta)
+    return pt.full_like(alpha_b, pt.nan)
+
+
+def std(alpha, beta):
+    alpha_b, _ = pt.broadcast_arrays(alpha, beta)
+    return pt.full_like(alpha_b, pt.nan)
+
+
+def skewness(alpha, beta):
+    alpha_b, _ = pt.broadcast_arrays(alpha, beta)
+    return pt.full_like(alpha_b, pt.nan)
+
+
+def kurtosis(alpha, beta):
+    alpha_b, _ = pt.broadcast_arrays(alpha, beta)
+    return pt.full_like(alpha_b, pt.nan)
+
+
+def entropy(alpha, beta):
+    _, beta_b = pt.broadcast_arrays(alpha, beta)
+    return pt.log(4 * pt.pi * beta_b)
+
+
+def cdf(x, alpha, beta):
+    return pt.exp(logcdf(x, alpha, beta))
+
+
+def isf(x, alpha, beta):
+    return ppf(1 - x, alpha, beta)
+
+
+def pdf(x, alpha, beta):
+    return pt.exp(logpdf(x, alpha, beta))
+
+
+def ppf(q, alpha, beta):
+    return ppf_bounds_cont(alpha + beta * pt.tan(pt.pi * (q - 0.5)), q, -pt.inf, pt.inf)
+
+
+def sf(x, alpha, beta):
+    return pt.exp(logsf(x, alpha, beta))
+
+
+def rvs(alpha, beta, size=None, random_state=None):
+    return pt.random.cauchy(alpha, beta, rng=random_state, size=size)
+
+
+def logcdf(x, alpha, beta):
+    alpha_b, beta_b = pt.broadcast_arrays(alpha, beta)
+    return pt.log((1 / pt.pi) * pt.arctan((x - alpha_b) / beta_b) + 0.5)
+
+
+def logpdf(x, alpha, beta):
+    return -pt.log(pt.pi) - pt.log(beta) - pt.log(1 + ((x - alpha) / beta) ** 2)
+
+
+def logsf(x, alpha, beta):
+    return pt.log(0.5 - (1 / pt.pi) * pt.arctan((x - alpha) / beta))

distributions/halfcauchy.py

@@ -0,0 +1,80 @@
+import pytensor.tensor as pt
+
+from .helper import ppf_bounds_cont, cdf_bounds
+
+
+def mean(beta):
+    return pt.full_like(beta, pt.inf)
+
+
+def mode(beta):
+    return pt.zeros_like(beta)
+
+
+def median(beta):
+    return ppf(0.5, beta)
+
+
+def var(beta):
+    return pt.full_like(beta, pt.inf)
+
+
+def std(beta):
+    return pt.full_like(beta, pt.inf)
+
+
+def skewness(beta):
+    return pt.full_like(beta, pt.nan)
+
+
+def kurtosis(beta):
+    return pt.full_like(beta, pt.nan)
+
+
+def entropy(beta):
+    return pt.log(2 * pt.pi * beta)
+
+
+def cdf(x, beta):
+    return pt.exp(logcdf(x, beta))
+
+
+def logcdf(x, beta):
+    return pt.switch(
+        pt.lt(x, 0),
+        -pt.inf,
+        pt.log(2 * pt.arctan(x / beta) / pt.pi),
+    )
+
+def isf(x, beta):
+    return ppf(1 - x, beta)
+
+
+def pdf(x, beta):
+    return pt.exp(logpdf(x, beta))
+
+
+def ppf(q, beta):
+    x_val = beta * pt.tan(pt.pi / 2 * q)
+    return ppf_bounds_cont(x_val, q, 0, pt.inf)
+
+
+def sf(x, beta):
+    return 1 - cdf(x, beta)
+
+
+def rvs(beta, size=None, random_state=None):
+    uniform_samples = pt.random.uniform(0, 1, rng=random_state, size=size)
+    return beta * pt.tan(pt.pi / 2 * uniform_samples)
+
+
+def logpdf(x, beta):
+    return pt.where(
+        pt.lt(x, 0),
+        -pt.inf,
+        pt.log(2) - pt.log(pt.pi * beta) - pt.log(1 + (x / beta)**2)
+    )
+
+
+def logsf(x, beta):
+    return pt.log1mexp(logcdf(x, beta))

distributions/halfnormal.py

@@ -0,0 +1,74 @@
+import pytensor.tensor as pt
+
+from .helper import ppf_bounds_cont
+
+
+def mean(sigma):
+    return sigma * pt.sqrt(2 / pt.pi)
+
+
+def mode(sigma):
+    return pt.zeros_like(sigma)
+
+
+def median(sigma):
+    return sigma * pt.sqrt(2) * pt.erfinv(0.5)
+
+
+def var(sigma):
+    return sigma**2 * (1 - 2 / pt.pi)
+
+
+def std(sigma):
+    return sigma * pt.sqrt(1 - 2 / pt.pi)
+
+
+def skewness(sigma):
+    return pt.full_like(sigma, 0.9952717)
+
+
+def kurtosis(sigma):
+    return pt.full_like(sigma, 0.8691773)
+
+
+def entropy(sigma):
+    return 0.5 * pt.log(pt.pi * sigma**2 / 2) + 0.5
+
+
+def cdf(x, sigma):
+    return pt.where(pt.lt(x, 0), 0.0, pt.erf(x / (sigma * pt.sqrt(2))))
+
+
+def logcdf(x, sigma):
+    return pt.where(pt.lt(x, 0), -pt.inf, pt.log(pt.erf(x / (sigma * pt.sqrt(2)))))
+
+
+def isf(x, sigma):
+    return ppf(1 - x, sigma)
+
+
+def pdf(x, sigma):
+    return pt.exp(logpdf(x, sigma))
+
+
+def ppf(q, sigma):
+    x_vals = sigma * pt.sqrt(2) * pt.erfinv(q)
+    return ppf_bounds_cont(x_vals, q, 0, pt.inf)
+
+
+def sf(x, sigma):
+    return 1 - cdf(x, sigma)
+
+
+def rvs(sigma, size=None, random_state=None):
+    return pt.abs(pt.random.normal(0, sigma, rng=random_state, size=size))
+
+
+def logpdf(x, sigma):
+    return pt.where(
+        pt.lt(x, 0), -pt.inf, pt.log(pt.sqrt(2 / pt.pi)) - pt.log(sigma) - 0.5 * (x / sigma) ** 2
+    )
+
+
+def logsf(x, sigma):
+    return pt.log1mexp(logcdf(x, sigma))

distributions/halfstudentt.py

@@ -0,0 +1,135 @@
+import pytensor.tensor as pt
+from pytensor.tensor.special import betaln
+from pytensor.tensor.math import betaincinv
+
+from .helper import cdf_bounds, ppf_bounds_cont
+from .halfnormal import entropy as halfnormal_entropy
+from .halfnormal import cdf as halfnormal_cdf
+from .halfnormal import logcdf as halfnormal_logpdf
+
+
+
+def mean(nu, sigma):
+    gamma0 = pt.exp(pt.gammaln((nu + 1) / 2))
+    gamma1 = pt.exp(pt.gammaln(nu / 2))
+    gamma_nonfinite = pt.isinf(gamma0) | pt.isinf(gamma1)
+
+    mean_finite = 2 * sigma * pt.sqrt(nu / pt.pi) * (gamma0 / (gamma1 * (nu - 1)))
+
+    mean_approx = sigma * pt.sqrt(2 / pt.pi)
+    mean_val = pt.where(gamma_nonfinite, mean_approx, mean_finite)
+
+    return pt.where(nu > 1, mean_val, pt.inf)
+
+
+def mode(nu, sigma):
+    _, sigma_b = pt.broadcast_arrays(nu, sigma)
+    return pt.zeros_like(sigma_b)
+
+
+def median(nu, sigma):
+    return ppf(0.5, nu, sigma)
+
+
+def var(nu, sigma):
+    gamma0 = pt.exp(pt.gammaln((nu + 1) / 2))
+    gamma1 = pt.exp(pt.gammaln(nu / 2))
+    gamma_nonfinite = pt.isinf(gamma0) | pt.isinf(gamma1)
+
+    var_finite = sigma**2 * (
+        (nu / (nu - 2)) - ((4 * nu) / (pt.pi * (nu - 1) ** 2)) * (gamma0 / gamma1) ** 2
+    )
+
+    var_approx = sigma**2 * (1 - 2.0 / pt.pi)
+    var_val = pt.where(gamma_nonfinite, var_approx, var_finite)
+
+    return pt.where(nu > 2, var_val, pt.inf)
+
+
+def std(nu, sigma):
+    variance = var(nu, sigma)
+    return pt.switch(pt.isinf(variance) | pt.isnan(variance), variance, pt.sqrt(variance))
+
+
+def skewness(nu, sigma):
+    nu_b, _ = pt.broadcast_arrays(nu, sigma)
+    return pt.full_like(nu_b, pt.nan)
+
+
+def kurtosis(nu, sigma):
+    nu_b, _ = pt.broadcast_arrays(nu, sigma)
+    return pt.full_like(nu_b, pt.nan)
+
+
+def entropy(nu, sigma):
+    # we use a halfnormal approximation for large nu
+    return pt.switch(
+        pt.gt(nu, 1e5),
+        halfnormal_entropy(sigma),
+        pt.log(sigma)
+        + 0.5 * (nu + 1) * (pt.psi(0.5 * (nu + 1)) - pt.psi(0.5 * nu))
+        + pt.log(pt.sqrt(nu))
+        + betaln(0.5 * nu, 0.5),
+    ) - pt.log(2)
+
+
+def cdf(x, nu, sigma):
+    # we use a halfnormal approximation for large nu
+    x = x / sigma
+    factor = 0.5 * pt.betainc(0.5 * nu, 0.5, nu / (x**2 + nu))
+    cdf_ = pt.switch(pt.lt(x, 0), factor, 1 - factor) * 2 - 1
+    halft_logcdf = cdf_bounds(cdf_, x, 0, pt.inf)
+
+    return pt.switch(pt.gt(nu, 1e5), halfnormal_cdf(x, sigma), halft_logcdf)
+
+
+def isf(x, nu, sigma):
+    return ppf(1 - x, nu, sigma)
+
+
+def pdf(x, nu, sigma):
+    return pt.exp(logpdf(x, nu, sigma))
+
+
+def ppf(q, nu, sigma):
+    q_factor = (q + 1) / 2
+    inv_factor = pt.switch(
+        pt.lt(q_factor, 0.5),
+        betaincinv(0.5 * nu, 0.5, 2 * q_factor),
+        pt.sqrt(nu / betaincinv(0.5 * nu, 0.5, 2 - 2 * q_factor) - nu),
+    )
+    return ppf_bounds_cont(inv_factor * sigma, q, 0, pt.inf)
+
+
+def sf(x, nu, sigma):
+    return 1 - cdf(x, nu, sigma)
+
+
+def rvs(nu, sigma, size=None, random_state=None):
+    t_samples = pt.random.t(nu, rng=random_state, size=size)
+    return pt.abs(t_samples * sigma)
+
+
+def logcdf(x, nu, sigma):
+    return pt.log(cdf(x, nu, sigma))
+
+
+def logpdf(x, nu, sigma):
+    # we use a halfnormal approximation for large nu
+    halft_logpdf = pt.where(
+        pt.lt(x, 0),
+        -pt.inf,
+        (
+            pt.gammaln((nu + 1) / 2)
+            - pt.gammaln(nu / 2)
+            - 0.5 * pt.log(nu * pt.pi * sigma**2)
+            - 0.5 * (nu + 1) * pt.log(1 + (x / sigma) ** 2 / nu)
+            + pt.log(2)
+        ),
+    )
+
+    return pt.switch(pt.gt(nu, 1e5), halfnormal_logpdf(x, sigma), halft_logpdf)
+
+
+def logsf(x, nu, sigma):
+    return pt.log1mexp(logcdf(x, nu, sigma))

distributions/helper.py

@@ -1,10 +1,17 @@
 import pytensor.tensor as pt
 
 
-def ppf_bounds_cont(value, q, lower, upper):
-    q = pt.as_tensor_variable(q)
+def cdf_bounds(cdf_val, x, lower, upper):
     return pt.switch(
-        pt.and_(q >= lower, q <= upper),
-        value,
+        pt.or_(pt.lt(x, lower), pt.gt(x, upper)),
+        0.0,
+        pt.switch(pt.eq(x, lower), 0.0, pt.switch(pt.eq(x, upper), 1.0, cdf_val)),
+    )
+
+
+def ppf_bounds_cont(x_val, q, lower, upper):
+    return pt.switch(
+        pt.or_(pt.lt(q, 0), pt.gt(q, 1)),
         pt.nan,
-    )
+        pt.switch(pt.eq(q, 0), lower, pt.switch(pt.eq(q, 1), upper, x_val)),
+    )