[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
one.
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``
exception?
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.
Friedrich
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