diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index 0980356432..12d93bab10 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -3727,6 +3727,13 @@ def dist(cls, mu=0, sigma=1.0, *args, **kwargs):
 
         return super().dist([mu, sigma], *args, **kwargs)
 
+    def get_moment(rv, size, mu, sigma):
+        mean = mu + sigma * (np.euler_gamma + at.log(2))
+
+        if not rv_size_is_none(size):
+            mean = at.full(size, mean)
+        return mean
+
     def logp(value, mu, sigma):
         """
         Calculate log-probability of Moyal distribution at specified value.
diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py
index f52ca8baa0..2572f3adf5 100644
--- a/pymc/tests/test_distributions_moments.py
+++ b/pymc/tests/test_distributions_moments.py
@@ -25,6 +25,7 @@
     Laplace,
     Logistic,
     LogNormal,
+    Moyal,
     NegativeBinomial,
     Normal,
     Pareto,
@@ -612,3 +613,18 @@ def test_discrete_uniform_moment(lower, upper, size, expected):
     with Model() as model:
         DiscreteUniform("x", lower=lower, upper=upper, size=size)
     assert_moment_is_expected(model, expected)
+
+
+@pytest.mark.parametrize(
+    "mu, sigma, size, expected",
+    [
+        (4.0, 3.0, None, 7.8110885363844345),
+        (4, np.full(5, 3), None, np.full(5, 7.8110885363844345)),
+        (np.arange(5), 1, None, np.arange(5) + 1.2703628454614782),
+        (np.arange(5), np.ones(5), (2, 5), np.full((2, 5), np.arange(5) + 1.2703628454614782)),
+    ],
+)
+def test_moyal_moment(mu, sigma, size, expected):
+    with Model() as model:
+        Moyal("x", mu=mu, sigma=sigma, size=size)
+    assert_moment_is_expected(model, expected)