Repository: scipy/scipy
Pull request: [ENH] ncx2 dispatch to chi2 when nc=0
- close #5441
- one test added

Problem statement:
scipy.stats.ncx2 fails for nc=0
For `nc = 0` it should reduce to the chi-squared distribution, but `scipy.stats.ncx2.rvs(1, 0)` gives:

``` python
C:\Anaconda\lib\site-packages\scipy\stats\_distn_infrastructure.pyc in rvs(self, *args, **kwds)
    882         cond = logical_and(self._argcheck(*args), (scale >= 0))
    883         if not np.all(cond):
--> 884             raise ValueError("Domain error in arguments.")
    885 
    886         # self._size is total size of all output values

ValueError: Domain error in arguments.
```

Similarly the pdf, etc. also fail.

Hints from issue discussion:
This specific error message comes from the fact that `ncx2` inherits a generic `_argcheck` which checks that all shape arguments are strictly positive. The fix is simple, just define `ncx2_gen._argcheck` to be `(df > 0) & (nc >=0)`.  

Once this is done, everything will keep failing because neither `np.random.noncentral_chisquare` (which `ncx2.rvs` defers to) nor formulas for pdf, cdf etc are written to handle zero non-centrality parameter. The easiest way of fixing it is probably to dispatch to chi2 methods for small enough nc:

```
In [20]: chi2.cdf(1, df) - ncx2.cdf(1, nc=1e-9, df=df)
Out[20]: 
array([  2.41974885e-10,   2.32778342e-14,   1.75990495e-20,
         1.26095641e-27,   1.84527014e-35,   8.11960902e-44,
         1.35776754e-52,   1.01017522e-61,   3.74443887e-71,
         7.53247129e-81])

In [21]: df
Out[21]: array([ 1, 11, 21, 31, 41, 51, 61, 71, 81, 91])
```

the dispatch should use `_lazywhere` to avoid division-by-zero warnings.

An alternative would be to patch up the `chndtr` and `chndtrix` functions. (which seem to have several accuracy issues anyway).

Would you be interested in working on this @roger-? 

Hi @ev-br , 

I'm new here and I wanted to contribute. I tried to understand rv_generic class, and then I defined a new method in the ncx2 class for argument checking as suggested by you. 
After that I passed non-centrality as zero. But this time I get a new error, pertaining to the numpy file mtrand.pyx. 

For `scipy.stats.ncx2.rvs(1, 0)` as given above:

Error: `File "mtrand.pyx", line 2162, in mtrand.RandomState.noncentral_chisquare (numpy/random/mtrand/mtrand.c:12806)
ValueError: df <= 0`

And for `scipy.stats.ncx2.rvs(2, 0)`:

Error: `File "mtrand.pyx", line 2164, in mtrand.RandomState.noncentral_chisquare (numpy/random/mtrand/mtrand.c:12830)
ValueError: nonc <= 0`

I am sorry for the naivety of my questions, I have invested very less time in understanding Scipy base, there are a lot of things that I need to learn.

Hi @ev-br and @Saxenauts ,

I just bumped into this bug as well, but it is not high priority for me / don't have time to fix it.

However, I think I have an easy fix. The original problem is that nc parameter appears in the denominator for standard expression used for this pdf. I thought this was quite weird as it seems there should be an equivalent formula which naturally simplifies to the chi2 when nc=1. Wikipedia to rescue - see last formula in the definition section https://en.wikipedia.org/wiki/Noncentral_chi-squared_distribution#Definition .
It's basically the chi2 formula but with the addition of the first two terms. This removes division by lambda, and the terms that depend on lambda behave fine even for nc=0 (the 0F1 function is available scipy.special.hyp0f1). So fixing the function to return the correct pdf value even for nc=0 should be a one-liner. I just don't know if this incurs some computational hit (and how important it is) as don't know how is hyp0f1 implemented.

Commit message: Add tests for pdf, logpdf, rvs

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 53ce63712ad8..d7af333d8a52 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -5204,19 +5204,25 @@ class ncx2_gen(rv_continuous):
     %(example)s
 
     """
+    def _argcheck(self, df, nc):
+        return (df > 0) & (nc >= 0)
+
     def _rvs(self, df, nc):
         return self._random_state.noncentral_chisquare(df, nc, self._size)
 
     def _logpdf(self, x, df, nc):
-        return _ncx2_log_pdf(x, df, nc)
+        cond = np.ones_like(x, dtype=bool) & (nc != 0)
+        return _lazywhere(cond, (x, df, nc), f=_ncx2_log_pdf, f2=chi2.logpdf)
 
     def _pdf(self, x, df, nc):
         # ncx2.pdf(x, df, nc) = exp(-(nc+x)/2) * 1/2 * (x/nc)**((df-2)/4)
         #                       * I[(df-2)/2](sqrt(nc*x))
-        return _ncx2_pdf(x, df, nc)
+        cond = np.ones_like(x, dtype=bool) & (nc != 0)
+        return _lazywhere(cond, (x, df, nc), f=_ncx2_pdf, f2=chi2.pdf)
 
     def _cdf(self, x, df, nc):
-        return _ncx2_cdf(x, df, nc)
+        cond = np.ones_like(x, dtype=bool) & (nc != 0)
+        return _lazywhere(cond, (x, df, nc), f=_ncx2_cdf, f2=chi2.cdf)
 
     def _ppf(self, q, df, nc):
         return sc.chndtrix(q, df, nc)
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 25e32b3240e8..ec83feb8a478 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -3258,6 +3258,26 @@ def test_ncx2_tails_pdf():
     assert_(np.isneginf(logval).all())
 
 
+@pytest.mark.parametrize('method, expected', [
+    ('cdf', np.array([1.72115630e-04, 2.74102033e-05])),
+    ('pdf', np.array([0.00078975, 0.00013013])),
+    ('logpdf', np.array([-7.14378973, -8.94700192]))
+])
+def test_ncx2_zero_nc(method, expected):
+    # gh-5441
+    # ncx2 with nc=0 is identical to chi2
+    result = getattr(stats.ncx2, method)(1, nc=[0, 4], df=10)
+    assert_allclose(result, expected, rtol=1e-3, atol=0)
+
+
+def test_ncx2_zero_nc_rvs():
+    # gh-5441
+    # ncx2 with nc=0 is identical to chi2
+    result = stats.ncx2.rvs(df=10, nc=[0,4], random_state=1)
+    expected = np.array([18.25726558,  8.24331569])
+    assert_allclose(result, expected, rtol=1e-3, atol=0)
+
+
 def test_foldnorm_zero():
     # Parameter value c=0 was not enabled, see gh-2399.
     rv = stats.foldnorm(0, scale=1)


Code review comments:
diff_hunk: @@ -3258,6 +3258,26 @@ def test_ncx2_tails_pdf():
     assert_(np.isneginf(logval).all())
 
 
+@pytest.mark.parametrize('method, expected', [
+    ('cdf', np.array([1.72115630e-04, 2.74102033e-05])),
+    ('pdf', np.array([0.00078975, 0.00013013])),
+    ('logpdf', np.array([-7.14378973, -8.94700192]))
+])
+def test_ncx2_zero_nc(method, expected):
+    # gh-5441
+    # ncx2 with nc=0 is identical to chi2
+    result = getattr(stats.ncx2, method)(1, nc=[0, 4], df=10)
path: scipy/stats/tests/test_distributions.py
text: @user1:
Can just parameterize with ncx2.pdf, ncx2.cdf, ncx2.logpdf and get rid of the getattr?

@author:
The parameter of the pytest decorator is regarded as `string` which doesn't have `__call__`.
Explicitly, if the parameters is set as follows:
```
@user2.mark.parametrize('method, expected', [
    ('ncx2.cdf', np.array([2.497951336e-09, 3.437288941e-10])),
    ('ncx2.pdf', np.array([1.238579980e-07, 1.710041145e-08])),
    ('ncx2.logpdf', np.array([-15.90413011, -17.88416331]))
])
```
`method(0.1, nc=[0, 4], df=10)` would run into `TypeError`.

diff_hunk: @@ -3258,6 +3258,26 @@ def test_ncx2_tails_pdf():
     assert_(np.isneginf(logval).all())
 
 
+@pytest.mark.parametrize('method, expected', [
+    ('cdf', np.array([1.72115630e-04, 2.74102033e-05])),
+    ('pdf', np.array([0.00078975, 0.00013013])),
+    ('logpdf', np.array([-7.14378973, -8.94700192]))
+])
+def test_ncx2_zero_nc(method, expected):
+    # gh-5441
+    # ncx2 with nc=0 is identical to chi2
+    result = getattr(stats.ncx2, method)(1, nc=[0, 4], df=10)
+    assert_allclose(result, expected, rtol=1e-3, atol=0)
path: scipy/stats/tests/test_distributions.py
text: @user1:
Why this combinbation of tolerances? 
Assuming that the numbers come from chi2 distribution for nc=0 and ncx2 for nc != 0, the tolerances could/should be made way tighter, close 1e-15 or some such.