Repository: scipy/scipy
File: scipy/stats/tests/test_continuous.py (Python)

Code before:
                        'sigma': {'endpoints': (0, np.inf), 'inclusive': (False, False)}}

            @property
            def support(self):
                def left(*, c, mu, sigma):
                    c, mu, sigma = np.broadcast_arrays(c, mu, sigma)
                    result = np.empty_like(c)
                    result[c >= 0] = -np.inf
                    result[c < 0] = mu[c < 0] + sigma[c < 0] / c[c < 0]
                    return result[()]

                def right(*, c, mu, sigma):
                    c, mu, sigma = np.broadcast_arrays(c, mu, sigma)
                    result = np.empty_like(c)
                    result[c <= 0] = np.inf
                    result[c > 0] = mu[c > 0] + sigma[c > 0] / c[c > 0]
                    return result[()]

                return {"endpoints": (left, right), "inclusive": (False, False)}

            def pdf(self, x, *, c, mu, sigma):
                x, c, mu, sigma = np.broadcast_arrays(x, c, mu, sigma)
                t = np.empty_like(x)
                mask = (c == 0)
                t[mask] = np.exp(-(x[mask] - mu[mask])/sigma[mask])
                t[~mask] = (1  - c[~mask]*(x[~mask] - mu[~mask])/sigma[~mask])**(1/c[~mask])
                result = 1/sigma * t**(1 - c)*np.exp(-t)
                return result[()]

            def cdf(self, x, *, c, mu, sigma):
                x, c, mu, sigma = np.broadcast_arrays(x, c, mu, sigma)
                t = np.empty_like(x)
                mask = (c == 0)
                t[mask] = np.exp(-(x[mask] - mu[mask])/sigma[mask])
                t[~mask] = (1  - c[~mask]*(x[~mask] - mu[~mask])/sigma[~mask])**(1/c[~mask])
                return np.exp(-t)[()]

        GenExtreme1 = stats.make_distribution(MyGenExtreme())
        GenExtreme2 = stats.make_distribution(stats.genextreme)

        X1 = GenExtreme1(c=c, mu=0, sigma=1)
        X2 = GenExtreme2(c=c)

        x = X1.sample(shape=10, rng=rng)
        p = X1.cdf(x)

        assert_allclose(X1.support(), X2.support())
        assert_allclose(X1.entropy(), X2.entropy(), rtol=5e-6)
        assert_allclose(X1.median(), X2.median())
        assert_allclose(X1.logpdf(x), X2.logpdf(x))
        assert_allclose(X1.pdf(x), X2.pdf(x))

Review comment (suggestion):
Not needed now because `_lazywhere` is being transitioned to `xpx.apply_where`, but just to show what it would look like:
```suggestion
                t = _lazywhere(c == 0, (x, c, mu, sigma),
                               lambda x, c, mu, sigma: np.exp(-(x - mu)/sigma),
                               lambda x, c, mu, sigma: (1 - c*(x - mu)/sigma)**(1 / c))
```
A bit easier on the eyes, and I think it does the broadcasting for you.
But maybe it's because I was looking at `apply_where` PR for an hour and got used to it.

Code after:
                    'mu': {'endpoints': (-np.inf, np.inf), 'inclusive': (False, False)},
                    'sigma': {'endpoints': (0, np.inf), 'inclusive': (False, False)}
                }

            @property
            def support(self):
                def left(*, c, mu, sigma):
                    c, mu, sigma = np.broadcast_arrays(c, mu, sigma)
                    result = np.empty_like(c)
                    result[c >= 0] = -np.inf
                    result[c < 0] = mu[c < 0] + sigma[c < 0] / c[c < 0]
                    return result[()]

                def right(*, c, mu, sigma):
                    c, mu, sigma = np.broadcast_arrays(c, mu, sigma)
                    result = np.empty_like(c)
                    result[c <= 0] = np.inf
                    result[c > 0] = mu[c > 0] + sigma[c > 0] / c[c > 0]
                    return result[()]

                return {"endpoints": (left, right), "inclusive": (False, False)}

            def pdf(self, x, *, c, mu, sigma):
                x, c, mu, sigma = np.broadcast_arrays(x, c, mu, sigma)
                t = np.empty_like(x)
                mask = (c == 0)
                t[mask] = np.exp(-(x[mask] - mu[mask])/sigma[mask])
                t[~mask] = (
                    1  - c[~mask]*(x[~mask] - mu[~mask])/sigma[~mask]
                )**(1/c[~mask])
                result = 1/sigma * t**(1 - c)*np.exp(-t)
                return result[()]

            def cdf(self, x, *, c, mu, sigma):
                x, c, mu, sigma = np.broadcast_arrays(x, c, mu, sigma)
                t = np.empty_like(x)
                mask = (c == 0)
                t[mask] = np.exp(-(x[mask] - mu[mask])/sigma[mask])
                t[~mask] = (
                    1  - c[~mask]*(x[~mask] - mu[~mask])/sigma[~mask]
                )**(1/c[~mask])
                return np.exp(-t)[()]

        GenExtreme1 = stats.make_distribution(MyGenExtreme())
        GenExtreme2 = stats.make_distribution(stats.genextreme)

        X1 = GenExtreme1(c=c, mu=0, sigma=1)
        X2 = GenExtreme2(c=c)

        x = X1.sample(shape=10, rng=rng)
        p = X1.cdf(x)