diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py
index d2d12fc1e5..76c1ad794f 100644
--- a/pymc3/distributions/continuous.py
+++ b/pymc3/distributions/continuous.py
@@ -3959,6 +3959,21 @@ def _distr_parameters_for_repr(self):
         return []
 
 
+class MoyalRV(RandomVariable):
+    name = "moyal"
+    ndim_supp = 0
+    ndims_params = [0, 0]
+    dtype = "floatX"
+    _print_name = ("Moyal", "\\operatorname{Moyal}")
+
+    @classmethod
+    def rng_fn(cls, rng, mu, sigma, size=None):
+        return stats.moyal.rvs(mu, sigma, size=size, random_state=rng)
+
+
+moyal = MoyalRV()
+
+
 class Moyal(Continuous):
     r"""
     Moyal log-likelihood.
@@ -4006,43 +4021,18 @@ class Moyal(Continuous):
     sigma: float
         Scale parameter (sigma > 0).
     """
+    rv_op = moyal
 
-    def __init__(self, mu=0, sigma=1.0, *args, **kwargs):
-        self.mu = at.as_tensor_variable(floatX(mu))
-        self.sigma = at.as_tensor_variable(floatX(sigma))
+    @classmethod
+    def dist(cls, mu=0, sigma=1.0, *args, **kwargs):
+        mu = at.as_tensor_variable(floatX(mu))
+        sigma = at.as_tensor_variable(floatX(sigma))
 
         assert_negative_support(sigma, "sigma", "Moyal")
 
-        self.mean = self.mu + self.sigma * (np.euler_gamma + at.log(2))
-        self.median = self.mu - self.sigma * at.log(2 * at.erfcinv(1 / 2) ** 2)
-        self.mode = self.mu
-        self.variance = (np.pi ** 2 / 2.0) * self.sigma ** 2
-
-        super().__init__(*args, **kwargs)
-
-    def random(self, point=None, size=None):
-        """
-        Draw random values from Moyal distribution.
-
-        Parameters
-        ----------
-        point: dict, optional
-            Dict of variable values on which random values are to be
-            conditioned (uses default point if not specified).
-        size: int, optional
-            Desired size of random sample (returns one sample if not
-            specified).
-
-        Returns
-        -------
-        array
-        """
-        # mu, sigma = draw_values([self.mu, self.sigma], point=point, size=size)
-        # return generate_samples(
-        #     stats.moyal.rvs, loc=mu, scale=sigma, dist_shape=self.shape, size=size
-        # )
+        return super().dist([mu, sigma], *args, **kwargs)
 
-    def logp(self, value):
+    def logp(value, mu, sigma):
         """
         Calculate log-probability of Moyal distribution at specified value.
 
@@ -4056,15 +4046,13 @@ def logp(self, value):
         -------
         TensorVariable
         """
-        mu = self.mu
-        sigma = self.sigma
         scaled = (value - mu) / sigma
         return bound(
             (-(1 / 2) * (scaled + at.exp(-scaled)) - at.log(sigma) - (1 / 2) * at.log(2 * np.pi)),
             0 < sigma,
         )
 
-    def logcdf(self, value):
+    def logcdf(value, mu, sigma):
         """
         Compute the log of the cumulative distribution function for Moyal distribution
         at the specified value.
@@ -4079,9 +4067,6 @@ def logcdf(self, value):
         -------
         TensorVariable
         """
-        mu = self.mu
-        sigma = self.sigma
-
         scaled = (value - mu) / sigma
         return bound(
             at.log(at.erfc(at.exp(-scaled / 2) * (2 ** -0.5))),
diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py
index a756281fdd..bad3670df0 100644
--- a/pymc3/tests/test_distributions.py
+++ b/pymc3/tests/test_distributions.py
@@ -2547,7 +2547,6 @@ def test_rice(self):
             lambda value, b, sigma: sp.rice.logpdf(value, b=b, loc=0, scale=sigma),
         )
 
-    @pytest.mark.xfail(reason="Distribution not refactored yet")
     def test_moyal_logp(self):
         # Using a custom domain, because the standard `R` domain undeflows with scipy in float64
         value_domain = Domain([-inf, -1.5, -1, -0.01, 0.0, 0.01, 1, 1.5, inf])
@@ -2558,7 +2557,6 @@ def test_moyal_logp(self):
             lambda value, mu, sigma: floatX(sp.moyal.logpdf(value, mu, sigma)),
         )
 
-    @pytest.mark.xfail(reason="Distribution not refactored yet")
     @pytest.mark.xfail(
         condition=(aesara.config.floatX == "float32"),
         reason="Pymc3 underflows earlier than scipy on float32",
diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py
index 2759f2d3b5..0bf820df8b 100644
--- a/pymc3/tests/test_distributions_random.py
+++ b/pymc3/tests/test_distributions_random.py
@@ -319,12 +319,6 @@ class TestZeroInflatedBinomial(BaseTestCases.BaseTestCase):
     params = {"n": 10, "p": 0.6, "psi": 0.3}
 
 
-@pytest.mark.xfail(reason="This distribution has not been refactored for v4")
-class TestMoyal(BaseTestCases.BaseTestCase):
-    distribution = pm.Moyal
-    params = {"mu": 0.0, "sigma": 1.0}
-
-
 class BaseTestDistribution(SeededTest):
     pymc_dist: Optional[Callable] = None
     pymc_dist_params = dict()
@@ -491,6 +485,19 @@ class TestStudentT(BaseTestDistribution):
     ]
 
 
+class TestMoyal(BaseTestDistribution):
+    pymc_dist = pm.Moyal
+    pymc_dist_params = {"mu": 0.0, "sigma": 1.0}
+    expected_rv_op_params = {"mu": 0.0, "sigma": 1.0}
+    reference_dist_params = {"loc": 0.0, "scale": 1.0}
+    reference_dist = seeded_scipy_distribution_builder("moyal")
+    tests_to_run = [
+        "check_pymc_params_match_rv_op",
+        "check_pymc_draws_match_reference",
+        "check_rv_size",
+    ]
+
+
 class TestStudentTLam(BaseTestDistribution):
     pymc_dist = pm.StudentT
     lam, sigma = get_tau_sigma(tau=2.0)
@@ -1391,7 +1398,6 @@ def ref_rand(size, mu, sigma):
 
         pymc3_random(pm.LogitNormal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand)
 
-    @pytest.mark.xfail(reason="This distribution has not been refactored for v4")
     def test_moyal(self):
         def ref_rand(size, mu, sigma):
             return st.moyal.rvs(loc=mu, scale=sigma, size=size)