Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

Rice distribution: extend to b=0, add an explicit rvs method. #2847

Merged
merged 3 commits into from

3 participants

@ev-br
Collaborator

closes gh-2164 and gh-2742.

As a bonus, rvs get a bit faster. With current master:

In [10]: scipy.__version__
Out[10]: '0.14.0.dev-efbde60'

In [11]: %timeit stats.rice.rvs(b=3.)
100 loops, best of 3: 4.18 ms per loop

In [12]: %timeit stats.rice.rvs(b=3., size=1000)
1 loops, best of 3: 4.13 s per loop

with this PR:

In [6]: scipy.__version__
Out[6]: '0.14.0.dev-bf7d8d9'

In [7]: %timeit stats.rice.rvs(b=3.)
10000 loops, best of 3: 45.5 µs per loop

In [8]: %timeit stats.rice.rvs(b=3., size=1000)
10000 loops, best of 3: 192 µs per loop
@josef-pkt
Collaborator

nice, ticket count is going down

looks good to me

scipy/stats/tests/test_continuous_extra.py
((7 lines not shown))
+ npt.assert_(np.isfinite(stats.rice.pdf(x, b=0.)).all())
+ npt.assert_(np.isfinite(stats.rice.logpdf(x, b=0.)).all())
+ npt.assert_(np.isfinite(stats.rice.cdf(x, b=0.)).all())
+ npt.assert_(np.isfinite(stats.rice.logcdf(x, b=0.)).all())
+
+ q = [0.1, 0.1, 0.5, 0.9]
+ npt.assert_(np.isfinite(stats.rice.ppf(q, b=0.)).all())
+
+ mvsk = stats.rice.stats(0, moments='mvsk')
+ npt.assert_(np.isfinite(mvsk).all())
+
+ # furthermore, pdf is continuous as b\to 0
+ # rice.pdf(x, b\to 0) = x exp(-x^2/2) + O(b^2)
+ # see e.g. Abramovich & Stegun 9.6.7 & 9.6.10
+ b = 1e-8
+ np.allclose(stats.rice.pdf(x, 0), stats.rice.pdf(x, b),
@rgommers Owner

I think you meant to use assert_allclose here?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/stats/tests/test_continuous_extra.py
((12 lines not shown))
+ q = [0.1, 0.1, 0.5, 0.9]
+ npt.assert_(np.isfinite(stats.rice.ppf(q, b=0.)).all())
+
+ mvsk = stats.rice.stats(0, moments='mvsk')
+ npt.assert_(np.isfinite(mvsk).all())
+
+ # furthermore, pdf is continuous as b\to 0
+ # rice.pdf(x, b\to 0) = x exp(-x^2/2) + O(b^2)
+ # see e.g. Abramovich & Stegun 9.6.7 & 9.6.10
+ b = 1e-8
+ np.allclose(stats.rice.pdf(x, 0), stats.rice.pdf(x, b),
+ atol = b, rtol=0)
+
+
+def test_rice_rvs():
+ # a supplement to the KS test (which is run separately)
@rgommers Owner

This test just checks that _rvs returns the right output shape. What does it have to do with a KS test?

@rgommers Owner

Ah, you meant check_distribution_rvs. Maybe remove the comment or be more explicit

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
Owner

Looks good. One thing to add would be to remove 'rice' from the list of slow distributions in test_continuous_basic.py

@ev-br
Collaborator

made changed suggested by @rgommers, rebased on master, pushed the rebased version.

@rgommers
Owner

Slight git mishap here....

@ev-br
Collaborator

ouch.
I screwed it up.
any suggestion of how to fix this mess? Or is it easier to close this, apply the changes to the master, and send a fresh PR?

@rgommers
Owner

You can create a new branch from current master, cherry-pick the relevant commits and then force push to this PR ($ git push origin newbranch:rice_enh).

@ev-br
Collaborator

@rgommers done just that, thanks a whole lot!
Just for me to understand (and avoid in the future): what was the error?
What I did was

  • fetch upstream master and merge it to the local master
  • rebased the feature branch on local master
  • tried pushing to the origin, which failed
  • pulled the feature branch from the origin, resolved merge conflicts,
  • pushed the feature branch back to origin As there was no force-pushing involved, I thought it'd to the right thing. Apparently not...
@josef-pkt
Collaborator

after rebase you need to force push to your github branch, to overwrite it

@rgommers
Owner

Yep. Step 1, 2 are correct, 3 needed -f.

@rgommers rgommers merged commit f1b8578 into from
@rgommers
Owner

Merged, thanks!

@ev-br ev-br deleted the branch
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
9 scipy/stats/distributions.py
@@ -5573,6 +5573,15 @@ class rice_gen(rv_continuous):
%(example)s
"""
+ def _argcheck(self, b):
+ return b >= 0
+
+ def _rvs(self, b):
+ # http://en.wikipedia.org/wiki/Rice_distribution
+ sz = self._size if self._size else 1
+ t = b/np.sqrt(2) + mtrand.standard_normal(size=(2, sz))
+ return np.sqrt((t*t).sum(axis=0))
+
def _pdf(self, x, b):
return x * exp(-(x-b)*(x-b)/2.0) * special.i0e(x*b)
View
2  scipy/stats/tests/test_continuous_basic.py
@@ -153,7 +153,7 @@
distmiss = [[dist,args] for dist,args in distcont if dist in distmissing]
distslow = ['rdist', 'gausshyper', 'recipinvgauss', 'ksone', 'genexpon',
- 'vonmises', 'vonmises_line', 'rice', 'mielke', 'semicircular',
+ 'vonmises', 'vonmises_line', 'mielke', 'semicircular',
'cosine', 'invweibull', 'powerlognorm', 'johnsonsu', 'kstwobign']
# distslow are sorted by speed (very slow to slow)
View
28 scipy/stats/tests/test_continuous_extra.py
@@ -132,5 +132,33 @@ def test_rdist_cdf_gh1285():
values, decimal=5)
+def test_rice_zero_b():
+ # rice distribution should work with b=0, cf gh-2164
+ x = [0.2, 1., 5.]
+ npt.assert_(np.isfinite(stats.rice.pdf(x, b=0.)).all())
+ npt.assert_(np.isfinite(stats.rice.logpdf(x, b=0.)).all())
+ npt.assert_(np.isfinite(stats.rice.cdf(x, b=0.)).all())
+ npt.assert_(np.isfinite(stats.rice.logcdf(x, b=0.)).all())
+
+ q = [0.1, 0.1, 0.5, 0.9]
+ npt.assert_(np.isfinite(stats.rice.ppf(q, b=0.)).all())
+
+ mvsk = stats.rice.stats(0, moments='mvsk')
+ npt.assert_(np.isfinite(mvsk).all())
+
+ # furthermore, pdf is continuous as b\to 0
+ # rice.pdf(x, b\to 0) = x exp(-x^2/2) + O(b^2)
+ # see e.g. Abramovich & Stegun 9.6.7 & 9.6.10
+ b = 1e-8
+ npt.assert_allclose(stats.rice.pdf(x, 0), stats.rice.pdf(x, b),
+ atol = b, rtol=0)
+
+
+def test_rice_rvs():
+ rvs = stats.rice.rvs
+ npt.assert_equal(rvs(b=3.).size, 1)
+ npt.assert_equal(rvs(b=3., size=(3, 5)).shape, (3, 5))
+
+
if __name__ == "__main__":
npt.run_module_suite()
Something went wrong with that request. Please try again.