diff --git a/pytensor_distributions/scaledinversechisquared.py b/pytensor_distributions/scaledinversechisquared.py
deleted file mode 100644
index dc08c83..0000000
--- a/pytensor_distributions/scaledinversechisquared.py
+++ /dev/null
@@ -1,109 +0,0 @@
-import pytensor.tensor as pt
-
-from pytensor_distributions.helper import cdf_bounds, ppf_bounds_cont
-from pytensor_distributions.lmoments import _lmoments
-
-
-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 lmoment1(nu, tau2):
-    return mean(nu, tau2)
-
-
-def lmoment2(nu, tau2):
-    return pt.switch(pt.gt(nu, 4), _lmoments(ppf, nu, tau2, r=2), pt.inf)
-
-
-def lmoment3(nu, tau2):
-    return pt.switch(pt.gt(nu, 6), _lmoments(ppf, nu, tau2, r=3), pt.nan)
-
-
-def lmoment4(nu, tau2):
-    return pt.switch(pt.gt(nu, 8), _lmoments(ppf, nu, tau2, r=4), 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, return_next_rng=True)[
-        1
-    ]
-
-
-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/helper_scipy.py b/tests/helper_scipy.py
index 30ffa21..d52745f 100644
--- a/tests/helper_scipy.py
+++ b/tests/helper_scipy.py
@@ -288,16 +288,6 @@ def run_lmoments_test(p_dist, p_params, name=None, rtol=5e-2, atol=5e-2, sample_
             if p_params[1].eval() <= 1:
                 s_l2 = s_tau3 = s_tau4 = float("inf")
 
-        if name == "scaledinversechisquared":
-            param = p_params[0].eval()
-            if param <= 4:
-                s_l2 = float("inf")
-                s_tau3 = s_tau4 = float("nan")
-            elif param <= 6:
-                s_tau3 = s_tau4 = float("nan")
-            elif param <= 8:
-                s_tau4 = float("nan")
-
         if name == "vonmises":
             s_tau3 = 0.0
 
diff --git a/tests/test_scaledinversechisquared.py b/tests/test_scaledinversechisquared.py
deleted file mode 100644
index 68057c8..0000000
--- a/tests/test_scaledinversechisquared.py
+++ /dev/null
@@ -1,32 +0,0 @@
-"""Test ScaledInverseChiSquared distribution against scipy implementation."""
-
-import pytest
-from scipy import stats
-
-from pytensor_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}),
-        ([3.5, 0.5], {"a": 1.75, "scale": 0.875}),
-        ([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,
-    )