[SciPy-User] rvs and broadcasting

Friedrich Romstedt friedrichromstedt@gmail....
Thu Nov 24 13:50:43 CST 2011

2011/11/24  <josef.pktd@gmail.com>:
> 1) broadcast shape parameters, loc and scale, if they are arrays
> produce rvs in that shape, and, if in this case size is not the same
> or 1, then raise a ValueError
> essentially
>    lower, upper, loc, scale = np.broadcast_arrays(lower, upper, loc, scale)
>    if (np.size(lower) > 1) and (size != (1,)) and (lower.shape != size):
>        raise ValueError('Do you really want this? Then do it yourself.')
> 2) broadcast shape parameters, loc and scale,  for each of these
> create random variables given by size, the return shape is essentially
> broadcasted shape concatenated with size, for example
> assert_equal(truncnorm_rvs(lower*np.arange(4)[:,None], upper,
>                                loc=np.arange(5), scale=1, size=(2,3)).shape,
>                     (4, 5, 2, 3))
> this version is attached.
> Any opinions about which version should be preferred?
> (As aside, truncnorm and other distribution with parameter dependent
> support might also have other problems than just the broadcasting of
> shape and scale.)

I don't know the context but the first solution is *a lot* more
readable.  I'm not even interested in understanding the cryptic second

I don't know if all of the checks done in the ``if`` scope are really
necessary, but I just believe you.

Alternatively, since you essentially want to check if things are
broadcastable AISI, why not just broadcasting them with a result array
of shape ``size`` you need to generate anyway along the way,
broadcasting all three and catching the ``ValueError: shape mismatch``

Broadcasting is cheap since it plays stride stricks IIRC:

>>> x = numpy.asarray([1, 2])
>>> x2, y = numpy.broadcast_arrays(x, [[1, 2], [3, 4]])
>>> x2.strides
(0, 4)

Would maybe make the algorithm more general?

Might be there are some side effects I'm overlooking atm.  I'm not
paying too much attention, since I already made my point that I would
prefer the first one, and my suggestion is not crucial since you'll
test it anyway.


P.S.: Maybe there is a numpy function around to just test
compatibility instead of testing + broadcasting at the same time in
direct succession?

More information about the SciPy-User mailing list