首页 文章

在Scipy中使用curve_fit拟合矢量函数

提问于
浏览 次
2

我想使用Scipy的curve_fit（或更合适的东西，如果可用的话）使用矢量输出拟合函数 . 例如，请考虑以下功能：

import numpy as np
def fmodel(x, a, b):
    return np.vstack([a*np.sin(b*x), a*x**2 - b*x, a*np.exp(b/x)])

每个组件都是不同的功能，但它们共享我希望适合的参数 . 理想情况下，我会做这样的事情：

x = np.linspace(1, 20, 50)
a = 0.1
b = 0.5
y = fmodel(x, a, b)
y_noisy = y + 0.2 * np.random.normal(size=y.shape)

from scipy.optimize import curve_fit
popt, pcov = curve_fit(f=fmodel, xdata=x, ydata=y_noisy, p0=[0.3, 0.1])

但 curve_fit 不能与带矢量输出的函数一起使用，并抛出错误 Result from function call is not a proper array of floats. . 我所做的是将输出变平如下：

def fmodel_flat(x, a, b):
    return fmodel(x[0:len(x)/3], a, b).flatten()

popt, pcov = curve_fit(f=fmodel_flat, xdata=np.tile(x, 3),
                       ydata=y_noisy.flatten(), p0=[0.3, 0.1])

这很有效 . 如果不是向量函数，我实际上也使用不同的输入来拟合几个函数但是共享模型参数，我可以连接输入和输出 .

有没有更合适的方法来适应Scipy的矢量函数或者一些额外的模块？对我来说主要考虑的是效率 - 适合的实际功能要复杂得多，适合可能需要一些时间，所以如果 curve_fit 的使用被破坏并导致运行时间过长，我想知道我应该做些什么 .

python numpy scipy curve-fitting

2 回答

  • 1

    从效率的角度来看，我认为你所做的事情是完美无缺的 . 我会试着看一下实现并提出一些更加量化的东西，但暂时还是我的推理 .

    您在曲线拟合期间所做的是优化参数 (a,b)

    res = sum_i |f(x_i; a,b)-y_i|^2

    是最小的 . 我的意思是你有任意维度的数据点 (x_i,y_i) ，两个参数 (a,b) 和近似查询点 x_i 的数据的拟合模型 .

    曲线拟合算法从一个起始 (a,b) 对开始，将其放入计算上述平方误差的黑盒子中，并尝试提出一个产生较小误差的新的 (a',b') 对 . 我的观点是上面的错误实际上是拟合算法的黑盒子：拟合的配置空间仅由 (a,b) 参数定义 . 如果你想象你如何实现一个简单的曲线拟合函数，你可以想象你试图做一个梯度下降，方形误差作为成本函数 .

    现在，黑盒子如何计算错误应该与拟合程序无关 . 很容易看出 x_i 的维度与标量函数无关，因为如果你有1000个1d查询点适合，或者在3d空间中有10x10x10网格则无关紧要 . 重要的是你有1000点 x_i ，你需要从模型中计算 f(x_i) ~ y_i .

    应该进一步注意的唯一微妙之处在于，在矢量值函数的情况下，误差的计算并不是微不足道的 . 在我看来，使用向量值函数的2范数在每个 x_i 点定义误差是很好的 . 但是嘿：在这种情况下， x_i 处的平方误差是

    |f(x_i; a,b)-y_i|^2 == sum_k (f(x_i; a,b)[k]-y_i[k])^2

    这意味着累积了每个组件的平方误差 . 这只是意味着你正在做的事情恰到好处：通过复制你的 x_i 点并单独考虑函数的每个组成部分，你的平方误差将包含每个点的误差的2范数 .

    所以我的观点是你所做的是在数学上是正确的，我不认为拟合程序的任何行为取决于多变量/向量值函数的处理方式 .

    回复于
  • 1

    如果我可以如此直率地推荐我自己的包 symfit ，我认为它确实是你所需要的 . 有关拟合共享参数的示例，请参见docs .

    您上面提到的具体问题将成为：

    from symfit import variables, parameters, Model, Fit, sin, exp

x, y_1, y_2, y_3 = variables('x, y_1, y_2, y_3')
a, b = parameters('a, b')
a.value = 0.3
b.value = 0.1

model = Model({
    y_1: a * sin(b * x), 
    y_2: a * x**2 - b * x, 
    y_3: a * exp(b / x),
})

xdata = np.linspace(1, 20, 50)
ydata = model(x=xdata, a=0.1, b=0.5)
y_noisy = ydata + 0.2 * np.random.normal(size=(len(model), len(xdata)))

fit = Fit(model, x=xdata, y_1=y_noisy[0], y_2=y_noisy[1], y_3=y_noisy[2])
fit_result = fit.execute()

    查看docs了解更多信息！

    回复于

相关问题