一次性对 scipy 的“curve_fit”进行多次迭代答案

【问题标题】：Doing many iterations of scipy's `curve_fit` in one go一次性对 scipy 的“curve_fit”进行多次迭代
【发布时间】：2022-01-21 16:10:53
【问题描述】：

考虑以下 MWE

import numpy as np
from scipy.optimize import curve_fit
X=np.arange(1,10,1)
Y=abs(X+np.random.randn(15,9))

def linear(x, a, b):
    return (x/b)**a

coeffs=[]
for ix in range(Y.shape[0]):
    print(ix)
    c0, pcov = curve_fit(linear, X, Y[ix])
    coeffs.append(c0)


XX=np.tile(X, Y.shape[0])
c0, pcov = curve_fit(linear, XX, Y.flatten())

我有一个问题，我基本上必须这样做，但不是 15 次迭代，而是数千次，而且速度很慢。

有没有办法用curve_fit 一次性完成所有这些迭代？我知道函数的结果应该是一维数组，所以只需像这样传递参数

c0, pcov = curve_fit(nlinear, X, Y)

不会工作。另外我认为答案必须是扁平化Y，所以我可以得到扁平化的结果，但我无法得到任何工作。

编辑

我知道如果我做类似的事情

XX=np.tile(X, Y.shape[0])
c0, pcov = curve_fit(nlinear, XX, Y.flatten())

然后我得到系数的“平均”值，但这不是我想要的。

编辑 2

为了记录，我使用 Jacques Kvam 的设置解决了问题，但使用 Numpy 实现（由于限制）

lX = np.log(X)
lY = np.log(Y)
A = np.vstack([lX, np.ones(len(lX))]).T
m, c=np.linalg.lstsq(A, lY.T)[0]

然后m 是a 并得到b：

b=np.exp(-c/m)

【问题讨论】：

函数linear 是您想要拟合的实际函数，还是更复杂函数的简化？
@WarrenWeckesser 这是实际功能。（是的，关于名字，我本来打算做一个线性例子，但在最后一刻改变了，我的错！）
好的，谢谢。还有一个迂腐的问题：在您的示例中，您通过将randn() 创建的噪声添加到X 来创建Y。这意味着Y 可能有负值。你的真实数据可能有负值吗？另外，真正的Y 可以包含 0 吗？
我撒了谎——还有一个问题（但它确实是上一个问题的一部分）：在你的真实数据中，X 是否包含 0？
可以接受对数并使用线性最小二乘吗？

标签： python numpy scipy

【解决方案1】：

最小二乘法不会给出相同的结果，因为在这种情况下噪声是由 log 转换的。如果噪声为零，则两种方法的结果相同。

import numpy as np
from numpy import random as rng
from scipy.optimize import curve_fit
rng.seed(0)
X=np.arange(1,7)
Y = np.zeros((4, 6))
for i in range(4):
    b = a = i + 1
    Y[i] = (X/b)**a + 0.01 * randn(6)

def linear(x, a, b):
    return (x/b)**a

coeffs=[]
for ix in range(Y.shape[0]):
    print(ix)
    c0, pcov = curve_fit(linear, X, Y[ix])
    coeffs.append(c0)

coefs是

[array([ 0.99309127,  0.98742861]),
 array([ 2.00197613,  2.00082722]),
 array([ 2.99130237,  2.99390585]),
 array([ 3.99644048,  3.9992937 ])]

我将使用 scikit-learn 的线性回归实现，因为我相信它可以很好地扩展。

from sklearn.linear_model import LinearRegression
lr = LinearRegression()

记录X和Y

lX = np.log(X)[None, :]
lY = np.log(Y)

现在拟合并检查系数是否与以前相同。

lr.fit(lX.T, lY.T)
lr.coef_

给出相似的指数。

array([ 0.98613517,  1.98643974,  2.96602892,  4.01718514])

现在检查除数。

np.exp(-lr.intercept_ / lr.coef_.ravel())

它给出了相似的系数，你可以看到这些方法在他们的答案中有所不同。

array([ 0.99199406,  1.98234916,  2.90677142,  3.73416501])

【讨论】：

您也可以使用scipy.linalg.lstsq 或numpy.linalg.lstsq。
加回噪声，并将使用日志计算的结果与curve_fit 产生的结果进行比较。在某些情况下，您会发现它们完全不同。使用线性回归拟合对数时，x 值较低的数据比 x 值较高的数据更重要。
[抱歉给你发垃圾邮件！] 我并不是说对 x 和 y 的对数使用线性最小二乘法是错误的，但通常不会给出与 @987654336 相同的结果@.
是的，我消除了噪音以获得相同的答案:) 并且因为他说可以采用日志并使用最小二乘法。不过最好指出这一点。
@WarrenWeckesser 和 JacquesKvam 我已经使用这个解决方案几个星期了，效果很好！但是，出现了类似的问题，该解决方案不起作用。您可能也有兴趣查看这个新问题，您可以查看here！在这个新案例中，我试图拟合的函数是不连续的，所以它变得更难了。干杯

【解决方案2】：

在某些情况下，将最适合的参数作为 numpy 数组用于进一步计算可能很有用。可以在 for 循环之后添加以下内容：

bestfit_par = np.asarray(coeffs)

【讨论】：