[SciPy-User] rvs and broadcasting
Thu Nov 24 13:50:43 CST 2011
> 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
> 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]])
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
More information about the SciPy-User