making a function that can take arguments in various shapes
<p>Q1)
Numpy functions can take arguments in different shapes. For instance, np.sum(V) can take either of two below and return outputs with different shapes.</p>

<pre><code>x1= np.array( [1,3] ) #(1)
x2= np.array([[[1,2],[3,4]], [[5,6],[7,8]]]) #(2)
</code></pre>

<p>I am making my own function something like below, which adds two values in an 1D vector with the length of two and return the real number.</p>

<pre><code>def foo(V):
    return V[0]+V[1];
</code></pre>

<p>However, this foo function can only take one 1D vector and cannot take any other shapes. It can only take x1 above as an argument but not x2. If I want to make my function work with either of two variables above(x1 and x2), or with any other shapes that has arrays with the length of 2 in their last dimension, how should I revise my foo function?</p>

<hr>

<p>---------------------------update------------------------------</p>

<p>My original function was a hardcoded negative gaussian pdf function.</p>

<pre><code>def nGauss(X, mu, cov):
    # multivariate negative gaussian.    
    # mu is a vector and cov is a covariance matrix.

    k = X.shape[0];
    dev = X-mu
    p1 = np.power( np.power(np.pi * 2, k) , -0.5);
    p2 = np.power( np.linalg.det(cov)  , -0.5)
    p3 = np.exp( -0.5 * np.dot( np.dot(dev.transpose(), np.linalg.inv(cov)), dev));

    return -1.0 * p1 * p2 * p3;
</code></pre>

<p>Now his function can return only one pdf value. For example, it can only take arguments like np.array([1,2]), but cannot take arguments X like np.array([[[1,2], [5,6]], [[7,8],[9,0]]]). Here my question was how to make my gaussian function takes arguments of arbitrary shapes and return the pdf value of each point maintaining the same structure except the last dimension, such as
<code>nGauss(np.array( [1,2] ), mu, cov)</code> returns [ 0.000023 ], and 
<code>nGauss(np.array([[[1,2], [5,6]], [[7,8],[9,0]]]), mu, cov)</code> returns [[ 0.000023, 0000014], [0.000012, 0.000042]].</p>

<p>I notice that scipy function 'multivariate_normal.pdf' can do this.</p>

<hr>

<p>Q2)
I am also having a difficulty in understanding np's basic array.</p>

<pre><code>t1=np.array([[1,2,3], [4,5,6]])
t2=np.array([1,2,3])
t3=np.array([[1,2,3], [4,5,6],5])
</code></pre>

<p>The shape of t1 is (2,3), and it seems legitimate in terms of matrix perspective; 2 rows and 3 columns. However, the shape of t2 is (3,), which I think has to be (3). What's the meaning of the empty space after "3,"? Also, the shape of t3 is (3,). In this case, is the meaning of the empty space that dimensions vary? </p>

<p>In advance, thanks for your help.</p>

<p>Your function works with both arrays:</p>

<pre><code>In [1]: def foo(V):
   ...:     return V[0]+V[1]
   ...: 
In [2]: foo(np.array([1,3]))
Out[2]: 4
In [3]: foo(np.array([[[1,2],[3,4]], [[5,6],[7,8]]]))
Out[3]: 
array([[ 6,  8],
       [10, 12]])
</code></pre>

<p>This answer is just the sum of these two arrays:</p>

<pre><code>In [4]: np.array([[[1,2],[3,4]], [[5,6],[7,8]]])[0]
Out[4]: 
array([[1, 2],
       [3, 4]])
In [5]: np.array([[[1,2],[3,4]], [[5,6],[7,8]]])[1]
Out[5]: 
array([[5, 6],
       [7, 8]])
</code></pre>

<p>If you expected something else, you'll have to show us.</p>

<p>As for your second question:</p>

<pre><code>In [6]: t1=np.array([[1,2,3], [4,5,6]])
   ...: t2=np.array([1,2,3])
   ...: t3=np.array([[1,2,3], [4,5,6],5])
   ...: 
In [7]: t1.shape
Out[7]: (2, 3)
In [8]: t2.shape
Out[8]: (3,)
In [9]: t3.shape
Out[9]: (3,)
</code></pre>

<p><code>(3,)</code> is a 1 element tuple. Compare these expressions.</p>

<pre><code>In [11]: (3)
Out[11]: 3
In [12]: (3,)
Out[12]: (3,)
</code></pre>

<p>There have been several recent questions about (3,) v (3,1) shape arrays, and <code>np.array([[1,2,3]])</code> v. <code>np.array([1,2,3])</code>.</p>

<p><code>t3</code> is an object dtype array, with 3 elements.  The 3 inputs are different length, so it can't create a 2d array.  Stay away from this type of array for now.  Focus on the simpler arrays.</p>

<pre><code>In [10]: t3
Out[10]: array([[1, 2, 3], [4, 5, 6], 5], dtype=object)
In [13]: t3[0]
Out[13]: [1, 2, 3]
In [14]: t3[2]
Out[14]: 5
</code></pre>

<p><a href="http://stackoverflow.com/questions/39706277/numpy-why-is-difference-of-a-2-1-array-and-a-vertical-matrix-slice-not-a-2-1">Numpy: Why is difference of a (2,1) array and a vertical matrix slice not a (2,1) array</a></p>

<p><a href="http://stackoverflow.com/questions/39694318/difference-between-single-and-double-bracket-numpy-array">Difference between single and double bracket Numpy array?</a></p>

<p>=====================</p>

<p>With the <code>nGauss</code>:</p>

<pre><code>In [53]: mu=np.array([0,0])
In [54]: cov=np.eye(2)
In [55]: xx=np.array([[[1,2], [5,6]], [[7,8],[9,0]]])
In [56]: np.apply_along_axis(nGauss, -1, xx, mu, cov)
Out[56]: 
array([[ -1.30642333e-02,  -9.03313360e-15],
       [ -4.61510838e-26,  -4.10103631e-19]])
</code></pre>

<p><code>apply_along_axis</code> iterates on the 1st 2 dim, passing each <code>xx[i,j,:]</code> to <code>nGauss</code>.  It's not fast, but is relatively easy to apply.</p>

<pre><code>k = X.shape[0];  # I assume you want
k = X.shape[[1]   # the last dimension
dev = X-mu     # works as long as mu has k terms
</code></pre>

<p>this is a scalar:</p>

<pre><code>p1 = np.power( np.power(np.pi * 2, k) , -0.5);
</code></pre>

<p>so is</p>

<pre><code>p2 = np.power( np.linalg.det(cov)  , -0.5)
</code></pre>

<p>So it comes down to generalizing this expression:</p>

<pre><code>p3 = np.exp( -0.5 * np.dot( np.dot(dev.transpose(), np.linalg.inv(cov)), dev));
</code></pre>

<p>In the simple (2,) <code>x</code> case, <code>dev</code> is 1d, and <code>dev.transpose()</code> does nothing.</p>

<p>It's easier to generalize <code>einsum</code> than <code>dot</code>; I think the equivalent is:</p>

<pre><code>p3 = np.einsum('j,j', np.einsum('i,ij', dev, np.linalg.inv(cov)), dev)
p3 = np.exp( -0.5 * p3)
</code></pre>

<p>which simplifies to</p>

<pre><code>p3 = np.einsum('i,ij,j', dev, np.linalg.inv(cov), dev)
</code></pre>

<p>generalizing to higher dim:</p>

<pre><code>p3 = np.einsum('...i,ij,...j', dev, np.linalg.inv(cov), dev)
</code></pre>

<p>So with:</p>

<pre><code>def nGaussA(X, mu, cov):
    # multivariate negative gaussian.    
    # mu is a vector and cov is a covariance matrix.

    k = X.shape[-1];
    dev = X-mu
    p1 = np.power( np.power(np.pi * 2, k) , -0.5);
    p2 = np.power( np.linalg.det(cov)  , -0.5)
    p3 = np.einsum('...i,ij,...j', dev, np.linalg.inv(cov), dev)
    p3 = np.exp( -0.5 * p3)
    return -1.0 * p1 * p2 * p3;
</code></pre>

<p>matching earlier values:</p>

<pre><code>In [85]: nGaussA(x,mu,cov)
Out[85]: -0.013064233284684921
In [86]: nGaussA(xx,mu,cov)
Out[86]: 
array([[ -1.30642333e-02,  -9.03313360e-15],
       [ -4.61510838e-26,  -4.10103631e-19]])
</code></pre>

<p>So the way to generalize the function is to check each step.  If it produces a scalar, keep it.  If operates with an <code>x</code> keep it.  But if it requires coordinating dimensions with other arrays, use a numpy operation that does that.  Often that involves broadcasting.  Sometimes it helps to study other numpy functions to see how they generalize (e.g. <code>apply_along_axis</code>, <code>apply_over_axes</code>, <code>cross</code>, etc).</p>

<p>An interactive numpy session is essential; allowing me to try ideas with small sample arrays.</p>

python numpy