在 Scipy 中用 curve_fit 拟合向量函数

Fitting a vector function with curve_fit in Scipy

我想使用 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 的这种使用被破坏并导致运行时间过长,我想知道我应该怎么做改为做。

我认为从效率的角度来看,您所做的非常好。我将尝试查看实现并提出更多量化的东西,但暂时这是我的推理。

你在曲线拟合期间所做的是优化参数 (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') 对来产生更小的误差。我的观点是,上面的错误实际上是拟合算法的黑匣子:拟合的配置 space 仅由 (a,b) 参数定义。如果您想象如何实现一个简单的曲线拟合函数,您可以想象您尝试做梯度下降,并将平方误差作为成本函数。

现在黑盒如何计算误差应该与拟合过程无关。很容易看出 x_i 的维数与标量函数无关,因为无论您有 1000 个 1d 查询点还是 3d space 中的 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 范数点.

所以我的观点是您所做的在数学上是正确的,我不希望拟合过程的任何行为取决于 multivariate/vector-valued 函数的处理方式。

如果我能直截了当地推荐我自己的软件包 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 了解更多信息!