BUG: stats: Improve genpareto stats calculations.

Implement the _stats method for genpareto, so we can
use the explicit formulas for the mean, variance,
skewness and excess kurtosis.

Closes scipygh-11168.

scipy/stats/_continuous_distns.py

@@ -2424,6 +2424,35 @@ def _ppf(self, q, c):
     def _isf(self, q, c):
         return -sc.boxcox(q, -c)
 
+    def _stats(self, c, moments='mv'):
+        if 'm' not in moments:
+            m = None
+        else:
+            m = _lazywhere(c < 1, (c,),
+                           lambda xi: 1/(1 - xi),
+                           np.inf)
+        if 'v' not in moments:
+            v = None
+        else:
+            v = _lazywhere(c < 1/2, (c,),
+                           lambda xi: 1 / (1 - xi)**2 / (1 - 2*xi),
+                           np.nan)
+        if 's' not in moments:
+            s = None
+        else:
+            s = _lazywhere(c < 1/3, (c,),
+                           lambda xi: 2 * (1 + xi) * np.sqrt(1 - 2*xi) /
+                                      (1 - 3*xi),
+                           np.nan)
+        if 'k' not in moments:
+            k = None
+        else:
+            k = _lazywhere(c < 1/4, (c,),
+                           lambda xi: 3 * (1 - 2*xi) * (2*xi**2 + xi + 3) /
+                                      (1 - 3*xi) / (1 - 4*xi) - 3,
+                           np.nan)
+        return m, v, s, k
+
     def _munp(self, n, c):
         def __munp(n, c):
             val = 0.0

scipy/stats/tests/test_distributions.py

@@ -1007,6 +1007,23 @@ def test_logsf(self):
         logp = stats.genpareto.logsf(1e10, .01, 0, 1)
         assert_allclose(logp, -1842.0680753952365)
 
+    # Values in 'expected_stats' are
+    # [mean, variance, skewness, excess kurtosis].
+    @pytest.mark.parametrize(
+        'c, expected_stats',
+        [(0, [1, 1, 2, 6]),
+         (1/4, [4/3, 32/9, 10/np.sqrt(2), np.nan]),
+         (1/9, [9/8, (81/64)*(9/7), (10/9)*np.sqrt(7), 754/45]),
+         (-1, [1/2, 1/12, 0, -6/5])])
+    def test_stats(self, c, expected_stats):
+        result = stats.genpareto.stats(c, moments='mvsk')
+        assert_allclose(result, expected_stats, rtol=1e-13, atol=1e-15)
+
+    def test_var(self):
+        # Regression test for gh-11168.
+        v = stats.genpareto.var(1e-8)
+        assert_allclose(v, 1.000000040000001, rtol=1e-13)
+
 
 class TestPearson3(object):
     def setup_method(self):