NumPy中的矢量化运算

Question

我有一个由5个numpy数组组成的numpy数组，每个数组有5个观测值。我需要使用积分的梯形规则来获得5个积分结果。我正在尝试使用for循环将函数应用于每个numpy数组，并寻找实现速度更快的方法。

def trapezoid_rule(a,b,n,f_nparray):
h = (b-a) / (n-1)
return (h/2) * (f_nparray[0] + 2 * sum(f_nparray[1:n-1]) + f_nparray[n-1])
simulation_res = np.array([[1,2,3,4,5], [1,3,5,7,9], [1,4,2,5,7], [1,5,2,6,4], [6,2,5,3,4]])
# List of integration results
I = []
for i in simulation_res:
I.append(trapezoid_rule(0,1,10,i))

预期输出格式=[a，b，c，d，e]

Answer 1

您可以在没有for循环的情况下执行此操作，但我将同时显示这两个循环。numpy有一个以矢量化方式工作的内置梯形规则。如果你想使用它：

import numpy as np
simulation_res = np.array([[1,2,3,4,5], 
[1,3,5,7,9], 
[1,4,2,5,7], 
[1,5,2,6,4], 
[6,2,5,3,4]])
result = np.trapz(simulation_res, axis=1)
print(result) 
# Gives [12.  20.  15.  15.5 15. ]

但是，如果由于某种原因需要循环，您可以在numpy数组中预先分配一些内存，并将其附加到

result = np.zeros(5)
for idx, i in enumerate(simulation_res):
result[idx] = np.trapz(i)
print(result)
# Gives [12.  20.  15.  15.5 15. ]

Answer 2

带有：

def foo(a,b,n,f_nparray):
h = (b-a) / (n-1)
return (h/2) * (f_nparray[:,0] + 2 * np.sum(f_nparray[:,1:n-1], axis=1) + f_nparray[:,n-1])
In [20]: foo(0,1,5,simulation_res)
Out[20]: array([3.   , 5.   , 3.75 , 3.875, 3.75 ])

与您的CCD_ 4相匹配。

我不得不将axis添加到np.sum；与我的评论相反，只使用0进行索引。我只需要2d数组的每行一个值。

与np.trapz:相比

In [29]: np.trapz(simulation_res, axis=1)/4
Out[29]: array([3.   , 5.   , 3.75 , 3.875, 3.75 ])
In [30]: timeit np.trapz(simulation_res, axis=1)/4
34 µs ± 1.66 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [31]: timeit foo(0,1,5,simulation_res)
25.7 µs ± 139 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)