Scipy 曲线拟合全局参数

Question

我正在尝试在 python 中使用 Scipy curve_fit 来拟合银河模型的两个全局参数。我有一个自变量数组和一个因变量数组。数据集的前 1/4 需要根据两个全局参数和两个局部参数拟合到一个函数，下一个季度根据两个全局参数和两个局部变量拟合到另一个函数等。

我是否可以编写一个函数，通过整个数组调用具有正确索引和全局参数的适当函数。

到目前为止我所拥有的是：

   def galaxy_func_inner(time,a,b,c,d):
     telescope_inner = lt.station(rot_angle=c,pol_angle=d)
     power = telescope_inner.calculate_gpowervslstarray(time)[0]
     return a*np.array(power)+b

   def galaxy_func_outer(time,a,b,c,d):
     telescope_outer = lt.station(rot_angle=c,pol_angle=d)
     power = telescope_outer.calculate_gpowervslstarray(time)[0]
     return a*np.array(power)+b

   def galaxy_func_global(time,R,P,a,b,c,d,e,f,g,h):
        for t_index in range(len(time)):
           if t_index in range(0,50):
                   return galaxy_func_outer(t_index,a,b,R,P)
           elif t_index in range(50,100):
                   return galaxy_func_outer(t_index,c,d,R,P)
           elif t_index in range(100,150):
                   return galaxy_func_inner(t_index,e,f,R,P)
           elif t_index in range(150,200):
                   return galaxy_func_inner(t_index,g,h,R,P)

问题是这个只适合第一次但是整个时间数组，并且单点只适合对应的模型点而不是整个数组。关于如何重新制定这个的任何帮助？我试图将其重新表述为：

 def galaxy_func_global(xdata,R,P,a,b,c,d,e,f,g,h):
     return galaxy_func_outer(xdata[0:50],a,b,R,P),galaxy_func_outer(xdata[50:100],c,d,R,P),galaxy_func_inner(xdata[100:150],e,f,R,P),galaxy_func_inner(xdata[150:200],g,h,R,P)

但我得到了错误：

File "galaxy_calibration.py", line 117, in <module>
    popt,pcov = curve_fit(galaxy_func_global,xdata,ydata)
File "/Library/Python/2.7/site-packages/scipy-0.14.0.dev_7cefb25-py2.7-macosx-10.9-intel.egg/scipy/optimize/minpack.py", line 555, in curve_fit
    res = leastsq(func, p0, args=args, full_output=1, **kw)
File "/Library/Python/2.7/site-packages/scipy-0.14.0.dev_7cefb25-py2.7-macosx-10.9-intel.egg/scipy/optimize/minpack.py", line 369, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
File "/Library/Python/2.7/site-packages/scipy-0.14.0.dev_7cefb25-py2.7-macosx-10.9-intel.egg/scipy/optimize/minpack.py", line 20, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
File "/Library/Python/2.7/site-packages/scipy-0.14.0.dev_7cefb25-py2.7-macosx-10.9-intel.egg/scipy/optimize/minpack.py", line 445, in _general_function
    return function(xdata, *params) - ydata
ValueError: operands could not be broadcast together with shapes (4,) (191,) 

任何帮助将不胜感激。

Answer 1

如果您想将输入数据分成 4 个批次（基于 time 点的索引）并根据批次处理数据，然后将结果返回到单个数组中，那么你可以这样做：

def galaxy_func_global(time,R,P,a,b,c,d,e,f,g,h):
    return np.concatenate([galaxy_func_outer(time[0:50],a,b,R,P),
                          galaxy_func_outer(time[50:100],c,d,R,P),
                          galaxy_func_inner(time[100:150],e,f,R,P),
                          galaxy_func_inner(time[150:200],g,h,R,P)])

这将切入您的time 数组以挑选出每个感兴趣的切片，然后为每个切片调用适当的函数。在我看来，这些函数返回简单的np.arrays，可以连接起来得到一个数组作为结果。

（我刚刚意识到我可以说“您尝试的几乎是完美的，但是您需要将结果数组连接成一个数组”:)

---

请注意，至少有两种方式可能会导致尺寸标注问题。

首先，您应该确保两个函数 (galaxy...inner/outer()) 的返回值都是一维 numpy 数组。否则你会遇到全局返回值的问题。

其次，出于显而易见的原因，每个拟合方法都需要一个函数，该函数的返回值与输入变量具有相同的大小（形状）。因此，如果 time 的长度不完全是 200 个元素，您当前的代码也会遇到问题，因为即使 time 更长，您的输出也会被截断为 200 个元素。至少你应该放

galaxy_func_inner(time[150:],g,h,R,P)

进入您的最后一个函数调用以捕获time 的所有剩余点，但如果您想正确执行此操作，请调用

def galaxy_func_global(time,R,P,a,b,c,d,e,f,g,h):
    inds=np.floor(np.linspace(0,len(time)-1,5))
    return np.concatenate([galaxy_func_outer(time[0:inds[1]],a,b,R,P),
                      galaxy_func_outer(time[inds[1]:inds[2]],c,d,R,P),
                      galaxy_func_inner(time[inds[2]:inds[3]],e,f,R,P),
                      galaxy_func_inner(time[inds[3]:],g,h,R,P)])

另请注意，您的原始错误正式属于这种类型：

File "/Library/Python/2.7/site-packages/scipy-0.14.0.dev_7cefb25-py2.7-macosx-10.9-intel.egg/scipy/optimize/minpack.py", line 445, in _general_function
    return function(xdata, *params) - ydata
ValueError: operands could not be broadcast together with shapes (4,) (191,) 

这告诉你python不能从function(xdata,*params)（即你的拟合模型）中减去ydata，因为一个长度为4而另一个长度为191。这是因为如果你的函数调用return a,b,c,d，然后它将返回一个 tuple (a,b,c,d)，因此返回值的长度为 4。更有趣的是，您的 ydata 的长度为 191，这可能意味着您仍然会遇到一个错误。