MAINT: stats.pearson3: fix ppf for negative skew

scipy/stats/_continuous_distns.py

@@ -7178,7 +7178,9 @@ def _ppf(self, q, skew):
         ans, q, _, mask, invmask, beta, alpha, zeta = (
             self._preprocess(q, skew))
         ans[mask] = _norm_ppf(q[mask])
-        ans[invmask] = sc.gammaincinv(alpha, q[invmask])/beta + zeta
+        q = q[invmask]
+        q[beta < 0] = 1 - q[beta < 0]  # for negative skew; see gh-17050
+        ans[invmask] = sc.gammaincinv(alpha, q)/beta + zeta
         return ans
 
     @_call_super_mom

scipy/stats/_distr_params.py

@@ -87,6 +87,7 @@
     ['norminvgauss', (1.25, 0.5)],
     ['pareto', (2.621716532144454,)],
     ['pearson3', (0.1,)],
+    ['pearson3', (-2,)],
     ['powerlaw', (1.6591133289905851,)],
     ['powerlaw', (0.6591133289905851,)],
     ['powerlognorm', (2.1413923530064087, 0.44639540782048337)],

scipy/stats/tests/test_distributions.py

@@ -1834,6 +1834,25 @@ def test_return_array_bug_11746(self):
         assert_equal(moment, 0)
         assert isinstance(moment, np.number)
 
+    def test_ppf_bug_17050(self):
+        # incorrect PPF for negative skews were reported in gh-17050
+        # Check that this is fixed (even in the array case)
+        skews = [-3, -1, 0, 0.5]
+        x_eval = 0.5
+        res = stats.pearson3.ppf(stats.pearson3.cdf(x_eval, skews), skews)
+        assert_allclose(res, x_eval)
+
+        # Negation of the skew flips the distribution about the origin, so
+        # the following should hold
+        skew = np.array([[-0.5], [1.5]])
+        x = np.linspace(-2, 2)
+        assert_allclose(stats.pearson3.pdf(x, skew),
+                        stats.pearson3.pdf(-x, -skew))
+        assert_allclose(stats.pearson3.cdf(x, skew),
+                        stats.pearson3.sf(-x, -skew))
+        assert_allclose(stats.pearson3.ppf(x, skew),
+                        -stats.pearson3.isf(x, -skew))
+
 
 class TestKappa4:
     def test_cdf_genpareto(self):

scipy/stats/tests/test_fit.py

@@ -378,8 +378,8 @@ def test_basic_fit(self, dist_name):
         dist = getattr(stats, dist_name)
         shapes = np.array(dist_data[dist_name])
         bounds = np.empty((len(shapes) + 2, 2), dtype=np.float64)
-        bounds[:-2, 0] = shapes/10**np.sign(shapes)
-        bounds[:-2, 1] = shapes*10**np.sign(shapes)
+        bounds[:-2, 0] = shapes/10.**np.sign(shapes)
+        bounds[:-2, 1] = shapes*10.**np.sign(shapes)
         bounds[-2] = (0, 10)
         bounds[-1] = (0, 10)
         loc = rng.uniform(*bounds[-2])