[SciPy-User] scipy.stats.norm strange doc string

[nicky van foreest]

Hi,

Thanks for your answers. Very helpful.

Nicky

On 13 March 2011 21:25, Warren Weckesser <warren.weckesser at enthought.com> wrote:
>
>
> On Sun, Mar 13, 2011 at 2:52 PM, nicky van foreest <vanforeest at gmail.com>
> wrote:
>>
>> Hi,
>>
>> The doc string of scipy.stats.norm tells me that the location and
>> scale parameters are array-like. However, when I try to pass arrays to
>> the loc and scale keywords I get an error. Specifically:
>>
>>
>> In [1]: from scipy.stats import norm
>>
>> In [2]: import numpy as np
>>
>> In [3]: mu = np.array([1,1])
>>
>> In [4]: simga = np.array([1,1])
>>
>> In [5]: x = [0,1,2,3]
>>
>> In [6]: norm.pdf(x, loc = mu, scale = simga)
>>
>>
>> results in :
>>
>> ValueError: shape mismatch: objects cannot be broadcast to a single shape
>>
>> I do understand how to resolve this problem, but for my specific
>> purpose I would have liked to pass mu and sigma as arrays, that is, I
>> would have liked to achieve
>>
>>        tau = np.zeros([g,m])
>>        for i in range(g):
>>            tau[i] = p[i]*norm.pdf(x, loc=mu[i], scale = sigma[i])
>>
>>
>> in one pass.
>
>
> All the arguments, including x, are broadcast, so you have ensure that their
> shapes are all compatible.   This can be accomplished with some judicious
> use of np.newaxis.  Here's a complete version of your snippet, with a "loop"
> version and a broadcasting version:
>
> -----
> import numpy as np
> from scipy.stats import norm
>
> mu = np.array([1.0, 1.25])
> sigma = np.array([4.0, 5.0])
> p = np.array([0.25, 0.75])
> x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>
> g = p.shape[0]
> m = x.shape[0]
>
> # Compute tau in a loop.
> tau = np.empty([g,m])
> for i in range(g):
>     tau[i] = p[i]*norm.pdf(x, loc=mu[i], scale = sigma[i])
>
> # Compute tau with broadcasting.
> tau2 = p[:,np.newaxis] * norm.pdf(x, loc=mu[:,np.newaxis],
> scale=sigma[:,np.newaxis])
>
> print "tau:"
> print tau
>
> print
> print "tau2:"
> print tau2
> -----
>
> Warren
>
>>
>> BTW:
>>
>> I am using this code to fit a set of normal distributions to a given
>> (quite) general distribution function by using the EM algorithm. Is
>> this already coded somewhere in scipy? If not, is somebody interested
>> in me making this available on the scipy cookbook?
>>
>> bye
>>
>> Nicky
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>