[SciPy-user] Creating a Matrix from a Sum

David Warde-Farley dwf@cs.toronto....
Fri Jul 3 14:01:44 CDT 2009

On 3-Jul-09, at 2:45 PM, David Warde-Farley wrote:

> On 2-Jul-09, at 9:25 AM, Lorenzo Isella wrote:
>> Dear All,
>> I need some help to efficiently write an array manipulation. I  
>> suspect
>> this could be a one-liner.
>> Assume that you have a histogram of observations that you store in a
>> vector x.
>> Let us say that its i-th entry, x_i, corresponds to the number of
>> observations in the i-th channel, for i=1,2...N. (or 0,1...N-1,  
>> please
>> let me know if there is any potential 0/1 pitfall in the following).
>> Now, for any two channels i and j, I want to calculate the  
>> probability
>> of having an observation  in any other channel k, where  k>=max(i,j).
>> That is to say
>> P(i,j)=sum_{k=max(i,j)}^N x_k /C,
>> where C=sum_{i=1}^N x_i is just a normalization factor.
> Hm, there might be a simpler/more efficient way but this is what I
> came up with as far as 1-liners:
> P = np.max(np.concatenate(np.broadcast_arrays(x[np.newaxis, :,
> np.newaxis], x[:, np.newaxis, np.newaxis]),axis=2),axis=2); P /=  
> x.sum()

Oops, nevermind, that's wrong.

Here's how I would do it.

You want the sum from k to the end (N) so I'd extract a vector like  
this (the cumulative sum of the reversed array, reversed again

	sums = np.cumsum(x[::-1])[::-1]

Then you want the max of i and j at each position in the matrix. For  
this, though there may be a more efficient way, you can use mgrid to  
get i and j indices at each position in the matrix and then use  
concatenate and max on them to concatenate along a third axis and then  
do a max along that same axis to recover a 2D array:

	idx = np.concatenate([M[:,:,np.newaxis] for M in mgrid[0:len(x),  

Finally, use fancy indexing and divide by your normalizer:

	P = sums[idx]  / x.sum()

Or, as a oneliner, much less readable:

	P = cumsum(x[::-1])[::-1][np.concatenate([M[:,:,np.newaxis] for M in  
mgrid[0:len(x), 0:len(x)]],axis=2).max(axis=2)] / x.sum()



More information about the SciPy-user mailing list