如何在numpy中創建一個n維網格來評估任意n的函數？

Question

我試圖創建一個天真的數值積分函數來說明蒙特卡洛積分在高維上的好處。我想是這樣的：

def quad_int(f, mins, maxs, numPoints=100): 
    ''' 
    Use the naive (Riemann sum) method to numerically integrate f on a box 
    defined by the mins and maxs. 

    INPUTS: 
    f   - A function handle. Should accept a 1-D NumPy array 
     as input. 
    mins  - A 1-D NumPy array of the minimum bounds on integration. 
    maxs  - A 1-D NumPy array of the maximum bounds on integration. 
    numPoints - An integer specifying the number of points to sample in 
     the Riemann-sum method. Defaults to 100. 
    ''' 
    n = len(mins) 

    # Create a grid of evenly spaced points to evaluate f on 
    # Evaluate f at each point in the grid; sum all these values up 
    dV = np.prod((maxs-mins/numPoints)) 
    # Multiply the sum by dV to get the approximate integral 

我知道我的dV將會導致與數值穩定性問題，但現在我在與正在創建的域什麼麻煩。如果維度的數量是固定的，這將是很容易只使用np.meshgrid這樣的：

# We don't want the last value since we are using left-hand Riemann sums 
x = np.linspace(mins[0],maxs[0],numPoints)[:-1] 
y = np.linspace(mins[1],maxs[1],numPoints)[:-1] 
z = np.linspace(mins[2],maxs[2],numPoints)[:-1] 

X, Y, Z = np.meshgrid(x,y,z) 
tot = 0 
for index, x in np.ndenumerate(X): 
    tot += f(x, Y[index], Z[index]) 

是否有一個模擬到np.meshgrid可以爲任意尺寸做到這一點，也許接受陣列的一個元組？或者是否有其他方式可以在更高維度上進行黎曼和？我想過遞歸地做這件事，但無法弄清楚它是如何工作的。

Answer 1

您可以使用列表理解來生成所有的linspaces，然後將該列表傳遞給meshgrid並使用*（將列表轉換爲參數元組）。

XXX = np.meshgrid(*[np.linspace(i,j,numPoints)[:-1] for i,j in zip(mins,maxs)]) 

XXX現在是n陣列（每個n維）的列表。

我正在使用直接的Python列表和參數操作。

np.lib.index_tricks有其他的索引和網格生成函數和可能有用的類。值得一讀，看看如何做。

---

索引未知維的數組時，各種numpy函數中使用的巧妙方法是構造爲所需索引的列表。它可以包括slice(None)，你通常會看到:。然後將其轉換爲元組並使用它。

In [606]: index=[2,3] 
In [607]: [slice(None)]+index 
Out[607]: [slice(None, None, None), 2, 3] 
In [609]: Y[tuple([slice(None)]+index)] 
Out[609]: array([ 0. , 0.5, 1. , 1.5]) 
In [611]: Y[:,2,3] 
Out[611]: array([ 0. , 0.5, 1. , 1.5]) 

他們使用一個列表，他們需要更改元素。轉換爲元組並不總是需要

index=[slice(None)]*3 
index[1]=0 
Y[index] # same as Y[:,0,:]