diff --git a/distributions/scaledinversechisquared.py b/distributions/scaledinversechisquared.py
new file mode 100644
index 0000000..9362ef8
--- /dev/null
+++ b/distributions/scaledinversechisquared.py
@@ -0,0 +1,90 @@
+import pytensor.tensor as pt
+
+from distributions.helper import cdf_bounds, ppf_bounds_cont
+
+
+def mean(nu, tau2):
+    return pt.switch(pt.gt(nu, 2), nu * tau2 / (nu - 2), pt.inf)
+
+
+def mode(nu, tau2):
+    return nu * tau2 / (nu + 2)
+
+
+def median(nu, tau2):
+    return ppf(0.5, nu, tau2)
+
+
+def var(nu, tau2):
+    return pt.switch(
+        pt.gt(nu, 4),
+        2 * nu**2 * tau2**2 / ((nu - 2) ** 2 * (nu - 4)),
+        pt.inf,
+    )
+
+
+def std(nu, tau2):
+    return pt.sqrt(var(nu, tau2))
+
+
+def skewness(nu, tau2):
+    return pt.switch(pt.gt(nu, 6), 4 / (nu - 6) * pt.sqrt(2 * (nu - 4)), pt.nan)
+
+
+def kurtosis(nu, tau2):
+    return pt.switch(pt.gt(nu, 8), (12 * (5 * nu - 22)) / ((nu - 6) * (nu - 8)), pt.nan)
+
+
+def entropy(nu, tau2):
+    h_nu = nu / 2
+    return h_nu + pt.log(h_nu * tau2) + pt.gammaln(h_nu) - (1 + h_nu) * pt.digamma(h_nu)
+
+
+def cdf(x, nu, tau2):
+    h_nu = nu / 2
+    return cdf_bounds(pt.gammaincc(h_nu, h_nu * tau2 / x), x, 0, pt.inf)
+
+
+def isf(x, nu, tau2):
+    return ppf(1 - x, nu, tau2)
+
+
+def pdf(x, nu, tau2):
+    return pt.exp(logpdf(x, nu, tau2))
+
+
+def ppf(q, nu, tau2):
+    h_nu = nu / 2
+    vals = h_nu * tau2 / pt.gammaincinv(h_nu, 1 - q)
+    return ppf_bounds_cont(vals, q, 0, pt.inf)
+
+
+def sf(x, nu, tau2):
+    return pt.exp(logsf(x, nu, tau2))
+
+
+def rvs(nu, tau2, size=None, random_state=None):
+    return (nu * tau2) / pt.random.chisquare(nu, rng=random_state, size=size)
+
+
+def logcdf(x, nu, tau2):
+    h_nu = nu / 2
+    return pt.switch(pt.le(x, 0), -pt.inf, pt.log(pt.gammaincc(h_nu, h_nu * tau2 / x)))
+
+
+def logpdf(x, nu, tau2):
+    h_nu = nu / 2
+    return pt.switch(
+        pt.le(x, 0),
+        -pt.inf,
+        -(pt.log(x) * (h_nu + 1))
+        - (h_nu * tau2) / x
+        + pt.log(tau2) * h_nu
+        - pt.gammaln(h_nu)
+        + pt.log(h_nu) * h_nu,
+    )
+
+
+def logsf(x, nu, tau2):
+    h_nu = nu / 2
+    return pt.switch(pt.le(x, 0), 0.0, pt.log(pt.gammainc(h_nu, h_nu * tau2 / x)))
diff --git a/tests/test_scaledinversechisquared.py b/tests/test_scaledinversechisquared.py
new file mode 100644
index 0000000..25a7c72
--- /dev/null
+++ b/tests/test_scaledinversechisquared.py
@@ -0,0 +1,32 @@
+"""Test ScaledInverseChiSquared distribution against scipy implementation."""
+
+import pytest
+from scipy import stats
+
+from distributions import scaledinversechisquared as ScaledInverseChiSquared
+from tests.helper_scipy import make_params, run_distribution_tests
+
+
+@pytest.mark.parametrize(
+    "params, sp_params",
+    [
+        ([4.0, 1.0], {"a": 2.0, "scale": 2.0}),
+        ([10.0, 2.0], {"a": 5.0, "scale": 10.0}),
+        ([2.5, 0.5], {"a": 1.25, "scale": 0.625}),
+        ([8.0, 3.0], {"a": 4.0, "scale": 12.0}),
+    ],
+)
+def test_scaledinversechisquared_vs_scipy(params, sp_params):
+    """Test ScaledInverseChiSquared distribution against scipy inverse gamma."""
+    p_params = make_params(*params, dtype="float64")
+    support = (0, float("inf"))
+
+    run_distribution_tests(
+        p_dist=ScaledInverseChiSquared,
+        sp_dist=stats.invgamma,
+        p_params=p_params,
+        sp_params=sp_params,
+        support=support,
+        name="scaledinversechisquared",
+        use_quantiles_for_rvs=True,
+    )