diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 95186e6b31c9..61d16193470f 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -4067,6 +4067,9 @@ def _cdf(self, x):
         # Equivalent to 2*norm.sf(np.sqrt(1/x))
         return sc.erfc(np.sqrt(0.5 / x))
 
+    def _sf(self, x):
+        return sc.erf(np.sqrt(0.5 / x))
+
     def _ppf(self, q):
         # Equivalent to 1.0/(norm.isf(q/2)**2) or 0.5/(erfcinv(q)**2)
         val = -sc.ndtri(q/2)
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 4453993b6aed..def499655d69 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -1797,13 +1797,6 @@ def test_cdf_ppf_next(self):
         assert_array_equal(rv.ppf(rv.cdf(rv.xk[:-1]) + 1e-8),
                            rv.xk[1:])
 
-    def test_expect(self):
-        xk = [1, 2, 4, 6, 7, 11]
-        pk = [0.1, 0.2, 0.2, 0.2, 0.2, 0.1]
-        rv = stats.rv_discrete(values=(xk, pk))
-
-        assert_allclose(rv.expect(), np.sum(rv.xk * rv.pk), atol=1e-14)
-
     def test_multidimension(self):
         xk = np.arange(12).reshape((3, 4))
         pk = np.array([[0.1, 0.1, 0.15, 0.05],
@@ -1843,8 +1836,15 @@ def test_shape_rv_sample(self):
         # same shapes => no error
         xk, pk = np.arange(6).reshape((3, 2)), np.full((3, 2), 1/6)
         assert_equal(stats.rv_discrete(values=(xk, pk)).pmf(0), 1/6)
-        
-    def test_expect(self):
+
+    def test_expect1(self):
+        xk = [1, 2, 4, 6, 7, 11]
+        pk = [0.1, 0.2, 0.2, 0.2, 0.2, 0.1]
+        rv = stats.rv_discrete(values=(xk, pk))
+
+        assert_allclose(rv.expect(), np.sum(rv.xk * rv.pk), atol=1e-14)
+
+    def test_expect2(self):
         # rv_sample should override _expect. Bug report from
         # https://stackoverflow.com/questions/63199792
         y = [200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0,
@@ -1857,13 +1857,14 @@ def test_expect(self):
         py = [0.0004, 0.0, 0.0033, 0.006500000000000001, 0.0, 0.0,
               0.004399999999999999, 0.6862, 0.0, 0.0, 0.0,
               0.00019999999999997797, 0.0006000000000000449,
-              0.024499999999999966, 0.006400000000000072, 0.0043999999999999595,
-              0.019499999999999962, 0.03770000000000007, 0.01759999999999995,
-              0.015199999999999991, 0.018100000000000005, 0.04500000000000004,
-              0.0025999999999999357, 0.0, 0.0041000000000001036,
-              0.005999999999999894, 0.0042000000000000925,
-              0.0050000000000000044, 0.0041999999999999815,
-              0.0004999999999999449, 0.009199999999999986, 0.008200000000000096,
+              0.024499999999999966, 0.006400000000000072,
+              0.0043999999999999595, 0.019499999999999962,
+              0.03770000000000007, 0.01759999999999995, 0.015199999999999991,
+              0.018100000000000005, 0.04500000000000004, 0.0025999999999999357,
+              0.0, 0.0041000000000001036, 0.005999999999999894,
+              0.0042000000000000925, 0.0050000000000000044,
+              0.0041999999999999815, 0.0004999999999999449,
+              0.009199999999999986, 0.008200000000000096,
               0.0, 0.0, 0.0046999999999999265, 0.0019000000000000128,
               0.0006000000000000449, 0.02510000000000001, 0.0,
               0.007199999999999984, 0.0, 0.012699999999999934, 0.0, 0.0,
@@ -3935,6 +3936,18 @@ def test_levy_cdf_ppf():
     assert_allclose(xx, x, rtol=1e-13)
 
 
+def test_levy_sf():
+    # Large values, far into the tail of the distribution.
+    x = np.array([1e15, 1e25, 1e35, 1e50])
+    # Expected values were calculated with mpmath.
+    expected = np.array([2.5231325220201597e-08,
+                         2.52313252202016e-13,
+                         2.52313252202016e-18,
+                         7.978845608028653e-26])
+    y = stats.levy.sf(x)
+    assert_allclose(y, expected, rtol=1e-14)
+
+
 def test_hypergeom_interval_1802():
     # these two had endless loops
     assert_equal(stats.hypergeom.interval(.95, 187601, 43192, 757),