使用numpy或scipy将3D数据数组拟合到1D函数上

Question

使用numpy或scipy将3D数据数组拟合到1D函数上

4

我目前正在尝试将大量数据拟合到正弦函数。在只有一个数据集（1D数组）的情况下，scipy.optimize.curve_fit（）可以很好地工作。但是，据我所见，如果函数本身仅为一维，则不允许更高维度的数据输入。我不想使用for循环迭代数组，因为在Python中这种方法速度非常慢。

到目前为止，我的代码应该类似于这样：

from scipy import optimize
import numpy as np    
def f(x,p1,p2,p3,p4): return p1 + p2*np.sin(2*np.pi*p3*x + p4)      #fit function

def fit(data,guess):
   n = data.shape[0] 
   leng = np.arange(n)
   param, pcov = optimize.curve_fit(f,leng,data,guess)
   return param, pcov

数据是一个三维数组(shape=(x,y,z))，我想将每一行data[:,a,b]拟合到函数中，其中param是一个形状为(4,y,z)的数组。当然，对于多维数据，这会导致以下错误：

ValueError: operands could not be broadcast together with shapes (2100,2100) (5)

也许有一个简单的解决方案，但我不确定该怎么做。有什么建议吗？

搜索答案证明很困难，因为大多数与这些关键字相关的主题都涉及更高维度的函数拟合。

- Linus-

不用担心for循环太小的问题。我相信曲线拟合才是你代码中的瓶颈。如果你怀疑循环是瓶颈，那就对代码进行分析！ - David Zwicker

好的，这个想法是如果曲线拟合可以使用整个数组而不是运行函数y*z次，那么它将会快得多。这就是我所说的for循环速度慢的意思。 - Linus-

也许你可以通过FFT获取sin的参数。 - HYRY

为什么曲线拟合应该更快？大部分时间可能会花在拟合程序上，而不是循环遍历数据上。 - David Zwicker

我原本以为它会更快，因为在比较使用for循环遍历数据数组和直接将整个数组输入函数时，我在其他程序中的经验是如此（高达100倍的因素）。我承认那些计算相对简单，也许当使用像curve_fit这样更耗时的方法时节省的时间可能不那么关键，但我仍然希望进程能够稍微加速一点。我也尝试过使用FFT，但由于采样点数量很少，结果并不是很令人满意。 - Linus-

尽管如此，在找不到其他解决方案的情况下，我会再次尝试使用FFT。 - Linus-

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Saullo G. P. Castro · Accepted Answer

使用 np.apply_along_axis() 可以解决您的问题。只需按照以下步骤操作：

func1d = lambda y, *args: optimize.curve_fit(f, xdata=x, ydata=y, *args)[0] #<-- [0] to get only popt
param = np.apply_along_axis( func1d, axis=2, arr=data )

请看下面的示例：

from scipy import optimize
import numpy as np
def f(x,p1,p2,p3,p4):
    return p1 + p2*np.sin(2*np.pi*p3*x + p4)
sx = 50  # size x
sy = 200 # size y
sz = 100 # size z
# creating the reference parameters
tmp = np.empty((4,sy,sz))
tmp[0,:,:] = (1.2-0.8) * np.random.random_sample((sy,sz)) + 0.8
tmp[1,:,:] = (1.2-0.8) * np.random.random_sample((sy,sz)) + 0.8
tmp[2,:,:] = np.ones((sy,sz))
tmp[3,:,:] = np.ones((sy,sz))*np.pi/4
param_ref = np.empty((4,sy,sz,sx))     # param_ref in this shape will allow an
for i in range(sx):                    # one-shot evaluation of f() to create 
    param_ref[:,:,:,i] = tmp           # the data sample
# creating the data sample
x = np.linspace(0,2*np.pi)
factor = (1.1-0.9)*np.random.random_sample((sy,sz,sx))+0.9
data = f(x, *param_ref) * factor       # the one-shot evalution is here
# finding the adjusted parameters
func1d = lambda y, *args: optimize.curve_fit(f, xdata=x, ydata=y, *args)[0] #<-- [0] to get only popt
param = np.apply_along_axis( func1d, axis=2, arr=data )