Side note: I would really like to eliminate the use of the _size attribute of the distribution classes. It is set only in the rvs method, and accessed only in the _rvs method, so it should simply be an argument to _rvs. Making it an attribute is pointless. However, the docstring for rv_continuous includes _rvs in the methods that subclasses can override, so we're stuck with the _size attribute for backwards compatibility.

Yep, sub-optimal design but no good way to get rid of self._size. It shouldn't get in the way of this fix though right?

Change is nicely local, called from only one place. So should be able to get it reviewed/merged without too large risks right before a release. Existing test coverage for non-scalar loc/scale/size is very thin at the moment though, needs expanding.

Yep, sub-optimal design but no good way to get rid of self._size. It shouldn't get in the way of this fix though right?

Right, the code in this PR leaves the use of ._size as it is. In an early draft, I got rid of it, but then realized _rvs is effectively a public part of the API, so I undid those changes (grumble, grumble).

Existing test coverage for non-scalar loc/scale/size is very thin at the moment though, needs expanding.

Working on it...

Updated with some tests. I'm still thinking about some more thorough tests. Suggestions are welcome!

This is still WIP. The pearson3 implementation needs some work. The Pearson Type III distribution can be expressed as a shifted gamma distribution, so I'm looking into rewriting pearson3 in terms of gamma (instead of trying to fix pearson3.rvs()).

I quickly drafted a test based on the stackoverflow thread, which is sensitive to the order of broadcasting and drawing variates: https://github.com/ev-br/scipy/tree/pr/5526
The relevant commit is ev-br@8c07855, the second commit tries to cover discrete distributions as well, but it likely doesn't do it.

In master, this fails all over. With this PR, I only see a couple of failures in truncnorm which does if self.a > 0 somewhere.
[and that is quite a Pandora's box].

This test is just a suggestion.

EDIT: Also array-valued shapes: ev-br@fda4518

Thanks! It occurred to me that another simple test would be to use np.vectorize to create a function that only calls the rvs method with scalars. Then the test would compare the results of the vectorize'd method with the actual method; if the random seed is set before each call, they should return the same array.

Then the test would compare the results of the vectorize'd method with the actual method; if the random seed is set before each call, they should return the same array.

Actually, this isn't true, because some of the rvs methods use two random numbers from the random number generator to generate one random variate of the distribution. For example, betaprime._rvs is

def _rvs(self, a, b):
    sz, rndm = self._size, self._random_state
    u1 = gamma.rvs(a, size=sz, random_state=rndm)
    u2 = gamma.rvs(b, size=sz, random_state=rndm)
    return (u1 / u2)

Calling that function several times with scalars will not produce the same variates as calling it once with broadcast arguments or with size != 1. That's not a bug in betaprime. It's just an implementation detail, but it means the idea of comparing np.vectorize'd results to actual broadcast results won't work in general. The good news is that all the continuous distributions in the dists list in test_distributions.py passed the test except for betaprime, dgamma, exponnorm, nct, dweibull, and rice (all of which use more than one random number per variate) and for pearson3 (which I am still skipping).

With this PR, I only see a couple of failures in truncnorm which does if self.a > 0 somewhere.

truncnorm uses if self.a > 0 twice. That can be fixed by using np.where or _lazywhere. I'll include a fix for truncnorm when I update the PR.

It would also be helpful to have a test explicitly comparing the output of rvs and a corresponding numpy generator for several distributions, especially with those which take non-empty shape parameters. Eg. beta or binom.

Regarding the size handling, I guess we could (i) add size parameter to the signature of _rvs method of distributions we have in scipy.stats, (ii) keep the old code path and deprecate it in the docs and release notes with the intent of removing it in a couple of releases, (iii) in code, inspect the signature of self._rvs and use either the old or the new code paths (this can be similar to how _stats_has_moments is handled).

This can be done in a separate PR, of course :-). I just think it would be nice to change both things in one release if possible.

Even without the python 3 test failures, this is still WIP. I haven't touched pearson3 yet.

@ev-br I didn't use your commit, but putting the test "kernel" in common_tests.py was definitely the way to go. I think the tests now cover what you had in your commits, and more.

@rgommers: In case it isn't obvious: this PR won't be ready in time for 0.17. Making the rvs() arguments work correctly might end up requiring that all the distribution methods properly broadcast. In theory, that's supposed to be the case already (see, for example, http://docs.scipy.org/doc/scipy/reference/tutorial/stats.html#broadcasting). In practice, not so much.

So we can remove the 0.17.0 label from gh-2069?

So we can remove the 0.17.0 label from gh-2069?

Yup, I just did.

Any time scale for completing this? Just hit the bug too.

@pwaller I've been on a bit of a hiatus from scipy development, but I plan to get this into 0.18.

In the meantime, what is the best workaround which won't break once this has landed? Is there any? I'm trying to generate many rvs with many values of parameters, e.g: beta.rvs([1, 2], [1, 2], size=(10, 2)) (I want 10 rvs each for 2 parameters). But I can't find any series of parameters which make that work. I guess I have to live with a non-vectorized implementation for now?

Edit: in the end I settled with array(list(beta(a, b).rvs(size=10000) for a, b in zip(As, Bs)))which is moderately fast, taking ~500ms for ~100 sets of parameters.

Despite my hiatus in scipy contributions, I haven't forgotten about this.

I think the tentative plan is to release 0.18 this month. Is there a more concrete timeline for 0.18?

The milestone page lists June 14, but judging based on the number of PRs marked for 0.18, I'm skeptical :) I don't know of anything more concrete than that, at least.

Thanks Eric. Time to dive back into this...

Updated.

This would need a mention in the release notes, in big bold letters.

@ev-br I added a note to the release notes about the bug fix. It is in the "Backwards incompatible changes" section, mainly because there are cases with arguments that are incompatible with a given size where the function used to return values and now it raises an error. My guess is that case will be pretty rare in the wild (but I admit it is just a guess).

I refrained from using "big bold letters". 😄

That's big enough and bold enough for me :-). This PR looks good to me now. I'm going to merge it before branching unless there are further comments.

Needs a rebase (again).

Here's something I don't understand:

In [22]: scipy.__version__
Out[22]: '0.18.0.dev0+36683cf'

In [23]:  stats.norm.rvs(size=4, random_state=np.random.RandomState(1234))
Out[23]: array([ 0.47143516, -1.19097569,  1.43270697, -0.3126519 ])

In [24]: stats.levy_stable.rvs(alpha=1, beta=0, size=4, random_state=np.random.RandomState(1234))
Out[24]: array([-0.12827986, -0.73108696, -2.39290307, -0.17495177])

In [25]: stats.levy_stable.rvs(alpha=1.5, beta=0.3, size=4, random_state=np.random.RandomState(1234))
Out[25]: array([-0.40182381,  1.65060746, -1.64155074,  0.49261358])

In scipy 0.17:

In [16]: scipy.__version__
Out[16]: '0.17.0'

In [17]: stats.norm.rvs(size=4, random_state=np.random.RandomState(1234))
Out[17]: array([ 0.47143516, -1.19097569,  1.43270697, -0.3126519 ])

In [18]:  stats.levy_stable.rvs(alpha=1, beta=0, size=4, random_state=np.random.RandomState(1234))
Out[18]: array([-1.13381723,  0.56711265,  3.39976097,  0.03294947])

In [19]:  stats.levy_stable.rvs(alpha=1.5, beta=0.3, size=4, random_state=np.random.RandomState(1234))
Out[19]: array([ 0.67579527, -0.76668594,  0.82701716,  1.36710987])

Notice that the random draws are the same for norm (I also spot-checked gamma and rice), as they should. For levy_stable however, they differ.

Looks like levy_stable._rvs() wasn't updated along with the rest of the distributions when the _random_state attribute was added. In scipy 0.17.1:

In [4]: stats.levy_stable.rvs(alpha=1, beta=0, size=4, random_state=np.random.RandomState(1234))
Out[4]: array([ 1.84349104, -2.86730347, -3.2989345 ,  0.44281769])

In [5]: stats.levy_stable.rvs(alpha=1, beta=0, size=4, random_state=np.random.RandomState(1234))
Out[5]: array([-0.07054039,  1.90011002, -3.46018551, -4.57797445])

Those outputs should be the same.

I'll fix that now.

Ouch. levy_stable is skipped in the loop over all distributions which, among other things, calls check_random_state_property, so it was not noticed.
This actually was my mess, thank you for cleaning it up!

Updated with fix for levy_stable.

LGTM, merging. Thanks Warren!