Merge pull request scipy#21335 from WarrenWeckesser/boost-cauchy

BUG/ENH: stats: updates for cauchy.

Author: dschmitz89

scipy/special/_add_newdocs.py

@@ -11922,6 +11922,50 @@ def add_newdoc(name, doc):
 
     """)
 
+add_newdoc(
+    "_cauchy_ppf",
+    """
+    _cauchy_ppf(p, loc, scale)
+
+    Percent point function (i.e. quantile) of the Cauchy distribution.
+
+    Parameters
+    ----------
+    p : array_like
+        Probabilities
+    loc : array_like
+        Location parameter of the distribution.
+    scale : array_like
+        Scale parameter of the distribution.
+
+    Returns
+    -------
+    scalar or ndarray
+
+    """)
+
+add_newdoc(
+    "_cauchy_isf",
+    """
+    _cauchy_isf(p, loc, scale)
+
+    Inverse survival function of the Cauchy distribution.
+
+    Parameters
+    ----------
+    p : array_like
+        Probabilities
+    loc : array_like
+        Location parameter of the distribution.
+    scale : array_like
+        Scale parameter of the distribution.
+
+    Returns
+    -------
+    scalar or ndarray
+
+    """)
+
 add_newdoc(
     "_ncx2_pdf",
     """

scipy/special/boost_special_functions.h

@@ -588,6 +588,46 @@ invgauss_isf_double(double x, double mu, double s)
     return invgauss_isf_wrap(x, mu, s);
 }
 
+template<typename Real>
+Real
+cauchy_ppf_wrap(const Real p, const Real loc, const Real scale)
+{
+    return boost::math::quantile(
+        boost::math::cauchy_distribution<Real, StatsPolicy>(loc, scale), p);
+}
+
+float
+cauchy_ppf_float(float p, float loc, float scale)
+{
+    return cauchy_ppf_wrap(p, loc, scale);
+}
+
+double
+cauchy_ppf_double(double p, double loc, double scale)
+{
+    return cauchy_ppf_wrap(p, loc, scale);
+}
+
+template<typename Real>
+Real
+cauchy_isf_wrap(const Real p, const Real loc, const Real scale)
+{
+    return boost::math::quantile(boost::math::complement(
+        boost::math::cauchy_distribution<Real, StatsPolicy>(loc, scale), p));
+}
+
+float
+cauchy_isf_float(float p, float loc, float scale)
+{
+    return cauchy_isf_wrap(p, loc, scale);
+}
+
+double
+cauchy_isf_double(double p, double loc, double scale)
+{
+    return cauchy_isf_wrap(p, loc, scale);
+}
+
 template<typename Real>
 Real
 ncx2_pdf_wrap(const Real x, const Real k, const Real l)

scipy/special/functions.json

@@ -1033,6 +1033,18 @@
             "invgauss_isf_double": "ddd->d"
         }
     },
+    "_cauchy_ppf": {
+        "boost_special_functions.h++": {
+            "cauchy_ppf_float": "fff->f",
+            "cauchy_ppf_double": "ddd->d"
+        }
+    },
+    "_cauchy_isf": {
+        "boost_special_functions.h++": {
+            "cauchy_isf_float": "fff->f",
+            "cauchy_isf_double": "ddd->d"
+        }
+    },
     "_ncx2_pdf": {
         "boost_special_functions.h++": {
             "ncx2_pdf_float": "fff->f",

scipy/stats/_continuous_distns.py

@@ -1454,16 +1454,16 @@ def _pdf(self, x):
         return 1.0/np.pi/(1.0+x*x)
 
     def _cdf(self, x):
-        return 0.5 + 1.0/np.pi*np.arctan(x)
+        return np.arctan2(1, -x)/np.pi
 
     def _ppf(self, q):
-        return np.tan(np.pi*q-np.pi/2.0)
+        return scu._cauchy_ppf(q, 0, 1)
 
     def _sf(self, x):
-        return 0.5 - 1.0/np.pi*np.arctan(x)
+        return np.arctan2(1, x)/np.pi
 
     def _isf(self, q):
-        return np.tan(np.pi/2.0-np.pi*q)
+        return scu._cauchy_isf(q, 0, 1)
 
     def _stats(self):
         return np.nan, np.nan, np.nan, np.nan

scipy/stats/tests/test_continuous_basic.py

@@ -75,7 +75,7 @@
 # Here 'fail' mean produce wrong results and/or raise exceptions, depending
 # on the implementation details of corresponding special functions.
 # cf https://github.com/scipy/scipy/pull/4979 for a discussion.
-fails_cmplx = {'argus', 'beta', 'betaprime', 'chi', 'chi2', 'cosine',
+fails_cmplx = {'argus', 'beta', 'betaprime', 'cauchy', 'chi', 'chi2', 'cosine',
                'dgamma', 'dweibull', 'erlang', 'f', 'foldcauchy', 'gamma',
                'gausshyper', 'gengamma', 'genhyperbolic',
                'geninvgauss', 'gennorm', 'genpareto',

scipy/stats/tests/test_distributions.py

@@ -507,6 +507,48 @@ def test_cdf_ppf(self):
         assert_allclose(x, xx)
 
 
+class TestCauchy:
+
+    # Reference values were computed with mpmath.
+    @pytest.mark.parametrize(
+        'x, ref',
+        [(-5e15, 6.366197723675814e-17),
+         (-5, 0.06283295818900118),
+         (-1, 0.25),
+         (0, 0.5),
+         (1, 0.75),
+         (5, 0.9371670418109989),
+         (5e15, 0.9999999999999999)]
+    )
+    @pytest.mark.parametrize(
+        'method, sgn',
+        [(stats.cauchy.cdf, 1),
+         (stats.cauchy.sf, -1)]
+    )
+    def test_cdf_sf(self, x, ref, method, sgn):
+        p = method(sgn*x)
+        assert_allclose(p, ref, rtol=1e-15)
+
+    # Reference values were computed with mpmath.
+    @pytest.mark.parametrize(
+        'p, ref',
+        [(1e-20, -3.1830988618379067e+19),
+         (1e-9, -318309886.1837906),
+         (0.25, -1.0),
+         (0.50, 0.0),
+         (0.75, 1.0),
+         (0.999999, 318309.88617359026),
+         (0.999999999999, 318316927901.77966)]
+    )
+    @pytest.mark.parametrize(
+        'method, sgn',
+        [(stats.cauchy.ppf, 1),
+         (stats.cauchy.isf, -1)])
+    def test_ppf_isf(self, p, ref, method, sgn):
+        x = sgn*method(p)
+        assert_allclose(x, ref, rtol=1e-15)
+
+
 class TestChi:
 
     # "Exact" value of chi.sf(10, 4), as computed by Wolfram Alpha with