Merge pull request scipy#6799 from thomaskeck/Argus

ENH: stats: Argus distribution

doc/source/tutorial/stats.rst

@@ -89,7 +89,7 @@ introspection:
     >>> dist_discrete = [d for d in dir(stats) if
     ...                  isinstance(getattr(stats,d), stats.rv_discrete)]
     >>> print 'number of continuous distributions:', len(dist_continu)
-    number of continuous distributions: 94
+    number of continuous distributions: 95
     >>> print 'number of discrete distributions:  ', len(dist_discrete)
     number of discrete distributions:   13
 

scipy/stats/__init__.py

@@ -26,6 +26,7 @@
    alpha             -- Alpha
    anglit            -- Anglit
    arcsine           -- Arcsine
+   argus             -- Argus
    beta              -- Beta
    betaprime         -- Beta Prime
    bradford          -- Bradford

scipy/stats/_continuous_distns.py

@@ -5206,6 +5206,63 @@ def _entropy(self, beta):
 halfgennorm = halfgennorm_gen(a=0, name='halfgennorm')
 
 
+def _argus_phi(chi):
+    """
+    Utility function for the argus distribution
+    used in the CDF and norm of the Argus Funktion
+    """
+    return _norm_cdf(chi) - chi * _norm_pdf(chi) - 0.5
+
+
+class argus_gen(rv_continuous):
+    """
+    Argus distribution
+
+    %(before_notes)s
+
+    Notes
+    -----
+    The probability density function for `argus` is::
+
+        argus.pdf(x, chi) = chi**3 / (sqrt(2*pi) * Psi(chi)) * x * sqrt(1-x**2) * exp(- 0.5 * chi**2 * (1 - x**2))
+
+        where:
+                 Psi(chi) = Phi(chi) - chi * phi(chi) - 1/2
+                 with Phi and phi being the CDF and PDF of a standard normal distribution, respectively.
+
+    `argus` takes ``chi`` as shape a parameter.
+
+    References
+    ----------
+
+    .. [1] "ARGUS distribution",
+           https://en.wikipedia.org/wiki/ARGUS_distribution
+
+    %(after_notes)s
+
+    %(example)s
+    """
+    def _pdf(self, x, chi):
+        """
+        Return PDF of the argus function
+        """
+        y = 1.0 - x**2
+        return chi**3 / (_norm_pdf_C * _argus_phi(chi)) * x * np.sqrt(y) * np.exp(-chi**2 * y / 2)
+
+    def _cdf(self, x, chi):
+        """
+        Return CDF of the argus function
+        """
+        return 1.0 - self._sf(x, chi)
+
+    def _sf(self, x, chi):
+        """
+        Return survival function of the argus function
+        """
+        return _argus_phi(chi * np.sqrt(1 - x**2)) / _argus_phi(chi)
+argus = argus_gen(name='argus', longname="An Argus Function", a=0.0, b=1.0)
+
+
 # Collect names of classes and objects in this module.
 pairs = list(globals().items())
 _distn_names, _distn_gen_names = get_distribution_names(pairs, rv_continuous)

scipy/stats/_distr_params.py

@@ -6,6 +6,7 @@
     ['alpha', (3.5704770516650459,)],
     ['anglit', ()],
     ['arcsine', ()],
+    ['argus', (1.0,)],
     ['beta', (2.3098496451481823, 0.62687954300963677)],
     ['betaprime', (5, 6)],
     ['bradford', (0.29891359763170633,)],

scipy/stats/tests/test_distributions.py

@@ -42,7 +42,7 @@
          'genlogistic', 'logistic', 'gumbel_l', 'gumbel_r', 'gompertz',
          'hypsecant', 'laplace', 'reciprocal', 'trapz', 'triang', 'tukeylambda',
          'vonmises', 'vonmises_line', 'pearson3', 'gennorm', 'halfgennorm',
-         'rice', 'kappa4', 'kappa3', 'truncnorm']
+         'rice', 'kappa4', 'kappa3', 'truncnorm', 'argus']
 
 
 def _assert_hasattr(a, b, msg=None):
@@ -3008,5 +3008,21 @@ def test_burr12_ppf_small_arg():
     assert_allclose(quantile, 5.7735026918962575e-09)
 
 
+def test_argus_function():
+    # There is no usable reference implementation.
+    # (RooFit implementation returns unreasonable results which are not normalized correctly)
+    # Instead we do some tests if the distribution behaves as expected for different shapes and scales
+    for i in range(1, 10):
+        for j in range(1, 10):
+            assert_equal(stats.argus.pdf(i + 0.001, chi=j, scale=i), 0.0)
+            assert_(stats.argus.pdf(i - 0.001, chi=j, scale=i) > 0.0)
+            assert_equal(stats.argus.pdf(-0.001, chi=j, scale=i), 0.0)
+            assert_(stats.argus.pdf(+0.001, chi=j, scale=i) > 0.0)
+
+    for i in range(1, 10):
+        assert_equal(stats.argus.cdf(1.0, chi=i), 1.0)
+        assert_equal(stats.argus.cdf(1.0, chi=i), 1.0 - stats.argus.sf(1.0, chi=i))
+
+
 if __name__ == "__main__":
     run_module_suite()